wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-02 16:48:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Cleaned up projection code and put it in LOM.

Browse Source
This commit is contained in:
Rowan Cockett
2013-08-05 17:16:01 -07:00
parent 83fe9df743
commit 4b545fb1b7
2 changed files with 30 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/LogicallyOrthogonalMesh.py
+12 -1
View File
@@ -207,7 +207,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
volTetra(self.gridN, A, E, F, H) + # cutted edge top
volTetra(self.gridN, A, H, F, C) + # middle
volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
self._vol = (vol1 + vol2)/2
return self._vol
@@ -319,6 +319,17 @@ NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M')
_tangents = None
tangents = property(**tangents())
def projectFaceVector(self, fV):
"""Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals"""
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)'
return mkvc(np.sum(fV*self.normals, 1), 2)
def projectEdgeVector(self, eV):
"""Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents"""
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)'
return mkvc(np.sum(eV*self.tangents, 1), 2)
if __name__ == '__main__':
nc = 5
SimPEG/tests/test_massMatrices.py
+18 -62
View File
@@ -70,26 +70,11 @@ class TestInnerProducts(OrderTest):
Ez = call(ez, self.M.gridEz)
E = np.matrix(np.r_[Ex, Ey, Ez]).T
elif self.M._meshType == 'LOM':
Tx = self.M.r(self.M.tangents, 'E', 'Ex', 'V')
Ty = self.M.r(self.M.tangents, 'E', 'Ey', 'V')
Tz = self.M.r(self.M.tangents, 'E', 'Ez', 'V')
EX_x = call(ex, self.M.gridEx)
EY_x = call(ey, self.M.gridEx)
EZ_x = call(ez, self.M.gridEx)
Ex = np.sum(np.c_[EX_x, EY_x, EZ_x]*np.c_[Tx[0], Tx[1], Tx[2]], 1)
EX_y = call(ex, self.M.gridEy)
EY_y = call(ey, self.M.gridEy)
EZ_y = call(ez, self.M.gridEy)
Ey = np.sum(np.c_[EX_y, EY_y, EZ_y]*np.c_[Ty[0], Ty[1], Ty[2]], 1)
EX_z = call(ex, self.M.gridEz)
EY_z = call(ey, self.M.gridEz)
EZ_z = call(ez, self.M.gridEz)
Ez = np.sum(np.c_[EX_z, EY_z, EZ_z]*np.c_[Tz[0], Tz[1], Tz[2]], 1)
E = np.matrix(np.r_[Ex, Ey, Ez]).T
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy),
cart(self.M.gridEz)))
E = np.matrix(self.M.projectEdgeVector(Ec))
A = self.M.getEdgeInnerProduct(sigma)
numeric = E.T*A*E
elif self.location == 'faces':
@@ -99,26 +84,11 @@ class TestInnerProducts(OrderTest):
Fz = call(ez, self.M.gridFz)
F = np.matrix(np.r_[Fx, Fy, Fz]).T
elif self.M._meshType == 'LOM':
Nx = self.M.r(self.M.normals, 'F', 'Fx', 'V')
Ny = self.M.r(self.M.normals, 'F', 'Fy', 'V')
Nz = self.M.r(self.M.normals, 'F', 'Fz', 'V')
FX_x = call(ex, self.M.gridFx)
FY_x = call(ey, self.M.gridFx)
FZ_x = call(ez, self.M.gridFx)
Fx = np.sum(np.c_[FX_x, FY_x, FZ_x]*np.c_[Nx[0], Nx[1], Nx[2]], 1)
FX_y = call(ex, self.M.gridFy)
FY_y = call(ey, self.M.gridFy)
FZ_y = call(ez, self.M.gridFy)
Fy = np.sum(np.c_[FX_y, FY_y, FZ_y]*np.c_[Ny[0], Ny[1], Ny[2]], 1)
FX_z = call(ex, self.M.gridFz)
FY_z = call(ey, self.M.gridFz)
FZ_z = call(ez, self.M.gridFz)
Fz = np.sum(np.c_[FX_z, FY_z, FZ_z]*np.c_[Nz[0], Nz[1], Nz[2]], 1)
F = np.matrix(np.r_[Fx, Fy, Fz]).T
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy),
cart(self.M.gridFz)))
F = np.matrix(self.M.projectFaceVector(Fc))
A = self.M.getFaceInnerProduct(sigma)
numeric = F.T*A*F
@@ -199,18 +169,11 @@ class TestInnerProducts2D(OrderTest):
Ey = call(ey, self.M.gridEy)
E = np.matrix(np.r_[Ex, Ey]).T
elif self.M._meshType == 'LOM':
Tx = self.M.r(self.M.tangents, 'E', 'Ex', 'V')
Ty = self.M.r(self.M.tangents, 'E', 'Ey', 'V')
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy)))
E = np.matrix(self.M.projectEdgeVector(Ec))
EX_x = call(ex, self.M.gridEx)
EY_x = call(ey, self.M.gridEx)
Ex = np.sum(np.c_[EX_x, EY_x]*np.c_[Tx[0], Tx[1]], 1)
EX_y = call(ex, self.M.gridEy)
EY_y = call(ey, self.M.gridEy)
Ey = np.sum(np.c_[EX_y, EY_y]*np.c_[Ty[0], Ty[1]], 1)
E = np.matrix(np.r_[Ex, Ey]).T
A = self.M.getEdgeInnerProduct(sigma)
numeric = E.T*A*E
elif self.location == 'faces':
@@ -219,18 +182,11 @@ class TestInnerProducts2D(OrderTest):
Fy = call(ey, self.M.gridFy)
F = np.matrix(np.r_[Fx, Fy]).T
elif self.M._meshType == 'LOM':
Nx = self.M.r(self.M.normals, 'F', 'Fx', 'V')
Ny = self.M.r(self.M.normals, 'F', 'Fy', 'V')
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy)))
F = np.matrix(self.M.projectFaceVector(Fc))
FX_x = call(ex, self.M.gridFx)
FY_x = call(ey, self.M.gridFx)
Fx = np.sum(np.c_[FX_x, FY_x]*np.c_[Nx[0], Nx[1]], 1)
FX_y = call(ex, self.M.gridFy)
FY_y = call(ey, self.M.gridFy)
Fy = np.sum(np.c_[FX_y, FY_y]*np.c_[Ny[0], Ny[1]], 1)
F = np.matrix(np.r_[Fx, Fy]).T
A = self.M.getFaceInnerProduct(sigma)
numeric = F.T*A*F