From 39b250c9ac8b42b94cbec8f3a09524feb2992ab0 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 14 Feb 2014 11:07:25 -0800 Subject: [PATCH] Split out BaseRectangularMesh (Issue #48) --- SimPEG/Mesh/BaseMesh.py | 723 +++++++++++++------------ SimPEG/Mesh/LogicallyOrthogonalMesh.py | 6 +- SimPEG/Mesh/TensorMesh.py | 6 +- SimPEG/Mesh/__init__.py | 2 +- SimPEG/Tests/test_basemesh.py | 6 +- 5 files changed, 375 insertions(+), 368 deletions(-) diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index bf6655c9..aa487415 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -38,6 +38,371 @@ class BaseMesh(object): """ return self._x0 + @property + def dim(self): + """ + The dimension of the mesh (1, 2, or 3). + + :rtype: int + :return: dim + """ + return len(self._n) + + @property + def nC(self): + """ + Total number of cells in the model. + + :rtype: int + :return: nC + + .. plot:: + :include-source: + + from SimPEG import Mesh, np + Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True) + """ + return self._n.prod() + + @property + def nN(self): + """ + Total number of nodes + + :rtype: int + :return: nN + + .. plot:: + :include-source: + + from SimPEG import Mesh, np + Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True) + """ + return (self._n+1).prod() + + @property + def nEx(self): + """ + Number of x-edges + + :rtype: int + :return: nEx + """ + return (self._n + np.r_[0,1,1][:self.dim]).prod() + + @property + def nEy(self): + """ + Number of y-edges + + :rtype: int + :return: nEy + """ + if self.dim < 2: return None + return (self._n + np.r_[1,0,1][:self.dim]).prod() + + @property + def nEz(self): + """ + Number of z-edges + + :rtype: int + :return: nEz + """ + if self.dim < 3: return None + return (self._n + np.r_[1,1,0][:self.dim]).prod() + + @property + def vnE(self): + """ + Total number of edges in each direction + + :rtype: numpy.array (dim, ) + :return: [nEx, nEy, nEz] + + .. plot:: + :include-source: + + from SimPEG import Mesh, np + Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True) + """ + return np.array([x for x in [self.nEx, self.nEy, self.nEz] if not x is None]) + + @property + def nE(self): + """ + Total number of edges. + + :rtype: int + :return: sum([nEx, nEy, nEz]) + + """ + return self.vnE.sum() + + @property + def nFx(self): + """ + Number of x-faces + + :rtype: int + :return: nFx + """ + return (self._n + np.r_[1,0,0][:self.dim]).prod() + + @property + def nFy(self): + """ + Number of y-faces + + :rtype: int + :return: nFy + """ + if self.dim < 2: return None + return (self._n + np.r_[0,1,0][:self.dim]).prod() + + @property + def nFz(self): + """ + Number of z-faces + + :rtype: int + :return: nFz + """ + if self.dim < 3: return None + return (self._n + np.r_[0,0,1][:self.dim]).prod() + + @property + def vnF(self): + """ + Total number of faces in each direction + + :rtype: numpy.array (dim, ) + :return: [nFx, nFy, nFz] + + .. plot:: + :include-source: + + from SimPEG import Mesh, np + Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True) + """ + return np.array([x for x in [self.nFx, self.nFy, self.nFz] if not x is None]) + + @property + def nF(self): + """ + Total number of faces. + + :rtype: int + :return: sum([nFx, nFy, nFz]) + + """ + return self.vnF.sum() + + @property + def normals(self): + """ + Face Normals + + :rtype: numpy.array (sum(nF), dim) + :return: normals + """ + if self.dim == 2: + nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)] + nY = np.c_[np.zeros(self.nFy), np.ones(self.nFy)] + return np.r_[nX, nY] + elif self.dim == 3: + nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx), np.zeros(self.nFx)] + nY = np.c_[np.zeros(self.nFy), np.ones(self.nFy), np.zeros(self.nFy)] + nZ = np.c_[np.zeros(self.nFz), np.zeros(self.nFz), np.ones(self.nFz)] + return np.r_[nX, nY, nZ] + + @property + def tangents(self): + """ + Edge Tangents + + :rtype: numpy.array (sum(nE), dim) + :return: normals + """ + if self.dim == 2: + tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)] + tY = np.c_[np.zeros(self.nEy), np.ones(self.nEy)] + return np.r_[tX, tY] + elif self.dim == 3: + tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx), np.zeros(self.nEx)] + tY = np.c_[np.zeros(self.nEy), np.ones(self.nEy), np.zeros(self.nEy)] + tZ = np.c_[np.zeros(self.nEz), np.zeros(self.nEz), np.ones(self.nEz)] + return np.r_[tX, tY, tZ] + + def projectFaceVector(self, fV): + """ + Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals + + :param numpy.array fV: face vector with shape (nF, dim) + :rtype: numpy.array with shape (nF, ) + :return: projected face vector + """ + 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 np.sum(fV*self.normals, 1) + + def projectEdgeVector(self, eV): + """ + Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents + + :param numpy.array eV: edge vector with shape (nE, dim) + :rtype: numpy.array with shape (nE, ) + :return: projected edge vector + """ + 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 np.sum(eV*self.tangents, 1) + + +class BaseRectangularMesh(BaseMesh): + """BaseRectangularMesh""" + def __init__(self, n, x0=None): + BaseMesh.__init__(self, n, x0) + + @property + def nCx(self): + """ + Number of cells in the x direction + + :rtype: int + :return: nCx + """ + return self._n[0] + + @property + def nCy(self): + """ + Number of cells in the y direction + + :rtype: int + :return: nCy or None if dim < 2 + """ + return None if self.dim < 2 else self._n[1] + + @property + def nCz(self): + """Number of cells in the z direction + + :rtype: int + :return: nCz or None if dim < 3 + """ + return None if self.dim < 3 else self._n[2] + + @property + def vnC(self): + """ + Total number of cells in each direction + + :rtype: numpy.array (dim, ) + :return: [nCx, nCy, nCz] + """ + return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None]) + + @property + def nNx(self): + """ + Number of nodes in the x-direction + + :rtype: int + :return: nNx + """ + return self.nCx + 1 + + @property + def nNy(self): + """ + Number of noes in the y-direction + + :rtype: int + :return: nNy or None if dim < 2 + """ + return None if self.dim < 2 else self.nCy + 1 + + @property + def nNz(self): + """ + Number of nodes in the z-direction + + :rtype: int + :return: nNz or None if dim < 3 + """ + return None if self.dim < 3 else self.nCz + 1 + + @property + def vnN(self): + """ + Total number of nodes in each direction + + :rtype: numpy.array (dim, ) + :return: [nNx, nNy, nNz] + """ + return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None]) + + @property + def vnEx(self): + """ + Number of x-edges in each direction + + :rtype: numpy.array (dim, ) + :return: vnEx + """ + return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None]) + + @property + def vnEy(self): + """ + Number of y-edges in each direction + + :rtype: numpy.array (dim, ) + :return: vnEy or None if dim < 2 + """ + return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None]) + + @property + def vnEz(self): + """ + Number of z-edges in each direction + + :rtype: numpy.array (dim, ) + :return: vnEz or None if dim < 3 + """ + return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None]) + + @property + def vnFx(self): + """ + Number of x-faces in each direction + + :rtype: numpy.array (dim, ) + :return: vnFx + """ + return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None]) + + @property + def vnFy(self): + """ + Number of y-faces in each direction + + :rtype: numpy.array (dim, ) + :return: vnFy or None if dim < 2 + """ + return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None]) + + @property + def vnFz(self): + """ + Number of z-faces in each direction + + :rtype: numpy.array (dim, ) + :return: vnFz or None if dim < 3 + """ + return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None]) + 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. @@ -144,361 +509,3 @@ class BaseMesh(object): return out else: return switchKernal(x) - - @property - def dim(self): - """ - The dimension of the mesh (1, 2, or 3). - - :rtype: int - :return: dim - """ - return len(self._n) - - @property - def nCx(self): - """ - Number of cells in the x direction - - :rtype: int - :return: nCx - """ - return self._n[0] - - @property - def nCy(self): - """ - Number of cells in the y direction - - :rtype: int - :return: nCy or None if dim < 2 - """ - return None if self.dim < 2 else self._n[1] - - @property - def nCz(self): - """Number of cells in the z direction - - :rtype: int - :return: nCz or None if dim < 3 - """ - return None if self.dim < 3 else self._n[2] - - @property - def vnC(self): - """ - Total number of cells in each direction - - :rtype: numpy.array (dim, ) - :return: [nCx, nCy, nCz] - """ - return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None]) - - @property - def nC(self): - doc = """ - Total number of cells in the model. - - :rtype: int - :return: nC - - .. plot:: - :include-source: - - from SimPEG import Mesh, np - Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True) - """ - return self.vnC.prod() - - @property - def nNx(self): - """ - Number of nodes in the x-direction - - :rtype: int - :return: nNx - """ - return self.nCx + 1 - - @property - def nNy(self): - """ - Number of noes in the y-direction - - :rtype: int - :return: nNy or None if dim < 2 - """ - return None if self.dim < 2 else self.nCy + 1 - - @property - def nNz(self): - """ - Number of nodes in the z-direction - - :rtype: int - :return: nNz or None if dim < 3 - """ - return None if self.dim < 3 else self.nCz + 1 - - @property - def vnN(self): - """ - Total number of nodes in each direction - - :rtype: numpy.array (dim, ) - :return: [nNx, nNy, nNz] - """ - return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None]) - - @property - def nN(self): - doc = """ - Total number of nodes - - :rtype: int - :return: nN - - .. plot:: - :include-source: - - from SimPEG import Mesh, np - Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True) - """ - 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 - """ - return self.vnEy.prod() - - @property - def nEz(self): - """ - Number of z-edges - - :rtype: int - :return: nEz - """ - return self.vnEz.prod() - - @property - def vnEx(self): - """ - Number of x-edges in each direction - - :rtype: numpy.array (dim, ) - :return: vnEx - """ - return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None]) - - @property - def vnEy(self): - """ - Number of y-edges in each direction - - :rtype: numpy.array (dim, ) - :return: vnEy or None if dim < 2 - """ - return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None]) - - @property - def vnEz(self): - """ - Number of z-edges in each direction - - :rtype: numpy.array (dim, ) - :return: vnEz or None if dim < 3 - """ - return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None]) - - @property - def vnE(self): - """ - Total number of edges in each direction - - :rtype: numpy.array (dim, ) - :return: [prod(vnEx), prod(vnEy), prod(vnEz)] - - .. plot:: - :include-source: - - from SimPEG import Mesh, np - Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True) - """ - return np.array([np.prod(x) for x in [self.vnEx, self.vnEy, self.vnEz] if not x is None]) - - - @property - def nE(self): - """ - Total number of edges. - - :rtype: int - :return: sum([prod(vnEx), prod(vnEy), prod(vnEz)]) - - """ - return self.vnE.sum() - - @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 - """ - return self.vnFy.prod() - - @property - def nFz(self): - """ - Number of z-faces - - :rtype: int - :return: nFz - """ - return self.vnFz.prod() - - @property - def vnFx(self): - """ - Number of x-faces in each direction - - :rtype: numpy.array (dim, ) - :return: vnFx - """ - return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None]) - - @property - def vnFy(self): - """ - Number of y-faces in each direction - - :rtype: numpy.array (dim, ) - :return: vnFy or None if dim < 2 - """ - return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None]) - - @property - def vnFz(self): - """ - Number of z-faces in each direction - - :rtype: numpy.array (dim, ) - :return: vnFz or None if dim < 3 - """ - return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None]) - - @property - def vnF(self): - """ - Total number of faces in each direction - - :rtype: numpy.array (dim, ) - :return: [prod(vnFx), prod(vnFy), prod(vnFz)] - - .. plot:: - :include-source: - - from SimPEG import Mesh, np - Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True) - """ - return np.array([np.prod(x) for x in [self.vnFx, self.vnFy, self.vnFz] if not x is None]) - - - @property - def nF(self): - """ - Total number of faces. - - :rtype: int - :return: sum([vnFx, vnFy, vnFz]) - - """ - return self.vnF.sum() - - @property - def normals(self): - """ - Face Normals - - :rtype: numpy.array (sum(nF), dim) - :return: normals - """ - if self.dim == 2: - nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)] - nY = np.c_[np.zeros(self.nFy), np.ones(self.nFy)] - return np.r_[nX, nY] - elif self.dim == 3: - nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx), np.zeros(self.nFx)] - nY = np.c_[np.zeros(self.nFy), np.ones(self.nFy), np.zeros(self.nFy)] - nZ = np.c_[np.zeros(self.nFz), np.zeros(self.nFz), np.ones(self.nFz)] - return np.r_[nX, nY, nZ] - - @property - def tangents(self): - """ - Edge Tangents - - :rtype: numpy.array (sum(nE), dim) - :return: normals - """ - if self.dim == 2: - tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)] - tY = np.c_[np.zeros(self.nEy), np.ones(self.nEy)] - return np.r_[tX, tY] - elif self.dim == 3: - tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx), np.zeros(self.nEx)] - tY = np.c_[np.zeros(self.nEy), np.ones(self.nEy), np.zeros(self.nEy)] - tZ = np.c_[np.zeros(self.nEz), np.zeros(self.nEz), np.ones(self.nEz)] - return np.r_[tX, tY, tZ] - - - def projectFaceVector(self, fV): - """ - Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals - - :param numpy.array fV: face vector with shape (nF, dim) - :rtype: numpy.array with shape (nF, ) - :return: projected face vector - """ - 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 np.sum(fV*self.normals, 1) - - def projectEdgeVector(self, eV): - """ - Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents - - :param numpy.array eV: edge vector with shape (nE, dim) - :rtype: numpy.array with shape (nE, ) - :return: projected edge vector - """ - 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 np.sum(eV*self.tangents, 1) diff --git a/SimPEG/Mesh/LogicallyOrthogonalMesh.py b/SimPEG/Mesh/LogicallyOrthogonalMesh.py index 94f34a71..480c4f06 100644 --- a/SimPEG/Mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/Mesh/LogicallyOrthogonalMesh.py @@ -1,5 +1,5 @@ from SimPEG import Utils, np -from BaseMesh import BaseMesh +from BaseMesh import BaseRectangularMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts from LomView import LomView @@ -11,7 +11,7 @@ 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 LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): +class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts, LomView): """ LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes. @@ -40,7 +40,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): assert len(nodes[0].shape) == len(nodes), "Dimension mismatch" assert len(nodes[0].shape) > 1, "Not worth using LOM for a 1D mesh." - super(LogicallyOrthogonalMesh, self).__init__(np.array(nodes[0].shape)-1, None) + BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None) # Save nodes to private variable _gridN as vectors self._gridN = np.ones((nodes[0].size, self.dim)) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 0c50ba08..12253bb5 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -1,10 +1,10 @@ from SimPEG import Utils, np, sp -from BaseMesh import BaseMesh +from BaseMesh import BaseRectangularMesh from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): +class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -48,7 +48,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) h[i] = h_i[:] # make a copy. - BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" # Ensure h contains 1D vectors diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 3da22a01..47b30243 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -2,7 +2,7 @@ from Cyl1DMesh import Cyl1DMesh from TensorMesh import TensorMesh from TreeMesh import TreeMesh from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh -from BaseMesh import BaseMesh +from BaseMesh import BaseMesh, BaseRectangularMesh from TensorView import TensorView from LomView import LomView from InnerProducts import InnerProducts diff --git a/SimPEG/Tests/test_basemesh.py b/SimPEG/Tests/test_basemesh.py index 6f4c3b5d..eb3acc21 100644 --- a/SimPEG/Tests/test_basemesh.py +++ b/SimPEG/Tests/test_basemesh.py @@ -1,13 +1,13 @@ import unittest import sys -from SimPEG.Mesh import BaseMesh +from SimPEG.Mesh import BaseRectangularMesh import numpy as np class TestBaseMesh(unittest.TestCase): def setUp(self): - self.mesh = BaseMesh([6, 2, 3]) + self.mesh = BaseRectangularMesh([6, 2, 3]) def test_meshDimensions(self): self.assertTrue(self.mesh.dim, 3) @@ -109,7 +109,7 @@ class TestBaseMesh(unittest.TestCase): class TestMeshNumbers2D(unittest.TestCase): def setUp(self): - self.mesh = BaseMesh([6, 2]) + self.mesh = BaseRectangularMesh([6, 2]) def test_meshDimensions(self): self.assertTrue(self.mesh.dim, 2)