diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index 5594b5b9..ae11647a 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -3,12 +3,20 @@ from BaseMesh import BaseMesh from DiffOperators import DiffOperators from utils import mkvc, ndgrid, volTetra, indexCube, faceInfo +# Some helper functions. +length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 +length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 +normalize2D = lambda x: x/np.kron(np.ones((1, 2)), mkvc(length2D(x), 2)) +normalize3D = lambda x: x/np.kron(np.ones((1, 3)), mkvc(length3D(x), 2)) + class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid """ LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes. """ + _meshType = 'LOM' + def __init__(self, nodes): assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray" @@ -117,12 +125,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid # Compute areas of cell faces if(self.dim == 2): xy = self.gridN - length = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 - A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) - area1 = length(xy[B, :] - xy[A, :]) + area1 = length2D(xy[B, :] - xy[A, :]) A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) - area2 = length(xy[D, :] - xy[A, :]) + area2 = length2D(xy[D, :] - xy[A, :]) self._area = np.r_[mkvc(area1), mkvc(area2)] elif(self.dim == 3): @@ -141,6 +147,45 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid _area = None area = property(**area()) + def edge(): + doc = "Edge legnths." + + def fget(self): + if(self._edge is None or self._tangents is None): + if(self.dim == 2): + xy = self.gridN + A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy])) + edge1 = xy[D, :] - xy[A, :] + A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy])) + edge2 = xy[B, :] - xy[A, :] + self._edge = np.r_[mkvc(length2D(edge1)), mkvc(length2D(edge2))] + self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge] + elif(self.dim == 3): + xyz = self.gridN + A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz])) + edge1 = xyz[D, :] - xyz[A, :] + A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz])) + edge2 = xyz[B, :] - xyz[A, :] + A, E = indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz])) + edge3 = xyz[E, :] - xyz[A, :] + self._edge = np.r_[mkvc(length3D(edge1)), mkvc(length3D(edge2)), mkvc(length3D(edge3))] + self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge] + return self._edge + return locals() + _edge = None + edge = property(**edge()) + + def tangents(): + doc = "Edge tangents." + + def fget(self): + if(self._tangents is None): + self.edge # calling .edge will create the tangents + return self._tangents + return locals() + _tangents = None + tangents = property(**tangents()) + if __name__ == '__main__': nc = 5 @@ -148,7 +193,7 @@ if __name__ == '__main__': nc = 7 h2 = np.cumsum(np.r_[0, np.ones(nc)/(nc)]) h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)]) - dee3 = True + dee3 = False if dee3: X, Y, Z = ndgrid(h1, h2, h3, vector=False) M = LogicallyOrthogonalMesh([X, Y, Z]) @@ -159,4 +204,4 @@ if __name__ == '__main__': # print M.r(M.gridCC, format='M') # print M.gridN[:, 0] print M.nE - print M.area + print M.r(M.tangents, 'E', 'Ex', 'M') diff --git a/SimPEG/tests/test_LogicallyOrthogonalMesh.py b/SimPEG/tests/test_LogicallyOrthogonalMesh.py new file mode 100644 index 00000000..6c6eade9 --- /dev/null +++ b/SimPEG/tests/test_LogicallyOrthogonalMesh.py @@ -0,0 +1,52 @@ +import numpy as np +import unittest +import sys +sys.path.append('../') +from TensorMesh import TensorMesh +from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh +from OrderTest import OrderTest +from scipy.sparse.linalg import dsolve +from utils import ndgrid + + +class BasicLOMTests(unittest.TestCase): + + def setUp(self): + a = np.array([1, 1, 1]) + b = np.array([1, 2]) + c = np.array([1, 4]) + gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h] + X, Y = ndgrid(gridIt([a, b]), vector=False) + self.TM2 = TensorMesh([a, b]) + self.LOM2 = LogicallyOrthogonalMesh([X, Y]) + X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) + self.TM3 = TensorMesh([a, b, c]) + self.LOM3 = LogicallyOrthogonalMesh([X, Y, Z]) + + 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]) + self.assertTrue(np.all(self.LOM3.area == test_area)) + + def test_vol_3D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) + np.testing.assert_almost_equal(self.LOM3.vol, test_vol) + self.assertTrue(True) # Pass if you get past the assertion. + + def test_vol_2D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2]) + t1 = np.all(self.LOM2.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.LOM3.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.LOM2.edge == test_edge) + self.assertTrue(t1) + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/utils.py b/SimPEG/utils.py index c4280457..55771729 100644 --- a/SimPEG/utils.py +++ b/SimPEG/utils.py @@ -259,7 +259,7 @@ def faceInfo(xyz, A, B, C, D, average=True): nD = cross(DA, CD) length = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 - normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(N), 2)) + normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(x), 2)) if average: # average the normals at each vertex. N = (nA + nB + nC + nD)/4 # this is intrinsically weighted by area