diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index cecf6d7b..9af638b8 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -399,6 +399,94 @@ class BaseRectangularMesh(BaseMesh): """ return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None]) + ################################## + # Redo the numbering so they are dependent of the vector numbers + ################################## + + @property + def nC(self): + """ + Total number of cells + + :rtype: int + :return: nC + """ + return self.vnC.prod() + + @property + def nN(self): + """ + Total number of nodes + + :rtype: int + :return: nN + """ + return self.vnN.prod() + + @property + def nEx(self): + """ + Number of x-edges + + :rtype: int + :return: nEx + """ + return self.vnEx.prod() + + @property + def nEy(self): + """ + Number of y-edges + + :rtype: int + :return: nEy + """ + if self.dim < 2: return + return self.vnEy.prod() + + @property + def nEz(self): + """ + Number of z-edges + + :rtype: int + :return: nEz + """ + if self.dim < 3: return + return self.vnEz.prod() + + @property + def nFx(self): + """ + Number of x-faces + + :rtype: int + :return: nFx + """ + return self.vnFx.prod() + + @property + def nFy(self): + """ + Number of y-faces + + :rtype: int + :return: nFy + """ + if self.dim < 2: return + return self.vnFy.prod() + + @property + def nFz(self): + """ + Number of z-faces + + :rtype: int + :return: nFz + """ + if self.dim < 3: return + return self.vnFz.prod() + def r(self, x, xType='CC', outType='CC', format='V'): """ Mesh.r is a quick reshape command that will do the best it can at giving you what you want. diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 226b8720..0502ed31 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -26,8 +26,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNx """ - if self.nCy == 1: - return self.nCx + if self.nCy == 1: return self.nCx return self.nCx + 1 @property @@ -38,30 +37,9 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNy """ - if self.nCy == 1: - return self.nCy - 1 + if self.nCy == 1: return 0 return self.nCy - @property - def nN(self): - """ - Total number of nodes - - :rtype: int - :return: nN - """ - return (np.r_[self.nNx, self.nNy, self.nNz]).prod() - - @property - def nFx(self): - """ - Number of x-faces - - :rtype: int - :return: nFx - """ - return self.nC - @property def vnFx(self): """ @@ -73,44 +51,15 @@ class CylMesh(BaseTensorMesh): return self.vnC @property - def nFy(self): + def vnEy(self): """ - Number of y-faces + Number of y-edges in each direction - :rtype: int - :return: nFy + :rtype: numpy.array (dim, ) + :return: vnEy or None if dim < 2 """ - return (self.vnC + np.r_[0,-1,0][:self.dim]).prod() - - @property - def nEx(self): - """ - Number of x-edges - - :rtype: int - :return: nEx - """ - return (self._n + np.r_[0,-1,1]).prod() - - @property - def nEy(self): - """ - Number of y-edges - - :rtype: int - :return: nEy - """ - return (self._n + np.r_[0,0,1]).prod() - - @property - def nEz(self): - """ - Number of z-edges - - :rtype: int - :return: nEz - """ - return (self._n + np.r_[0,-1,0]).prod() + nNx = self.nNx if self.nCy == 1 else self.nNx - 1 + return np.r_[nNx, self.nCy, self.nNz] @property def vectorCCx(self): @@ -232,8 +181,8 @@ class CylMesh(BaseTensorMesh): @property def nodalGrad(self): """Construct gradient operator (nodes to edges).""" - if self.nCy == 1: - raise Exception('Nodal grad does not make sense for cylindrically symmetric mesh.') + # Nodal grad does not make sense for cylindrically symmetric mesh. + if self.nCy == 1: return None raise NotImplementedError('nodalGrad not yet implemented') @property diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 70c34c88..18db1686 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -11,21 +11,24 @@ class TestCyl2DMesh(unittest.TestCase): hz = np.r_[2,1] self.mesh = Mesh.CylMesh([hx, 1,hz]) - def test_cylMesh_numbers(self): + def test_dim(self): self.assertTrue(self.mesh.dim == 3) + def test_nC(self): self.assertTrue(self.mesh.nC == 6) self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 1) self.assertTrue(self.mesh.nCz == 2) self.assertTrue(np.all(self.mesh.vnC == [3, 1, 2])) + def test_nN(self): self.assertTrue(self.mesh.nN == 0) self.assertTrue(self.mesh.nNx == 3) self.assertTrue(self.mesh.nNy == 0) self.assertTrue(self.mesh.nNz == 3) self.assertTrue(np.all(self.mesh.vnN == [3, 0, 3])) + def test_nF(self): self.assertTrue(self.mesh.nFx == 6) self.assertTrue(np.all(self.mesh.vnFx == [3, 1, 2])) self.assertTrue(self.mesh.nFy == 0) @@ -35,6 +38,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(self.mesh.nF == 15) self.assertTrue(np.all(self.mesh.vnF == [6, 0, 9])) + def test_nE(self): self.assertTrue(self.mesh.nEx == 0) self.assertTrue(np.all(self.mesh.vnEx == [3, 0, 3])) self.assertTrue(self.mesh.nEy == 9) @@ -126,6 +130,9 @@ class TestCyl2DMesh(unittest.TestCase): G = np.c_[x,y,z] self.assertTrue(np.linalg.norm((G-self.mesh.gridEy).ravel()) == 0) + def test_lightOperators(self): + self.assertTrue(self.mesh.nodalGrad is None) + MESHTYPES = ['uniformCylMesh'] @@ -178,18 +185,42 @@ class TestCyl3DMesh(unittest.TestCase): hz = np.r_[2,1] self.mesh = Mesh.CylMesh([hx, hy,hz]) - def test_cylMesh_numbers(self): + def test_dim(self): + self.assertTrue(self.mesh.dim == 3) + + def test_nC(self): self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 2) self.assertTrue(self.mesh.nCz == 2) self.assertTrue(np.all(self.mesh.vnC == [3, 2, 2])) + def test_nN(self): self.assertTrue(self.mesh.nN == 24) self.assertTrue(self.mesh.nNx == 4) self.assertTrue(self.mesh.nNy == 2) self.assertTrue(self.mesh.nNz == 3) self.assertTrue(np.all(self.mesh.vnN == [4, 2, 3])) + def test_nF(self): + self.assertTrue(self.mesh.nFx == 12) + self.assertTrue(np.all(self.mesh.vnFx == [3, 2, 2])) + self.assertTrue(self.mesh.nFy == 12) + self.assertTrue(np.all(self.mesh.vnFy == [3, 2, 2])) + self.assertTrue(self.mesh.nFz == 18) + self.assertTrue(np.all(self.mesh.vnFz == [3, 2, 3])) + self.assertTrue(self.mesh.nF == 42) + self.assertTrue(np.all(self.mesh.vnF == [12, 12, 18])) + + def test_nE(self): + self.assertTrue(self.mesh.nEx == 18) + self.assertTrue(np.all(self.mesh.vnEx == [3, 2, 3])) + self.assertTrue(self.mesh.nEy == 18) + self.assertTrue(np.all(self.mesh.vnEy == [3, 2, 3])) + # self.assertTrue(self.mesh.nEz == 0) + # self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) + # self.assertTrue(self.mesh.nE == 9) + # self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0])) + def test_vectorsCC(self): v = np.r_[0.5, 1.5, 2.25] self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0)