diff --git a/SimPEG/Mesh/LogicallyRectMesh.py b/SimPEG/Mesh/CurvilinearMesh.py similarity index 97% rename from SimPEG/Mesh/LogicallyRectMesh.py rename to SimPEG/Mesh/CurvilinearMesh.py index b2dc6f55..e4eeda3e 100644 --- a/SimPEG/Mesh/LogicallyRectMesh.py +++ b/SimPEG/Mesh/CurvilinearMesh.py @@ -10,9 +10,9 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)) normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) -class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): +class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): """ - LogicallyRectMesh is a mesh class that deals with logically rectangular meshes. + CurvilinearMesh is a mesh class that deals with logically rectangular meshes. Example of a logically rectangular mesh: @@ -21,13 +21,13 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): from SimPEG import Mesh, Utils X, Y = Utils.exampleLrmGrid([3,3],'rotate') - M = Mesh.LogicallyRectMesh([X, Y]) + M = Mesh.CurvilinearMesh([X, Y]) M.plotGrid(showIt=True) """ __metaclass__ = Utils.SimPEGMetaClass - _meshType = 'LRM' + _meshType = 'Curv' def __init__(self, nodes): assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray" @@ -38,7 +38,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i) assert len(nodes[0].shape) == len(nodes), "Dimension mismatch" - assert len(nodes[0].shape) > 1, "Not worth using LRM for a 1D mesh." + assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh." BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None) @@ -343,7 +343,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): from SimPEG import Mesh, Utils X, Y = Utils.exampleLrmGrid([3,3],'rotate') - M = Mesh.LogicallyRectMesh([X, Y]) + M = Mesh.CurvilinearMesh([X, Y]) M.plotGrid(showIt=True) """ @@ -435,9 +435,9 @@ if __name__ == '__main__': dee3 = True if dee3: X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False) - M = LogicallyRectMesh([X, Y, Z]) + M = CurvilinearMesh([X, Y, Z]) else: X, Y = Utils.ndgrid(h1, h2, vector=False) - M = LogicallyRectMesh([X, Y]) + M = CurvilinearMesh([X, Y]) print M.r(M.normals, 'F', 'Fx', 'V') diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 2cf48a7c..42a08702 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -328,7 +328,7 @@ class InnerProducts(object): iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if M._meshType == 'LRM': + if M._meshType == 'Curv': fN1 = M.r(M.normals, 'F', 'Fx', 'M') fN2 = M.r(M.normals, 'F', 'Fy', 'M') @@ -353,7 +353,7 @@ class InnerProducts(object): PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF)) - if M._meshType == 'LRM': + if M._meshType == 'Curv': I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]), getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy])) PXX = I2x2 * PXX @@ -376,7 +376,7 @@ class InnerProducts(object): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if M._meshType == 'LRM': + if M._meshType == 'Curv': fN1 = M.r(M.normals, 'F', 'Fx', 'M') fN2 = M.r(M.normals, 'F', 'Fy', 'M') fN3 = M.r(M.normals, 'F', 'Fz', 'M') @@ -410,7 +410,7 @@ class InnerProducts(object): PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr() - if M._meshType == 'LRM': + if M._meshType == 'Curv': I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]), getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]), getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ])) @@ -432,7 +432,7 @@ class InnerProducts(object): iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if M._meshType == 'LRM': + if M._meshType == 'Curv': eT1 = M.r(M.tangents, 'E', 'Ex', 'M') eT2 = M.r(M.tangents, 'E', 'Ey', 'M') @@ -452,7 +452,7 @@ class InnerProducts(object): PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr() - if M._meshType == 'LRM': + if M._meshType == 'Curv': I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]), getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j])) PXX = I2x2 * PXX @@ -466,7 +466,7 @@ class InnerProducts(object): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if M._meshType == 'LRM': + if M._meshType == 'Curv': eT1 = M.r(M.tangents, 'E', 'Ex', 'M') eT2 = M.r(M.tangents, 'E', 'Ey', 'M') eT3 = M.r(M.tangents, 'E', 'Ez', 'M') @@ -495,7 +495,7 @@ class InnerProducts(object): PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr() - if M._meshType == 'LRM': + if M._meshType == 'Curv': 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]]), getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]), getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k])) diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index deb11321..9aadfcba 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,5 +1,5 @@ from TensorMesh import TensorMesh from CylMesh import CylMesh -from LogicallyRectMesh import LogicallyRectMesh +from CurvilinearMesh import CurvilinearMesh from TreeMesh import TreeMesh from BaseMesh import BaseMesh diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index dbcb92b6..c6aa5bc4 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag, diagEst from SimPEG import Utils -from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh +from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh import numpy as np import scipy.sparse as sp import unittest @@ -115,7 +115,7 @@ class OrderTest(unittest.TestCase): max_h = max([np.max(hi) for hi in self.M.h]) return max_h - elif 'LRM' in self._meshType: + elif 'Curv' in self._meshType: if 'uniform' in self._meshType: kwrd = 'rect' elif 'rotate' in self._meshType: @@ -126,10 +126,10 @@ class OrderTest(unittest.TestCase): raise Exception('Lom not supported for 1D') elif self.meshDimension == 2: X, Y = Utils.exampleLrmGrid([nc, nc], kwrd) - self.M = LogicallyRectMesh([X, Y]) + self.M = CurvilinearMesh([X, Y]) elif self.meshDimension == 3: X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd) - self.M = LogicallyRectMesh([X, Y, Z]) + self.M = CurvilinearMesh([X, Y, Z]) return 1./nc def getError(self): diff --git a/SimPEG/Tests/test_CurvilinearMesh.py b/SimPEG/Tests/test_CurvilinearMesh.py new file mode 100644 index 00000000..42e3d877 --- /dev/null +++ b/SimPEG/Tests/test_CurvilinearMesh.py @@ -0,0 +1,104 @@ +import numpy as np +import unittest +from SimPEG.Mesh import TensorMesh, CurvilinearMesh +from SimPEG.Utils import ndgrid + + +class BasicCurvTests(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.Curv2 = CurvilinearMesh([X, Y]) + X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) + self.TM3 = TensorMesh([a, b, c]) + self.Curv3 = CurvilinearMesh([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.Curv3.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.Curv3.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.Curv2.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.Curv3.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.Curv2.edge == test_edge) + self.assertTrue(t1) + + def test_tangents(self): + T = self.Curv2.tangents + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy))) + + T = self.Curv3.tangents + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx))) + + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy))) + + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz))) + + def test_normals(self): + N = self.Curv2.normals + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy))) + + N = self.Curv3.normals + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx))) + + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy))) + + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz))) + + def test_grid(self): + self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC)) + self.assertTrue(np.all(self.Curv2.gridN == self.TM2.gridN)) + self.assertTrue(np.all(self.Curv2.gridFx == self.TM2.gridFx)) + self.assertTrue(np.all(self.Curv2.gridFy == self.TM2.gridFy)) + self.assertTrue(np.all(self.Curv2.gridEx == self.TM2.gridEx)) + self.assertTrue(np.all(self.Curv2.gridEy == self.TM2.gridEy)) + + self.assertTrue(np.all(self.Curv3.gridCC == self.TM3.gridCC)) + self.assertTrue(np.all(self.Curv3.gridN == self.TM3.gridN)) + self.assertTrue(np.all(self.Curv3.gridFx == self.TM3.gridFx)) + self.assertTrue(np.all(self.Curv3.gridFy == self.TM3.gridFy)) + self.assertTrue(np.all(self.Curv3.gridFz == self.TM3.gridFz)) + self.assertTrue(np.all(self.Curv3.gridEx == self.TM3.gridEx)) + self.assertTrue(np.all(self.Curv3.gridEy == self.TM3.gridEy)) + self.assertTrue(np.all(self.Curv3.gridEz == self.TM3.gridEz)) + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/Tests/test_LogicallyRectMesh.py b/SimPEG/Tests/test_LogicallyRectMesh.py deleted file mode 100644 index 5d223f68..00000000 --- a/SimPEG/Tests/test_LogicallyRectMesh.py +++ /dev/null @@ -1,104 +0,0 @@ -import numpy as np -import unittest -from SimPEG.Mesh import TensorMesh, LogicallyRectMesh -from SimPEG.Utils import ndgrid - - -class BasicLRMTests(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.LRM2 = LogicallyRectMesh([X, Y]) - X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) - self.TM3 = TensorMesh([a, b, c]) - self.LRM3 = LogicallyRectMesh([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.LRM3.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.LRM3.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.LRM2.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.LRM3.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.LRM2.edge == test_edge) - self.assertTrue(t1) - - def test_tangents(self): - T = self.LRM2.tangents - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM2.nEx))) - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM2.nEx))) - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM2.nEy))) - self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM2.nEy))) - - T = self.LRM3.tangents - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM3.nEx))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM3.nEx))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LRM3.nEx))) - - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM3.nEy))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM3.nEy))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LRM3.nEy))) - - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LRM3.nEz))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LRM3.nEz))) - self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LRM3.nEz))) - - def test_normals(self): - N = self.LRM2.normals - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM2.nFx))) - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM2.nFx))) - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM2.nFy))) - self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM2.nFy))) - - N = self.LRM3.normals - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM3.nFx))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM3.nFx))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LRM3.nFx))) - - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM3.nFy))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM3.nFy))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LRM3.nFy))) - - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LRM3.nFz))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LRM3.nFz))) - self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LRM3.nFz))) - - def test_grid(self): - self.assertTrue(np.all(self.LRM2.gridCC == self.TM2.gridCC)) - self.assertTrue(np.all(self.LRM2.gridN == self.TM2.gridN)) - self.assertTrue(np.all(self.LRM2.gridFx == self.TM2.gridFx)) - self.assertTrue(np.all(self.LRM2.gridFy == self.TM2.gridFy)) - self.assertTrue(np.all(self.LRM2.gridEx == self.TM2.gridEx)) - self.assertTrue(np.all(self.LRM2.gridEy == self.TM2.gridEy)) - - self.assertTrue(np.all(self.LRM3.gridCC == self.TM3.gridCC)) - self.assertTrue(np.all(self.LRM3.gridN == self.TM3.gridN)) - self.assertTrue(np.all(self.LRM3.gridFx == self.TM3.gridFx)) - self.assertTrue(np.all(self.LRM3.gridFy == self.TM3.gridFy)) - self.assertTrue(np.all(self.LRM3.gridFz == self.TM3.gridFz)) - self.assertTrue(np.all(self.LRM3.gridEx == self.TM3.gridEx)) - self.assertTrue(np.all(self.LRM3.gridEy == self.TM3.gridEy)) - self.assertTrue(np.all(self.LRM3.gridEz == self.TM3.gridEz)) - - -if __name__ == '__main__': - unittest.main() diff --git a/SimPEG/Tests/test_innerProduct.py b/SimPEG/Tests/test_innerProduct.py index bbaf4b14..0a5ff809 100644 --- a/SimPEG/Tests/test_innerProduct.py +++ b/SimPEG/Tests/test_innerProduct.py @@ -7,7 +7,7 @@ from SimPEG import Utils class TestInnerProducts(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] + meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv'] meshDimension = 3 meshSizes = [16, 32] @@ -154,7 +154,7 @@ class TestInnerProducts(OrderTest): class TestInnerProducts2D(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] + meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv'] meshDimension = 2 meshSizes = [4, 8, 16, 32, 64, 128] diff --git a/SimPEG/Tests/test_innerProductDerivs.py b/SimPEG/Tests/test_innerProductDerivs.py index 4b7d63cd..6eb561c1 100644 --- a/SimPEG/Tests/test_innerProductDerivs.py +++ b/SimPEG/Tests/test_innerProductDerivs.py @@ -7,9 +7,9 @@ from TestUtils import checkDerivative class TestInnerProductsDerivs(unittest.TestCase): def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False): - if meshType == 'LRM': + if meshType == 'Curv': hRect = Utils.exampleLrmGrid(h,'rotate') - mesh = Mesh.LogicallyRectMesh(hRect) + mesh = Mesh.CurvilinearMesh(hRect) elif meshType == 'Tree': mesh = Mesh.TreeMesh(h) elif meshType == 'Tensor': @@ -24,9 +24,9 @@ class TestInnerProductsDerivs(unittest.TestCase): return checkDerivative(fun, sig, num=5, plotIt=False) def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False): - if meshType == 'LRM': + if meshType == 'Curv': hRect = Utils.exampleLrmGrid(h,'rotate') - mesh = Mesh.LogicallyRectMesh(hRect) + mesh = Mesh.CurvilinearMesh(hRect) elif meshType == 'Tree': mesh = Mesh.TreeMesh(h) elif meshType == 'Tensor': @@ -137,65 +137,65 @@ class TestInnerProductsDerivs(unittest.TestCase): - def test_FaceIP_2D_float_LRM(self): - self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM')) - def test_FaceIP_3D_float_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM')) - def test_FaceIP_2D_isotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM')) - def test_FaceIP_3D_isotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM')) - def test_FaceIP_2D_anisotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM')) - def test_FaceIP_3D_anisotropic_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM')) - def test_FaceIP_2D_tensor_LRM(self): - self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM')) - def test_FaceIP_3D_tensor_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM')) + def test_FaceIP_2D_float_Curv(self): + self.assertTrue(self.doTestFace([10, 4],0, False, 'Curv')) + def test_FaceIP_3D_float_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Curv')) + def test_FaceIP_2D_isotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4],1, False, 'Curv')) + def test_FaceIP_3D_isotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Curv')) + def test_FaceIP_2D_anisotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4],2, False, 'Curv')) + def test_FaceIP_3D_anisotropic_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Curv')) + def test_FaceIP_2D_tensor_Curv(self): + self.assertTrue(self.doTestFace([10, 4],3, False, 'Curv')) + def test_FaceIP_3D_tensor_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Curv')) - def test_FaceIP_2D_float_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM')) - def test_FaceIP_3D_float_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM')) - def test_FaceIP_2D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM')) - def test_FaceIP_3D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM')) - def test_FaceIP_2D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM')) - def test_FaceIP_3D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM')) + def test_FaceIP_2D_float_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4],0, True, 'Curv')) + def test_FaceIP_3D_float_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Curv')) + def test_FaceIP_2D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4],1, True, 'Curv')) + def test_FaceIP_3D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Curv')) + def test_FaceIP_2D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4],2, True, 'Curv')) + def test_FaceIP_3D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Curv')) - def test_EdgeIP_2D_float_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM')) - def test_EdgeIP_3D_float_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM')) - def test_EdgeIP_2D_isotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM')) - def test_EdgeIP_3D_isotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM')) - def test_EdgeIP_2D_anisotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM')) - def test_EdgeIP_3D_anisotropic_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM')) - def test_EdgeIP_2D_tensor_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM')) - def test_EdgeIP_3D_tensor_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM')) + def test_EdgeIP_2D_float_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],0, False, 'Curv')) + def test_EdgeIP_3D_float_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Curv')) + def test_EdgeIP_2D_isotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],1, False, 'Curv')) + def test_EdgeIP_3D_isotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Curv')) + def test_EdgeIP_2D_anisotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],2, False, 'Curv')) + def test_EdgeIP_3D_anisotropic_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Curv')) + def test_EdgeIP_2D_tensor_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],3, False, 'Curv')) + def test_EdgeIP_3D_tensor_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Curv')) - def test_EdgeIP_2D_float_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM')) - def test_EdgeIP_3D_float_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM')) - def test_EdgeIP_2D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM')) - def test_EdgeIP_3D_isotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM')) - def test_EdgeIP_2D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM')) - def test_EdgeIP_3D_anisotropic_fast_LRM(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM')) + def test_EdgeIP_2D_float_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],0, True, 'Curv')) + def test_EdgeIP_3D_float_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Curv')) + def test_EdgeIP_2D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],1, True, 'Curv')) + def test_EdgeIP_3D_isotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Curv')) + def test_EdgeIP_2D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4],2, True, 'Curv')) + def test_EdgeIP_3D_anisotropic_fast_Curv(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Curv')) diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index 9a465fcd..416689bc 100644 --- a/SimPEG/Tests/test_operators.py +++ b/SimPEG/Tests/test_operators.py @@ -4,7 +4,7 @@ from TestUtils import OrderTest import matplotlib.pyplot as plt #TODO: 'randomTensorMesh' -MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] +MESHTYPES = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -38,7 +38,7 @@ class TestCurl(OrderTest): curlE_ana = self.M.projectFaceVector(Fc) curlE = self.M.edgeCurl.dot(E) - if self._meshType == 'rotateLRM': + if self._meshType == 'rotateCurv': # Really it is the integration we should be caring about: # So, let us look at the l2 norm. err = np.linalg.norm(self.M.area*(curlE - curlE_ana), 2) @@ -229,7 +229,7 @@ class TestFaceDiv3D(OrderTest): divF = self.M.faceDiv.dot(F) divF_ana = call3(sol, self.M.gridCC) - if self._meshType == 'rotateLRM': + if self._meshType == 'rotateCurv': # Really it is the integration we should be caring about: # So, let us look at the l2 norm. err = np.linalg.norm(self.M.vol*(divF-divF_ana), 2) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 6b04e042..6637a138 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,7 +1,7 @@ from matutils import * from codeutils import * from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile -from lrmutils import volTetra, faceInfo, indexCube +from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate from CounterUtils import * diff --git a/SimPEG/Utils/lrmutils.py b/SimPEG/Utils/curvutils.py similarity index 100% rename from SimPEG/Utils/lrmutils.py rename to SimPEG/Utils/curvutils.py diff --git a/docs/SimPEGFrameworkRevised.png b/docs/SimPEGFrameworkRevised.png new file mode 100644 index 00000000..18ce7a30 Binary files /dev/null and b/docs/SimPEGFrameworkRevised.png differ diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index c3fbf4d6..bc9426c6 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -75,7 +75,7 @@ We multiply by square-root of volume on each side of the tensor conductivity to .. math:: \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\ - \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM} + \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv} Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above. We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like: @@ -114,7 +114,7 @@ Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of .. math:: - \mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)} + \mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{Curv only}} \mathbf{Q}_{(i)} .. note:: diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index c698cf67..7bf398b1 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -29,7 +29,7 @@ the implementations. tM = Mesh.TensorMesh(sz) qM = Mesh.TreeMesh(sz) qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0) - rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) fig, axes = plt.subplots(1,3,figsize=(14,5)) opts = {} @@ -38,7 +38,7 @@ the implementations. qM.plotGrid(ax=axes[1], **opts) axes[1].set_title('TreeMesh') rM.plotGrid(ax=axes[2], **opts) - axes[2].set_title('LogicallyRectMesh') + axes[2].set_title('CurvilinearMesh') plt.show() diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst index 13ceffba..9fb52e53 100644 --- a/docs/api_MeshCode.rst +++ b/docs/api_MeshCode.rst @@ -21,7 +21,7 @@ Tree Mesh Logically Rectangular Mesh ========================== -.. automodule:: SimPEG.Mesh.LogicallyRectMesh +.. automodule:: SimPEG.Mesh.Curvilinear :show-inheritance: :members: :undoc-members: diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 7903cdc4..042aef67 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -23,10 +23,10 @@ Solver Utilities :members: :undoc-members: -LRM Utilities +Curv Utilities ============= -.. automodule:: SimPEG.Utils.lrmutils +.. automodule:: SimPEG.Utils.curvutils :members: :undoc-members: