diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index 1c29c947..71788230 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -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 diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index bfcc80c5..07d67e9c 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -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