Split out BaseRectangularMesh (Issue #48)

This commit is contained in:
rowanc1
2014-02-14 11:07:25 -08:00
parent ea92cf3fc0
commit 39b250c9ac
5 changed files with 375 additions and 368 deletions
+365 -358
View File
@@ -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)
+3 -3
View File
@@ -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))
+3 -3
View File
@@ -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
+1 -1
View File
@@ -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
+3 -3
View File
@@ -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)