wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-04 21:22:10 +08:00
Code Issues Packages Projects Releases Wiki Activity

FaceInnerProduct for 1D

Browse Source
This commit is contained in:
rowanc1
2014-02-16 18:31:30 -08:00
parent 5fac160c38
commit 530de6cb5c
4 changed files with 60 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/BaseMesh.py
+2 -2
View File
@@ -243,7 +243,7 @@ class BaseMesh(object):
:return: projected face vector
"""
assert type(fV) == np.ndarray, 'fV must be an ndarray'
assert len(fV.shape) == 2 and fV.shape[0] == np.sum(self.nF) and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
return np.sum(fV*self.normals, 1)
def projectEdgeVector(self, eV):
@@ -255,7 +255,7 @@ class BaseMesh(object):
:return: projected edge vector
"""
assert type(eV) == np.ndarray, 'eV must be an ndarray'
assert len(eV.shape) == 2 and eV.shape[0] == np.sum(self.nE) and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
return np.sum(eV*self.tangents, 1)
SimPEG/Mesh/InnerProducts.py
+55 -10
View File
@@ -142,7 +142,14 @@ class InnerProducts(object):
Note that this is completed for each cell in the mesh at the same time.
"""
if M.dim == 2:
if M.dim == 1:
v = np.sqrt(0.5*M.vol)
V1 = sdiag(v) # We will multiply on each side to keep symmetry
Px = _getFacePx(M)
P000 = V2*Px('fXm')
P100 = V2*Px('fXp')
elif M.dim == 2:
# Square root of cell volume multiplied by 1/4
v = np.sqrt(0.25*M.vol)
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
@@ -168,9 +175,13 @@ class InnerProducts(object):
P111 = V3*Pxxx('fXp', 'fYp', 'fZp')
Mu = _makeTensor(M, mu)
A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
P = [P000, P100, P010, P110]
if M.dim == 3:
A = P000.T*Mu*P000 + P100.T*Mu*P100
P = [P000, P100]
if M.dim > 1:
A = A + P010.T*Mu*P010 + P110.T*Mu*P110
P += [P010, P110]
if M.dim > 2:
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
P += [P001, P101, P011, P111]
if returnP:
@@ -246,8 +257,10 @@ class InnerProducts(object):
Note that this is completed for each cell in the mesh at the same time.
"""
if M.dim == 1:
raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
# We will multiply by V on each side to keep symmetry
if M.dim == 2:
elif M.dim == 2:
# Square root of cell volume multiplied by 1/4
v = np.sqrt(0.25*M.vol)
V = sdiag(np.r_[v, v])
@@ -304,7 +317,13 @@ def _makeTensor(M, sigma):
if sigma is None: # default is ones
sigma = np.ones((M.nC, 1))
if M.dim == 2:
if M.dim == 1:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(sigma)
else:
raise Exception('Unexpected shape of sigma')
elif M.dim == 2:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma])
@@ -314,6 +333,8 @@ def _makeTensor(M, sigma):
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
Sigma = sp.vstack((row1, row2))
else:
raise Exception('Unexpected shape of sigma')
elif M.dim == 3:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
@@ -325,8 +346,15 @@ def _makeTensor(M, sigma):
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
Sigma = sp.vstack((row1, row2, row3))
else:
raise Exception('Unexpected shape of sigma')
return Sigma
def _getFacePx(M):
assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
return _getFacePx_Rectangular(M)
def _getFacePxx(M):
if M._meshType == 'TREE':
return M._getFacePxx
@@ -351,6 +379,23 @@ def _getEdgePxxx(M):
return _getEdgePxxx_Rectangular(M)
def _getFacePx_Rectangular(M):
"""Returns a function for creating projection matrices
"""
i = np.int64(range(M.nCx))
def Px(xFace):
"""
xFace is 'fXp' or 'fXm'
"""
posFx = 0 if xFace == 'fXm' else 1
IND = ii + posFx
PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF))
return PX
return Px
def _getFacePxx_Rectangular(M):
"""returns a function for creating projection matrices
@@ -408,7 +453,7 @@ def _getFacePxx_Rectangular(M):
IND = np.r_[ind1, ind2].flatten()
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nF)))
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
if M._meshType == 'LOM':
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
@@ -465,7 +510,7 @@ def _getFacePxxx_Rectangular(M):
IND = np.r_[ind1, ind2, ind3].flatten()
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nF))).tocsr()
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
if M._meshType == 'LOM':
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
@@ -500,7 +545,7 @@ def _getEdgePxx_Rectangular(M):
IND = np.r_[ind1, ind2].flatten()
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nE))).tocsr()
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
if M._meshType == 'LOM':
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
@@ -543,7 +588,7 @@ def _getEdgePxxx_Rectangular(M):
IND = np.r_[ind1, ind2, ind3].flatten()
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nE))).tocsr()
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
if M._meshType == 'LOM':
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
SimPEG/Mesh/LogicallyOrthogonalMesh.py
+1
View File
@@ -32,6 +32,7 @@ class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
def __init__(self, nodes):
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
assert len(nodes) > 1, "len(node) must be greater than 1"
for i, nodes_i in enumerate(nodes):
assert type(nodes_i) == np.ndarray, ("nodes[%i] is not a numpy array." % i)
SimPEG/Mesh/TreeMesh.py
+2
View File
@@ -679,6 +679,8 @@ class TreeMesh(InnerProducts, BaseMesh):
def __init__(self, h_in, x0=None):
assert type(h_in) is list, 'h_in must be a list'
assert len(h_in) > 1, "len(h_in) must be greater than 1"
h = range(len(h_in))
for i, h_i in enumerate(h_in):
if type(h_i) in [int, long, float]: