mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 17:25:03 +08:00
Merged in diffOperators (pull request #4)
Differential Operators are now working.
This commit is contained in:
Lars Ruthotto
@@ -0,0 +1,260 @@
|
||||
import numpy as np
|
||||
from scipy import sparse as sp
|
||||
from sputils import sdiag, speye, kron3, spzeros
|
||||
from utils import mkvc
|
||||
|
||||
|
||||
def ddx(n):
|
||||
"""Define 1D derivatives, inner, this means we go from n+1 to n+1"""
|
||||
return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
|
||||
|
||||
|
||||
def checkBC(bc):
|
||||
""" Checks if boundary condition 'bc' is valid. """
|
||||
if(type(bc) is str):
|
||||
bc = [bc, bc]
|
||||
assert type(bc) is list, 'bc must be a list'
|
||||
assert len(bc) == 2, 'bc must have two elements'
|
||||
|
||||
for bc_i in bc:
|
||||
assert type(bc_i) is str, "each bc must be a string"
|
||||
assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'"
|
||||
return bc
|
||||
|
||||
|
||||
def ddxCellGrad(n, bc):
|
||||
"""Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1"""
|
||||
bc = checkBC(bc)
|
||||
|
||||
D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
|
||||
# Set the first side
|
||||
if(bc[0] == 'dirichlet'):
|
||||
D[0, 0] = 2
|
||||
elif(bc[0] == 'neumann'):
|
||||
D[0, 0] = 0
|
||||
# Set the second side
|
||||
if(bc[1] == 'dirichlet'):
|
||||
D[-1, -1] = -2
|
||||
elif(bc[1] == 'neumann'):
|
||||
D[-1, -1] = 0
|
||||
return D
|
||||
|
||||
|
||||
def av(n):
|
||||
"""Define 1D averaging operator from cell-centres to nodes."""
|
||||
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
|
||||
|
||||
|
||||
class DiffOperators(object):
|
||||
"""
|
||||
Class creates the differential operators that you need!
|
||||
"""
|
||||
def __init__(self):
|
||||
raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
|
||||
|
||||
def faceDiv():
|
||||
doc = "Construct divergence operator (face-stg to cell-centres)."
|
||||
|
||||
def fget(self):
|
||||
if(self._faceDiv is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
# Compute faceDivergence operator on faces
|
||||
if(self.dim == 1):
|
||||
D = ddx(n[0])
|
||||
elif(self.dim == 2):
|
||||
D1 = sp.kron(speye(n[1]), ddx(n[0]))
|
||||
D2 = sp.kron(ddx(n[1]), speye(n[0]))
|
||||
D = sp.hstack((D1, D2), format="csr")
|
||||
elif(self.dim == 3):
|
||||
D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
|
||||
D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
|
||||
D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
|
||||
D = sp.hstack((D1, D2, D3), format="csr")
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.area
|
||||
V = self.vol
|
||||
self._faceDiv = sdiag(1/V)*D*sdiag(S)
|
||||
|
||||
return self._faceDiv
|
||||
return locals()
|
||||
_faceDiv = None
|
||||
faceDiv = property(**faceDiv())
|
||||
|
||||
def nodalGrad():
|
||||
doc = "Construct gradient operator (nodes to edges)."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalGrad is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
# Compute divergence operator on faces
|
||||
if(self.dim == 1):
|
||||
G = ddx(n[0])
|
||||
elif(self.dim == 2):
|
||||
D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
|
||||
D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
|
||||
G = sp.vstack((D1, D2), format="csr")
|
||||
elif(self.dim == 3):
|
||||
D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
|
||||
D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
|
||||
D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
|
||||
G = sp.vstack((D1, D2, D3), format="csr")
|
||||
# Compute lengths of cell edges
|
||||
L = self.edge
|
||||
self._nodalGrad = sdiag(1/L)*G
|
||||
return self._nodalGrad
|
||||
return locals()
|
||||
_nodalGrad = None
|
||||
nodalGrad = property(**nodalGrad())
|
||||
|
||||
def setCellGradBC(self, BC):
|
||||
"""
|
||||
Function that sets the boundary conditions for cell-centred derivative operators.
|
||||
|
||||
Examples:
|
||||
|
||||
BC = 'neumann' # Neumann in all directions
|
||||
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
|
||||
# Dirichlet else
|
||||
|
||||
"""
|
||||
if(type(BC) is str):
|
||||
BC = [BC for _ in self.n] # Repeat the str self.dim times
|
||||
elif(type(BC) is list):
|
||||
assert len(BC) == self.dim, 'BC list must be the size of your mesh'
|
||||
else:
|
||||
raise Exception("BC must be a str or a list.")
|
||||
|
||||
for i, bc_i in enumerate(BC):
|
||||
BC[i] = checkBC(bc_i)
|
||||
|
||||
self._cellGrad = None # ensure we create a new gradient next time we call it
|
||||
self._cellGradBC = BC
|
||||
return BC
|
||||
_cellGradBC = 'neumann'
|
||||
|
||||
def cellGrad():
|
||||
doc = "The cell centered Gradient, takes you to cell faces."
|
||||
|
||||
def fget(self):
|
||||
if(self._cellGrad is None):
|
||||
BC = self.setCellGradBC(self._cellGradBC)
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
G = ddxCellGrad(n[0], BC[0])
|
||||
elif(self.dim == 2):
|
||||
G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
|
||||
G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
|
||||
G = sp.vstack((G1, G2), format="csr")
|
||||
elif(self.dim == 3):
|
||||
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
|
||||
G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
|
||||
G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
G = sp.vstack((G1, G2, G3), format="csr")
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.area
|
||||
V = self.vol
|
||||
self._cellGrad = sdiag(S)*G*sdiag(1/V)
|
||||
return self._cellGrad
|
||||
return locals()
|
||||
_cellGrad = None
|
||||
cellGrad = property(**cellGrad())
|
||||
|
||||
def edgeCurl():
|
||||
doc = "Construct the 3D curl operator."
|
||||
|
||||
def fget(self):
|
||||
if(self._edgeCurl is None):
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(self.hx)
|
||||
n2 = np.size(self.hy)
|
||||
n3 = np.size(self.hz)
|
||||
|
||||
# Compute lengths of cell edges
|
||||
L = self.edge
|
||||
|
||||
# Compute areas of cell faces
|
||||
S = self.area
|
||||
|
||||
# Compute divergence operator on faces
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
|
||||
D32 = kron3(d3, speye(n2), speye(n1+1))
|
||||
D23 = kron3(speye(n3), d2, speye(n1+1))
|
||||
D31 = kron3(d3, speye(n2+1), speye(n1))
|
||||
D13 = kron3(speye(n3), speye(n2+1), d1)
|
||||
D21 = kron3(speye(n3+1), d2, speye(n1))
|
||||
D12 = kron3(speye(n3+1), speye(n2), d1)
|
||||
|
||||
O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
|
||||
O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
|
||||
O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
|
||||
|
||||
C = sp.vstack((sp.hstack((O1, -D32, D23)),
|
||||
sp.hstack((D31, O2, -D13)),
|
||||
sp.hstack((-D21, D12, O3))), format="csr")
|
||||
|
||||
self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
|
||||
return self._edgeCurl
|
||||
return locals()
|
||||
_edgeCurl = None
|
||||
edgeCurl = property(**edgeCurl())
|
||||
|
||||
def faceAve():
|
||||
doc = "Construct the averaging operator on cell faces to cell centers."
|
||||
|
||||
def fget(self):
|
||||
if(self._faceAve is None):
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
self._faceAve = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._faceAve = sp.hstack((sp.kron(speye(n[1]), av(n[0])),
|
||||
sp.kron(av(n[1]), speye(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._faceAve = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._faceAve
|
||||
return locals()
|
||||
_faceAve = None
|
||||
faceAve = property(**faceAve())
|
||||
|
||||
def edgeAve():
|
||||
doc = "Construct the averaging operator on cell edges to cell centers."
|
||||
|
||||
def fget(self):
|
||||
if(self._edgeAve is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
elif(self.dim == 2):
|
||||
self._edgeAve = sp.hstack((sp.kron(av(n[1]), speye(n[0])),
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._edgeAve = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._edgeAve
|
||||
return locals()
|
||||
_edgeAve = None
|
||||
edgeAve = property(**edgeAve())
|
||||
|
||||
def getEdgeMass(self, materialProp=None):
|
||||
"""mass matix for products of edge functions w'*M(materialProp)*e"""
|
||||
if(materialProp is None):
|
||||
materialProp = np.ones(self.nC)
|
||||
Av = self.edgeAve
|
||||
return sdiag(Av.T * (self.vol * mkvc(materialProp)))
|
||||
|
||||
def getFaceMass(self, materialProp=None):
|
||||
"""mass matix for products of edge functions w'*M(materialProp)*e"""
|
||||
if(materialProp is None):
|
||||
materialProp = np.ones(self.nC)
|
||||
Av = self.faceAve
|
||||
return sdiag(Av.T*(self.vol*mkvc(materialProp)))
|
||||
+79
-4
@@ -1,10 +1,11 @@
|
||||
import numpy as np
|
||||
from BaseMesh import BaseMesh
|
||||
from TensorView import TensorView
|
||||
from utils import ndgrid
|
||||
from DiffOperators import DiffOperators
|
||||
from utils import ndgrid, mkvc
|
||||
|
||||
|
||||
class TensorMesh(BaseMesh, TensorView):
|
||||
class TensorMesh(BaseMesh, TensorView, DiffOperators):
|
||||
"""
|
||||
TensorMesh is a mesh class that deals with tensor product meshes.
|
||||
|
||||
@@ -21,15 +22,16 @@ class TensorMesh(BaseMesh, TensorView):
|
||||
|
||||
"""
|
||||
def __init__(self, h, x0=None):
|
||||
super(TensorMesh, self).__init__(np.array([len(x) for x in h]), x0)
|
||||
super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0)
|
||||
|
||||
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
|
||||
|
||||
for i, h_i in enumerate(h):
|
||||
assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
|
||||
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
|
||||
|
||||
# Ensure h contains 1D vectors
|
||||
self._h = [x.ravel() for x in h]
|
||||
self._h = [mkvc(x) for x in h]
|
||||
|
||||
def h():
|
||||
doc = "h is a list containing the cell widths of the tensor mesh in each dimension."
|
||||
@@ -186,6 +188,79 @@ class TensorMesh(BaseMesh, TensorView):
|
||||
def getCellNumbering(self):
|
||||
pass
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
def vol():
|
||||
doc = "Construct cell volumes of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._vol is None):
|
||||
vh = self.h
|
||||
# Compute cell volumes
|
||||
if(self.dim == 1):
|
||||
self._vol = mkvc(vh[0])
|
||||
elif(self.dim == 2):
|
||||
# Cell sizes in each direction
|
||||
self._vol = mkvc(np.outer(vh[0], vh[1]))
|
||||
elif(self.dim == 3):
|
||||
# Cell sizes in each direction
|
||||
self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2]))
|
||||
return self._vol
|
||||
return locals()
|
||||
_vol = None
|
||||
vol = property(**vol())
|
||||
|
||||
def area():
|
||||
doc = "Construct face areas of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._area is None):
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 1):
|
||||
self._area = np.ones(n[0]+1)
|
||||
elif(self.dim == 2):
|
||||
area1 = np.outer(np.ones(n[0]+1), vh[1])
|
||||
area2 = np.outer(vh[0], np.ones(n[1]+1))
|
||||
self._area = np.r_[mkvc(area1), mkvc(area2)]
|
||||
elif(self.dim == 3):
|
||||
area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2])))
|
||||
area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2])))
|
||||
area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1))))
|
||||
self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)]
|
||||
return self._area
|
||||
return locals()
|
||||
_area = None
|
||||
area = property(**area())
|
||||
|
||||
def edge():
|
||||
doc = "Construct edge legnths of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._edge is None):
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.n
|
||||
# Compute edge lengths
|
||||
if(self.dim == 1):
|
||||
self._edge = mkvc(vh[0])
|
||||
elif(self.dim == 2):
|
||||
l1 = np.outer(vh[0], np.ones(n[1]+1))
|
||||
l2 = np.outer(np.ones(n[0]+1), vh[1])
|
||||
self._edge = np.r_[mkvc(l1), mkvc(l2)]
|
||||
elif(self.dim == 3):
|
||||
l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
|
||||
l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1))))
|
||||
l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2])))
|
||||
self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)]
|
||||
return self._edge
|
||||
return locals()
|
||||
_edge = None
|
||||
edge = property(**edge())
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
print('Welcome to tensor mesh!')
|
||||
|
||||
@@ -1 +1,2 @@
|
||||
from TensorMesh import TensorMesh
|
||||
import utils
|
||||
|
||||
@@ -1,75 +0,0 @@
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
from utils import mkvc
|
||||
from sputils import ddx, sdiag, speye, kron3
|
||||
|
||||
def getvol(h):
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# Compute cell volumes
|
||||
v12 = h1.T*h2
|
||||
V = mkvc(v12.reshape(-1,1)*h3)
|
||||
|
||||
return V
|
||||
|
||||
def getarea(h):
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
# Compute areas of cell faces
|
||||
area1 = np.ones((n1+1,1))*mkvc(h2.T*h3)
|
||||
area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
|
||||
area3 = h1.T*mkvc(h2.T*np.ones(n3+1))
|
||||
area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
|
||||
return area
|
||||
|
||||
def getDivMatrix(h):
|
||||
"""Consturct the 3D divergence operator on Faces."""
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
|
||||
# Compute areas of cell faces
|
||||
#area1 = np.ones((n1+1,1))*mkvc(h2.T*h3)
|
||||
#area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
|
||||
#area3 = h1.T*mkvc(h2.T*np.ones(n3+1))
|
||||
#area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
|
||||
area = getarea(h)
|
||||
|
||||
S = sdiag(area)
|
||||
|
||||
# Compute cell volumes
|
||||
#v12 = h1.T*h2
|
||||
#V = mkvc(v12.reshape(-1,1)*h3)
|
||||
V = getvol(h)
|
||||
|
||||
# Compute divergence operator on faces
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
D1 = kron3(speye(n3), speye(n2), d1)
|
||||
D2 = kron3(speye(n3), d2, speye(n1))
|
||||
D3 = kron3(d3, speye(n2), speye(n1))
|
||||
|
||||
D = sparse.hstack((sparse.hstack((D1, D2)), D3))
|
||||
return sdiag(1/V)*D*S
|
||||
|
||||
@@ -1,207 +0,0 @@
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
from utils import mkvc
|
||||
from sputils import *
|
||||
#from sputils import ddx, sdiag, speye, kron3, spzeros, appendBottom3,
|
||||
|
||||
def getVol(h):
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# Compute cell volumes
|
||||
V = mkvc(np.outer(mkvc(np.outer(h1,h2)),h3))
|
||||
|
||||
return V
|
||||
|
||||
def getArea(h):
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
# Compute areas of cell faces
|
||||
area1 = mkvc(np.outer(np.ones(n1+1),np.outer(h2,h3)))
|
||||
area2 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),h3))))
|
||||
area3 = mkvc(np.outer(h1,mkvc(np.outer(h2,np.ones(n3+1)))))
|
||||
area = np.hstack((np.hstack((area1, area2)), area3))
|
||||
|
||||
return area
|
||||
|
||||
def getLength(h):
|
||||
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
|
||||
# compute the length of each edge
|
||||
Length1 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),np.ones(n3+1)))))
|
||||
Length2 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(h2,np.ones(n3+1)))))
|
||||
Length3 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(np.ones(n2+1),h3))))
|
||||
|
||||
Length = np.hstack((np.hstack((Length1, Length2)), Length3))
|
||||
|
||||
return Length
|
||||
|
||||
|
||||
def getDivMatrix(h):
|
||||
"""Consturct the 3D divergence operator on Faces."""
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]
|
||||
h2 = h[1]
|
||||
h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1)
|
||||
n2 = np.size(h2)
|
||||
n3 = np.size(h3)
|
||||
|
||||
area = getArea(h)
|
||||
S = sdiag(area)
|
||||
|
||||
# Compute cell volumes
|
||||
V = getVol(h)
|
||||
|
||||
# Compute divergence operator on faces
|
||||
d1 = ddx(n1)
|
||||
d2 = ddx(n2)
|
||||
d3 = ddx(n3)
|
||||
D1 = kron3(speye(n3), speye(n2), d1)
|
||||
D2 = kron3(speye(n3), d2, speye(n1))
|
||||
D3 = kron3(d3, speye(n2), speye(n1))
|
||||
|
||||
D = sparse.hstack((sparse.hstack((D1, D2)), D3))
|
||||
return sdiag(1/V)*D*S
|
||||
|
||||
|
||||
def getCurlMatrix(h):
|
||||
"""Edge CURL """
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]; h2 = h[1]; h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3)
|
||||
|
||||
d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3)
|
||||
# derivatives on x-edge variables
|
||||
D32 = kron3(d3, speye(n2), speye(n1+1))
|
||||
D23 = kron3(speye(n3), d2, speye(n1+1))
|
||||
D31 = kron3(d3, speye(n2+1), speye(n1))
|
||||
D13 = kron3(speye(n3), speye(n2+1), d1)
|
||||
D21 = kron3(speye(n3+1), d2, speye(n1))
|
||||
D12 = kron3(speye(n3+1), speye(n2), d1)
|
||||
|
||||
O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
|
||||
O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
|
||||
O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
|
||||
|
||||
CURL = appendBottom3(
|
||||
appendRight3(O1, -D32, D23),
|
||||
appendRight3(D31, O2, -D13),
|
||||
appendRight3(-D21, D12, O3))
|
||||
|
||||
|
||||
area = getArea(h)
|
||||
S = sdiag(1/area)
|
||||
|
||||
# Compute edge length
|
||||
lngth = getLength(h)
|
||||
L = sdiag(lngth)
|
||||
|
||||
return S*(CURL*L)
|
||||
|
||||
|
||||
def getNodalGradient(h):
|
||||
"""Nodal Gradients"""
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]; h2 = h[1]; h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3)
|
||||
|
||||
|
||||
D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1))
|
||||
D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1))
|
||||
D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1))
|
||||
|
||||
# topological gradient
|
||||
GRAD = appendBottom3(D1, D2, D3)
|
||||
|
||||
# scale for non-uniform mesh
|
||||
# Compute edge length
|
||||
lngth = getLength(h)
|
||||
L = sdiag(1/lngth)
|
||||
|
||||
return L*GRAD
|
||||
|
||||
|
||||
def getEdgeToCellAverge(h):
|
||||
|
||||
"""Average from Edge to Cell center """
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]; h2 = h[1]; h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3)
|
||||
|
||||
a1 = av(n1); a2 = av(n2); a3 = av(n3)
|
||||
# derivatives on x-edge variables
|
||||
A1 = kron3(a3, a2, speye(n1))
|
||||
A2 = kron3(a3, speye(n2), a1)
|
||||
A3 = kron3(speye(n3), a2, a1)
|
||||
|
||||
return appendRight3(A1, A2, A3)
|
||||
|
||||
|
||||
def getFaceToCellAverge(h):
|
||||
|
||||
"""Average from Edge to Cell center """
|
||||
|
||||
# Cell sizes in each direction
|
||||
h1 = h[0]; h2 = h[1]; h3 = h[2]
|
||||
|
||||
# The number of cell centers in each direction
|
||||
n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3)
|
||||
|
||||
a1 = av(n1); a2 = av(n2); a3 = av(n3)
|
||||
# derivatives on x-edge variables
|
||||
A1 = kron3(speye(n3), speye(n2), a1)
|
||||
A2 = kron3(speye(n3), a2, speye(n1))
|
||||
A3 = kron3(a3, speye(n2), speye(n1))
|
||||
|
||||
return appendRight3(A1, A2, A3)
|
||||
|
||||
|
||||
def getEdgeMassMatrix(h,sigma):
|
||||
# mass matix for products of edge functions w'*M(sigma)*e
|
||||
|
||||
Av = getEdgeToCellAverge(h)
|
||||
v = getVol(h)
|
||||
sigma = mkvc(sigma)
|
||||
|
||||
return sdiag(Av.T*(v*sigma))
|
||||
|
||||
def getFaceMassMatrix(h,sigma):
|
||||
# mass matix for products of edge functions w'*M(sigma)*e
|
||||
|
||||
Av = getFaceToCellAverge(h)
|
||||
v = getVol(h)
|
||||
sigma = mkvc(sigma)
|
||||
|
||||
return sdiag(Av.T*(v*sigma))
|
||||
+17
-27
@@ -1,35 +1,30 @@
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
from numpy import ones
|
||||
from scipy import sparse as sp
|
||||
|
||||
|
||||
def ddx(n):
|
||||
"""Define 1D derivatives"""
|
||||
return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1)
|
||||
|
||||
def sdiag(h):
|
||||
"""Sparse diagonal matrix"""
|
||||
return sparse.spdiags(h, 0, np.size(h), np.size(h))
|
||||
return sp.spdiags(h, 0, np.size(h), np.size(h), format="csr")
|
||||
|
||||
|
||||
def speye(n):
|
||||
"""Sparse identity"""
|
||||
return sparse.identity(n)
|
||||
return sp.identity(n, format="csr")
|
||||
|
||||
|
||||
def kron3(A, B, C):
|
||||
"""Two kron prods"""
|
||||
return sparse.kron(sparse.kron(A, B), C)
|
||||
|
||||
def av(n):
|
||||
"""Define 1D average"""
|
||||
return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1))
|
||||
|
||||
"""Three kron prods"""
|
||||
return sp.kron(sp.kron(A, B), C, format="csr")
|
||||
|
||||
|
||||
def spzeros(n1, n2):
|
||||
"""spzeros"""
|
||||
return sparse.coo_matrix((n1, n2))
|
||||
|
||||
return sp.coo_matrix((n1, n2)).tocsr()
|
||||
|
||||
|
||||
def appendBottom(A, B):
|
||||
"""append on bottom"""
|
||||
C = sparse.vstack((A, B))
|
||||
C = sp.vstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
@@ -43,7 +38,7 @@ def appendBottom3(A, B, C):
|
||||
|
||||
def appendRight(A, B):
|
||||
"""append on right"""
|
||||
C = sparse.hstack((A, B))
|
||||
C = sp.hstack((A, B))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
@@ -57,9 +52,9 @@ def appendRight3(A, B, C):
|
||||
|
||||
def blkDiag(A, B):
|
||||
"""blockdigonal"""
|
||||
O12 = sparse.coo_matrix((np.shape(A)[0], np.shape(B)[1]))
|
||||
O21 = sparse.coo_matrix((np.shape(B)[0], np.shape(A)[1]))
|
||||
C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B))))
|
||||
O12 = sp.coo_matrix((np.shape(A)[0], np.shape(B)[1]))
|
||||
O21 = sp.coo_matrix((np.shape(B)[0], np.shape(A)[1]))
|
||||
C = sp.vstack((sp.hstack((A, O12)), sp.hstack((O21, B))))
|
||||
C = C.tocsr()
|
||||
return C
|
||||
|
||||
@@ -69,8 +64,3 @@ def blkDiag3(A, B, C):
|
||||
ABC = blkDiag(blkDiag(A, B), C)
|
||||
ABC = ABC.tocsr()
|
||||
return ABC
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,113 @@
|
||||
import sys
|
||||
sys.path.append('../../')
|
||||
from SimPEG import TensorMesh
|
||||
import numpy as np
|
||||
import unittest
|
||||
|
||||
|
||||
class OrderTest(unittest.TestCase):
|
||||
"""
|
||||
|
||||
OrderTest is a base class for testing convergence orders with respect to mesh
|
||||
sizes of integral/differential operators.
|
||||
|
||||
Mathematical Problem:
|
||||
|
||||
Given are an operator A and its discretization A[h]. For a given test function f
|
||||
and h --> 0 we compare:
|
||||
|
||||
error(h) = \| A[h](f) - A(f) \|_{\infty}
|
||||
|
||||
Note that you can provide any norm.
|
||||
|
||||
Test is passed when estimated rate order of convergence is at least 90% of the
|
||||
estimated rate supplied by the user.
|
||||
|
||||
Minimal example for a curl operator:
|
||||
|
||||
class TestCURL(OrderTest):
|
||||
name = "Curl"
|
||||
|
||||
def getError(self):
|
||||
# For given Mesh, generate A[h], f and A(f) and return norm of error.
|
||||
|
||||
|
||||
fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x))
|
||||
sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y))
|
||||
|
||||
|
||||
Ex = fun(self.M.gridEx[:, 1])
|
||||
Ey = fun(self.M.gridEy[:, 2])
|
||||
Ez = fun(self.M.gridEz[:, 0])
|
||||
f = np.concatenate((Ex, Ey, Ez))
|
||||
|
||||
Fx = sol(self.M.gridFx[:, 2])
|
||||
Fy = sol(self.M.gridFy[:, 0])
|
||||
Fz = sol(self.M.gridFz[:, 1])
|
||||
Af = np.concatenate((Fx, Fy, Fz))
|
||||
|
||||
# Generate DIV matrix
|
||||
Ah = self.M.edgeCurl
|
||||
|
||||
curlE = Ah*E
|
||||
err = np.linalg.norm((Ah*f -Af), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
# runs the test
|
||||
self.orderTest()
|
||||
|
||||
See also: test_operatorOrder.py
|
||||
|
||||
"""
|
||||
|
||||
name = "Order Test"
|
||||
expectedOrder = 2
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
|
||||
def setupMesh(self, nc):
|
||||
"""
|
||||
For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.
|
||||
"""
|
||||
h1 = np.ones(nc)/nc
|
||||
h2 = np.ones(nc)/nc
|
||||
h3 = np.ones(nc)/nc
|
||||
h = [h1, h2, h3]
|
||||
self.M = TensorMesh(h)
|
||||
|
||||
def getError(self):
|
||||
"""For given h, generate A[h], f and A(f) and return norm of error."""
|
||||
return 1.
|
||||
|
||||
def orderTest(self):
|
||||
"""
|
||||
For number of cells specified in meshSizes setup mesh, call getError
|
||||
and prints mesh size, error, ratio between current and previous error,
|
||||
and estimated order of convergence.
|
||||
|
||||
|
||||
"""
|
||||
order = []
|
||||
err_old = 0.
|
||||
nc_old = 0.
|
||||
for ii, nc in enumerate(self.meshSizes):
|
||||
self.setupMesh(nc)
|
||||
err = self.getError()
|
||||
if ii == 0:
|
||||
print ''
|
||||
print 'Testing order of: ' + self.name
|
||||
print '_____________________________________________'
|
||||
print ' h | error | e(i-1)/e(i) | order'
|
||||
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
|
||||
print '%4i | %8.2e |' % (nc, err)
|
||||
else:
|
||||
order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc)))
|
||||
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1])
|
||||
err_old = err
|
||||
nc_old = nc
|
||||
print '---------------------------------------------'
|
||||
self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order)))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -1,45 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from TensorMesh import TensorMesh
|
||||
from getDIV import getDivMatrix, getarea, getvol
|
||||
|
||||
# Define the mesh
|
||||
err=0.
|
||||
for i in range(4):
|
||||
icount=i+1;
|
||||
nc = 2*icount;
|
||||
h1 = np.pi/nc*np.ones((1,nc))
|
||||
h2 = np.pi/nc*np.ones((1,nc))
|
||||
h3 = np.pi/nc*np.ones((1,nc))
|
||||
h = [h1, h2, h3]
|
||||
x0 = -np.pi/2*np.ones((3, 1))
|
||||
M = TensorMesh(h, x0)
|
||||
#n = M.plotGrid()
|
||||
|
||||
# Generate DIV matrix
|
||||
DIV = getDivMatrix(h)
|
||||
|
||||
#Test function
|
||||
fun = lambda x: np.sin(x)
|
||||
Fx = fun(M.gridFx[:,0])
|
||||
Fy = fun(M.gridFy[:,1])
|
||||
Fz = fun(M.gridFz[:,2])
|
||||
|
||||
F = np.concatenate((Fx,Fy,Fz))
|
||||
divF = DIV*F
|
||||
sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2])
|
||||
|
||||
area = getarea(h)
|
||||
vol = getvol(h)
|
||||
err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
|
||||
if icount == 1:
|
||||
err1 = err
|
||||
print 'h | 2 norm | error ratio'
|
||||
print '---------------------------------------'
|
||||
print '%6.4f | %8.2e |'% (h1[0,0], err)
|
||||
else:
|
||||
print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err1/err)
|
||||
|
||||
@@ -3,15 +3,18 @@ import unittest
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from TensorMesh import TensorMesh
|
||||
from OrderTest import OrderTest
|
||||
from scipy.sparse.linalg import dsolve
|
||||
|
||||
|
||||
class TestSequenceFunctions(unittest.TestCase):
|
||||
class BasicTensorMeshTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
a = np.array([1, 1, 1])
|
||||
b = np.array([1, 2])
|
||||
x0 = np.array([3, 5])
|
||||
self.mesh2 = TensorMesh([a, b], x0)
|
||||
c = np.array([1, 4])
|
||||
self.mesh2 = TensorMesh([a, b], np.array([3, 5]))
|
||||
self.mesh3 = TensorMesh([a, b, c])
|
||||
|
||||
def test_vectorN_2D(self):
|
||||
testNx = np.array([3, 4, 5, 6])
|
||||
@@ -29,6 +32,140 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
ytest = np.all(self.mesh2.vectorCCy == testNy)
|
||||
self.assertTrue(xtest and ytest)
|
||||
|
||||
def test_area_3D(self):
|
||||
test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
|
||||
t1 = np.all(self.mesh3.area == test_area)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_vol_3D(self):
|
||||
test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
|
||||
t1 = np.all(self.mesh3.vol == test_vol)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_vol_2D(self):
|
||||
test_vol = np.array([1, 1, 1, 2, 2, 2])
|
||||
t1 = np.all(self.mesh2.vol == test_vol)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_edge_3D(self):
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
|
||||
t1 = np.all(self.mesh3.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_edge_2D(self):
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
|
||||
t1 = np.all(self.mesh2.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
class TestCurl(OrderTest):
|
||||
name = "Curl"
|
||||
|
||||
def getError(self):
|
||||
fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x))
|
||||
sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y))
|
||||
|
||||
Ex = fun(self.M.gridEx[:, 1])
|
||||
Ey = fun(self.M.gridEy[:, 2])
|
||||
Ez = fun(self.M.gridEz[:, 0])
|
||||
E = np.concatenate((Ex, Ey, Ez))
|
||||
|
||||
Fx = sol(self.M.gridFx[:, 2])
|
||||
Fy = sol(self.M.gridFy[:, 0])
|
||||
Fz = sol(self.M.gridFz[:, 1])
|
||||
curlE_anal = np.concatenate((Fx, Fy, Fz))
|
||||
|
||||
# Generate DIV matrix
|
||||
CURL = self.M.edgeCurl
|
||||
|
||||
curlE = CURL*E
|
||||
err = np.linalg.norm((curlE-curlE_anal), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestFaceDiv(OrderTest):
|
||||
name = "Face Divergence"
|
||||
|
||||
def getError(self):
|
||||
DIV = self.M.faceDiv
|
||||
|
||||
#Test function
|
||||
fun = lambda x: np.sin(x)
|
||||
Fx = fun(self.M.gridFx[:, 0])
|
||||
Fy = fun(self.M.gridFy[:, 1])
|
||||
Fz = fun(self.M.gridFz[:, 2])
|
||||
|
||||
F = np.concatenate((Fx, Fy, Fz))
|
||||
divF = DIV*F
|
||||
sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2])
|
||||
|
||||
err = np.linalg.norm((divF-divF_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad(OrderTest):
|
||||
name = "Nodal Gradient"
|
||||
|
||||
def getError(self):
|
||||
GRAD = self.M.nodalGrad
|
||||
#Test function
|
||||
fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z))
|
||||
|
||||
phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2])
|
||||
gradE = GRAD*phi
|
||||
|
||||
Ex = sol(self.M.gridEx[:, 0])
|
||||
Ey = sol(self.M.gridEy[:, 1])
|
||||
Ez = sol(self.M.gridEz[:, 2])
|
||||
|
||||
gradE_anal = np.concatenate((Ex, Ey, Ez))
|
||||
err = np.linalg.norm((gradE-gradE_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestPoissonEqn(OrderTest):
|
||||
name = "Poisson Equation"
|
||||
meshSizes = [16, 20, 24]
|
||||
|
||||
def getError(self):
|
||||
# Create some functions to integrate
|
||||
fun = lambda x: np.sin(2*np.pi*x[:, 0])*np.sin(2*np.pi*x[:, 1])*np.sin(2*np.pi*x[:, 2])
|
||||
sol = lambda x: -3.*((2*np.pi)**2)*fun(x)
|
||||
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
|
||||
D = self.M.faceDiv
|
||||
G = self.M.cellGrad
|
||||
if self.forward:
|
||||
sA = sol(self.M.gridCC)
|
||||
sN = D*G*fun(self.M.gridCC)
|
||||
err = np.linalg.norm((sA - sN), np.inf)
|
||||
else:
|
||||
fA = fun(self.M.gridCC)
|
||||
fN = dsolve.spsolve(D*G, sol(self.M.gridCC))
|
||||
err = np.linalg.norm((fA - fN), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderForward(self):
|
||||
self.name = "Poisson Equation - Forward"
|
||||
self.forward = True
|
||||
self.orderTest()
|
||||
|
||||
def test_orderBackward(self):
|
||||
self.name = "Poisson Equation - Backward"
|
||||
self.forward = False
|
||||
self.orderTest()
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user