mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 09:24:25 +08:00
Merge pull request #55 from simpeg/nCv2vnC
Refactoring of BaseMesh Change nCv etc. to vnC
This commit is contained in:
Rowan Cockett
+372
-305
@@ -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.
|
||||
@@ -100,13 +465,13 @@ class BaseMesh(object):
|
||||
elif xType in ['F', 'E']:
|
||||
# This will only deal with components of fields, not full 'F' or 'E'
|
||||
xx = Utils.mkvc(xx) # unwrap it in case it is a matrix
|
||||
nn = self.nFv if xType == 'F' else self.nEv
|
||||
nn = self.vnF if xType == 'F' else self.vnE
|
||||
nn = np.r_[0, nn]
|
||||
|
||||
nx = [0, 0, 0]
|
||||
nx[0] = self.nFx if xType == 'F' else self.nEx
|
||||
nx[1] = self.nFy if xType == 'F' else self.nEy
|
||||
nx[2] = self.nFz if xType == 'F' else self.nEz
|
||||
nx[0] = self.vnFx if xType == 'F' else self.vnEx
|
||||
nx[1] = self.vnFy if xType == 'F' else self.vnEy
|
||||
nx[2] = self.vnFz if xType == 'F' else self.vnEz
|
||||
|
||||
for dim, dimName in enumerate(['x', 'y', 'z']):
|
||||
if dimName in outType:
|
||||
@@ -118,11 +483,11 @@ class BaseMesh(object):
|
||||
elif xTypeIsFExyz:
|
||||
# This will deal with partial components (x, y or z) lying on edges or faces
|
||||
if 'x' in xType:
|
||||
nn = self.nFx if 'F' in xType else self.nEx
|
||||
nn = self.vnFx if 'F' in xType else self.vnEx
|
||||
elif 'y' in xType:
|
||||
nn = self.nFy if 'F' in xType else self.nEy
|
||||
nn = self.vnFy if 'F' in xType else self.vnEy
|
||||
elif 'z' in xType:
|
||||
nn = self.nFz if 'F' in xType else self.nEz
|
||||
nn = self.vnFz if 'F' in xType else self.vnEz
|
||||
assert xx.size == np.prod(nn), 'Vector is not the right size.'
|
||||
return outKernal(xx, nn)
|
||||
|
||||
@@ -144,301 +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 nCv(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.nCv.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 nNv(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.nNv.prod()
|
||||
|
||||
@property
|
||||
def nEx(self):
|
||||
"""
|
||||
Number of x-edges in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: nEx
|
||||
"""
|
||||
return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
|
||||
|
||||
@property
|
||||
def nEy(self):
|
||||
"""
|
||||
Number of y-edges in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: nEy 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 nEz(self):
|
||||
"""
|
||||
Number of z-edges in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: nEz 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 nEv(self):
|
||||
"""
|
||||
Total number of edges in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: [prod(nEx), prod(nEy), prod(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([np.prod(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([prod(nEx), prod(nEy), prod(nEz)])
|
||||
|
||||
"""
|
||||
return self.nEv.sum()
|
||||
|
||||
@property
|
||||
def nFx(self):
|
||||
"""
|
||||
Number of x-faces in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: nFx
|
||||
"""
|
||||
return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
|
||||
|
||||
@property
|
||||
def nFy(self):
|
||||
"""
|
||||
Number of y-faces in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: nFy 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 nFz(self):
|
||||
"""
|
||||
Number of z-faces in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: nFz 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 nFv(self):
|
||||
"""
|
||||
Total number of faces in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: [prod(nFx), prod(nFy), prod(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([np.prod(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.nFv.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.nFv[0]), np.zeros(self.nFv[0])]
|
||||
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])]
|
||||
return np.r_[nX, nY]
|
||||
elif self.dim == 3:
|
||||
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])]
|
||||
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])]
|
||||
nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])]
|
||||
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.nEv[0]), np.zeros(self.nEv[0])]
|
||||
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])]
|
||||
return np.r_[tX, tY]
|
||||
elif self.dim == 3:
|
||||
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])]
|
||||
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])]
|
||||
tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])]
|
||||
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)
|
||||
|
||||
@@ -83,11 +83,11 @@ class Cyl1DMesh(object):
|
||||
return locals()
|
||||
nC = property(**nC())
|
||||
|
||||
def nCv():
|
||||
def vnC():
|
||||
doc = "Total number of cells in each direction"
|
||||
fget = lambda self: np.array([self.nCx, self.nCz])
|
||||
return locals()
|
||||
nCv = property(**nCv())
|
||||
vnC = property(**vnC())
|
||||
|
||||
def nNr():
|
||||
doc = "Number of nodes in the radial direction"
|
||||
@@ -113,21 +113,21 @@ class Cyl1DMesh(object):
|
||||
return locals()
|
||||
nFr = property(**nFr())
|
||||
|
||||
def nFz():
|
||||
def vnFz():
|
||||
doc = "Number of z faces"
|
||||
fget = lambda self: self.nNz * self.nCx
|
||||
return locals()
|
||||
nFz = property(**nFz())
|
||||
vnFz = property(**vnFz())
|
||||
|
||||
def nFv():
|
||||
def vnF():
|
||||
doc = "Total number of faces in each direction"
|
||||
fget = lambda self: np.array([self.nFr, self.nFz])
|
||||
fget = lambda self: np.array([self.nFr, self.vnFz])
|
||||
return locals()
|
||||
nFv = property(**nFv())
|
||||
vnF = property(**vnF())
|
||||
|
||||
def nF():
|
||||
doc = "Total number of faces"
|
||||
fget = lambda self: self.nFr + self.nFz
|
||||
fget = lambda self: self.nFr + self.vnFz
|
||||
return locals()
|
||||
nF = property(**nF())
|
||||
|
||||
|
||||
@@ -129,7 +129,7 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if(self._faceDiv is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
if(self.dim == 1):
|
||||
D = ddx(n[0])
|
||||
@@ -158,7 +158,7 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if(self._faceDivx is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
if(self.dim == 1):
|
||||
D1 = ddx(n[0])
|
||||
@@ -183,7 +183,7 @@ class DiffOperators(object):
|
||||
if(self.dim < 2): return None
|
||||
if(self._faceDivy is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
if(self.dim == 2):
|
||||
D2 = sp.kron(ddx(n[1]), speye(n[0]))
|
||||
@@ -225,7 +225,7 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if(self._nodalGrad is None):
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute divergence operator on faces
|
||||
if(self.dim == 1):
|
||||
G = ddx(n[0])
|
||||
@@ -253,7 +253,7 @@ class DiffOperators(object):
|
||||
if(self._nodalLaplacian is None):
|
||||
print 'Warning: Laplacian has not been tested rigorously.'
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute divergence operator on faces
|
||||
if(self.dim == 1):
|
||||
D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
|
||||
@@ -291,7 +291,7 @@ class DiffOperators(object):
|
||||
|
||||
"""
|
||||
if(type(BC) is str):
|
||||
BC = [BC for _ in self.nCv] # Repeat the str self.dim times
|
||||
BC = [BC for _ in self.vnC] # Repeat the str self.dim times
|
||||
elif(type(BC) is list):
|
||||
assert len(BC) == self.dim, 'BC list must be the size of your mesh'
|
||||
else:
|
||||
@@ -313,7 +313,7 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if(self._cellGrad is None):
|
||||
BC = self.setCellGradBC(self._cellGradBC_list)
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
G = ddxCellGrad(n[0], BC[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -340,7 +340,7 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if(self._cellGradBC is None):
|
||||
BC = self.setCellGradBC(self._cellGradBC_list)
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
G = ddxCellGradBC(n[0], BC[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -367,7 +367,7 @@ class DiffOperators(object):
|
||||
def fget(self):
|
||||
if getattr(self, '_cellGradx', None) is None:
|
||||
BC = ['neumann', 'neumann']
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
G1 = ddxCellGrad(n[0], BC)
|
||||
elif(self.dim == 2):
|
||||
@@ -388,7 +388,7 @@ class DiffOperators(object):
|
||||
if self.dim < 2: return None
|
||||
if getattr(self, '_cellGrady', None) is None:
|
||||
BC = ['neumann', 'neumann']
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 2):
|
||||
G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
|
||||
elif(self.dim == 3):
|
||||
@@ -466,7 +466,7 @@ class DiffOperators(object):
|
||||
def aveF2CC(self):
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
if getattr(self, '_aveF2CC', None) is None:
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveF2CC = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -483,7 +483,7 @@ class DiffOperators(object):
|
||||
def aveF2CCV(self):
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
if getattr(self, '_aveF2CCV', None) is None:
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveF2CCV = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -499,7 +499,7 @@ class DiffOperators(object):
|
||||
def aveCC2F(self):
|
||||
"Construct the averaging operator on cell cell centers to faces."
|
||||
if getattr(self, '_aveCC2F', None) is None:
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveCC2F = avExtrap(n[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -516,7 +516,7 @@ class DiffOperators(object):
|
||||
"Construct the averaging operator on cell edges to cell centers."
|
||||
if getattr(self, '_aveE2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
elif(self.dim == 2):
|
||||
@@ -533,7 +533,7 @@ class DiffOperators(object):
|
||||
"Construct the averaging operator on cell edges to cell centers."
|
||||
if getattr(self, '_aveE2CCV', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
elif(self.dim == 2):
|
||||
@@ -550,7 +550,7 @@ class DiffOperators(object):
|
||||
"Construct the averaging operator on cell nodes to cell centers."
|
||||
if getattr(self, '_aveN2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveN2CC = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -565,7 +565,7 @@ class DiffOperators(object):
|
||||
|
||||
if getattr(self, '_aveN2E', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveN2E = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
@@ -582,7 +582,7 @@ class DiffOperators(object):
|
||||
"Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
|
||||
if getattr(self, '_aveN2F', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveN2F = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
|
||||
@@ -403,8 +403,8 @@ def _getFacePxx_Rectangular(M):
|
||||
posFx = 0 if xFace == 'fXm' else 1
|
||||
posFy = 0 if yFace == 'fYm' else 1
|
||||
|
||||
ind1 = sub2ind(M.nFx, np.c_[ii + posFx, jj])
|
||||
ind2 = sub2ind(M.nFy, np.c_[ii, jj + posFy]) + M.nFv[0]
|
||||
ind1 = sub2ind(M.vnFx, np.c_[ii + posFx, jj])
|
||||
ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posFy]) + M.nFx
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
@@ -459,9 +459,9 @@ def _getFacePxxx_Rectangular(M):
|
||||
posY = 0 if yFace == 'fYm' else 1
|
||||
posZ = 0 if zFace == 'fZm' else 1
|
||||
|
||||
ind1 = sub2ind(M.nFx, np.c_[ii + posX, jj, kk])
|
||||
ind2 = sub2ind(M.nFy, np.c_[ii, jj + posY, kk]) + M.nFv[0]
|
||||
ind3 = sub2ind(M.nFz, np.c_[ii, jj, kk + posZ]) + M.nFv[0] + M.nFv[1]
|
||||
ind1 = sub2ind(M.vnFx, np.c_[ii + posX, jj, kk])
|
||||
ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posY, kk]) + M.nFx
|
||||
ind3 = sub2ind(M.vnFz, np.c_[ii, jj, kk + posZ]) + M.nFx + M.nFy
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
@@ -495,8 +495,8 @@ def _getEdgePxx_Rectangular(M):
|
||||
posX = 0 if xEdge == 'eX0' else 1
|
||||
posY = 0 if yEdge == 'eY0' else 1
|
||||
|
||||
ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX])
|
||||
ind2 = sub2ind(M.nEy, np.c_[ii + posY, jj]) + M.nEv[0]
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
@@ -537,9 +537,9 @@ def _getEdgePxxx_Rectangular(M):
|
||||
posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1]
|
||||
posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1]
|
||||
|
||||
ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX[0], kk + posX[1]])
|
||||
ind2 = sub2ind(M.nEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEv[0]
|
||||
ind3 = sub2ind(M.nEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEv[0] + M.nEv[1]
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX[0], kk + posX[1]])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEx
|
||||
ind3 = sub2ind(M.vnEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEx + M.nEy
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
|
||||
@@ -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))
|
||||
@@ -199,13 +199,13 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
|
||||
def fget(self):
|
||||
if(self._vol is None):
|
||||
if self.dim == 2:
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.nCv+1)
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
|
||||
normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
|
||||
self._vol = area
|
||||
elif self.dim == 3:
|
||||
# Each polyhedron can be decomposed into 5 tetrahedrons
|
||||
# However, this presents a choice so we may as well divide in two ways and average.
|
||||
A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.nCv+1)
|
||||
A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC+1)
|
||||
|
||||
vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top
|
||||
Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top
|
||||
@@ -233,11 +233,11 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 2):
|
||||
xy = self.gridN
|
||||
A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy]))
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
|
||||
edge1 = xy[B, :] - xy[A, :]
|
||||
normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
|
||||
area1 = length2D(edge1)
|
||||
A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy]))
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
|
||||
# Note that we are doing A-D to make sure the normal points the right way.
|
||||
# Think about it. Look at the picture. Normal points towards C iff you do this.
|
||||
edge2 = xy[A, :] - xy[D, :]
|
||||
@@ -247,13 +247,13 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
|
||||
self._normals = [normalize2D(normal1), normalize2D(normal2)]
|
||||
elif(self.dim == 3):
|
||||
|
||||
A, E, F, B = Utils.indexCube('AEFB', self.nCv+1, np.array([self.nNx, self.nCy, self.nCz]))
|
||||
A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz]))
|
||||
normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
|
||||
|
||||
A, D, H, E = Utils.indexCube('ADHE', self.nCv+1, np.array([self.nCx, self.nNy, self.nCz]))
|
||||
A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz]))
|
||||
normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
|
||||
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.nCv+1, np.array([self.nCx, self.nCy, self.nNz]))
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz]))
|
||||
normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
|
||||
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
|
||||
@@ -296,19 +296,19 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
|
||||
if(self._edge is None or self._tangents is None):
|
||||
if(self.dim == 2):
|
||||
xy = self.gridN
|
||||
A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy]))
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
|
||||
edge1 = xy[D, :] - xy[A, :]
|
||||
A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy]))
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
|
||||
edge2 = xy[B, :] - xy[A, :]
|
||||
self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
|
||||
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
|
||||
elif(self.dim == 3):
|
||||
xyz = self.gridN
|
||||
A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy, self.nNz]))
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz]))
|
||||
edge1 = xyz[D, :] - xyz[A, :]
|
||||
A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy, self.nNz]))
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz]))
|
||||
edge2 = xyz[B, :] - xyz[A, :]
|
||||
A, E = Utils.indexCube('AE', self.nCv+1, np.array([self.nNx, self.nNy, self.nCz]))
|
||||
A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz]))
|
||||
edge3 = xyz[E, :] - xyz[A, :]
|
||||
self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
|
||||
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
|
||||
|
||||
@@ -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
|
||||
@@ -282,7 +282,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 1):
|
||||
self._area = np.ones(n[0]+1)
|
||||
@@ -308,7 +308,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.nCv
|
||||
n = self.vnC
|
||||
# Compute edge lengths
|
||||
if(self.dim == 1):
|
||||
self._edge = Utils.mkvc(vh[0])
|
||||
@@ -409,7 +409,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
|
||||
|
||||
ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1
|
||||
if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind:
|
||||
nF_nE = self.nFv if 'F' in locType else self.nEv
|
||||
nF_nE = self.vnF if 'F' in locType else self.vnE
|
||||
components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE]
|
||||
components[ind] = Utils.interpmat(loc, *self.getTensor(locType))
|
||||
Q = sp.hstack(components)
|
||||
|
||||
+21
-21
@@ -61,17 +61,17 @@ class TensorView(object):
|
||||
elif imageType == 'N':
|
||||
assert I.size == self.nN, "Incorrect dimensions for N."
|
||||
elif imageType == 'Fx':
|
||||
if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M')
|
||||
if I.size != np.prod(self.vnFx): I, fy, fz = self.r(I,'F','F','M')
|
||||
elif imageType == 'Fy':
|
||||
if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M')
|
||||
if I.size != np.prod(self.vnFy): fx, I, fz = self.r(I,'F','F','M')
|
||||
elif imageType == 'Fz':
|
||||
if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M')
|
||||
if I.size != np.prod(self.vnFz): fx, fy, I = self.r(I,'F','F','M')
|
||||
elif imageType == 'Ex':
|
||||
if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M')
|
||||
if I.size != np.prod(self.vnEx): I, ey, ez = self.r(I,'E','E','M')
|
||||
elif imageType == 'Ey':
|
||||
if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M')
|
||||
if I.size != np.prod(self.vnEy): ex, I, ez = self.r(I,'E','E','M')
|
||||
elif imageType == 'Ez':
|
||||
if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M')
|
||||
if I.size != np.prod(self.vnEz): ex, ey, I = self.r(I,'E','E','M')
|
||||
elif imageType[0] == 'E':
|
||||
plotAll = len(imageType) == 1
|
||||
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
|
||||
@@ -135,21 +135,21 @@ class TensorView(object):
|
||||
ax.axis('tight')
|
||||
elif self.dim == 2:
|
||||
if imageType == 'CC':
|
||||
C = I[:].reshape(self.nCv, order='F')
|
||||
C = I[:].reshape(self.vnC, order='F')
|
||||
elif imageType == 'N':
|
||||
C = I[:].reshape(self.nNv, order='F')
|
||||
C = I[:].reshape(self.vnN, order='F')
|
||||
C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:])
|
||||
elif imageType == 'Fx':
|
||||
C = I[:].reshape(self.nFx, order='F')
|
||||
C = I[:].reshape(self.vnFx, order='F')
|
||||
C = 0.5*(C[:-1, :] + C[1:, :] )
|
||||
elif imageType == 'Fy':
|
||||
C = I[:].reshape(self.nFy, order='F')
|
||||
C = I[:].reshape(self.vnFy, order='F')
|
||||
C = 0.5*(C[:, :-1] + C[:, 1:] )
|
||||
elif imageType == 'Ex':
|
||||
C = I[:].reshape(self.nEx, order='F')
|
||||
C = I[:].reshape(self.vnEx, order='F')
|
||||
C = 0.5*(C[:,:-1] + C[:,1:] )
|
||||
elif imageType == 'Ey':
|
||||
C = I[:].reshape(self.nEy, order='F')
|
||||
C = I[:].reshape(self.vnEy, order='F')
|
||||
C = 0.5*(C[:-1,:] + C[1:,:] )
|
||||
|
||||
if clim is None:
|
||||
@@ -164,27 +164,27 @@ class TensorView(object):
|
||||
|
||||
# get copy of image and average to cell-centres is necessary
|
||||
if imageType == 'CC':
|
||||
Ic = I[:].reshape(self.nCv, order='F')
|
||||
Ic = I[:].reshape(self.vnC, order='F')
|
||||
elif imageType == 'N':
|
||||
Ic = I[:].reshape(self.nNv, order='F')
|
||||
Ic = I[:].reshape(self.vnN, order='F')
|
||||
Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] )
|
||||
elif imageType == 'Fx':
|
||||
Ic = I[:].reshape(self.nFx, order='F')
|
||||
Ic = I[:].reshape(self.vnFx, order='F')
|
||||
Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:])
|
||||
elif imageType == 'Fy':
|
||||
Ic = I[:].reshape(self.nFy, order='F')
|
||||
Ic = I[:].reshape(self.vnFy, order='F')
|
||||
Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:])
|
||||
elif imageType == 'Fz':
|
||||
Ic = I[:].reshape(self.nFz, order='F')
|
||||
Ic = I[:].reshape(self.vnFz, order='F')
|
||||
Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:])
|
||||
elif imageType == 'Ex':
|
||||
Ic = I[:].reshape(self.nEx, order='F')
|
||||
Ic = I[:].reshape(self.vnEx, order='F')
|
||||
Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1])
|
||||
elif imageType == 'Ey':
|
||||
Ic = I[:].reshape(self.nEy, order='F')
|
||||
Ic = I[:].reshape(self.vnEy, order='F')
|
||||
Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1])
|
||||
elif imageType == 'Ez':
|
||||
Ic = I[:].reshape(self.nEz, order='F')
|
||||
Ic = I[:].reshape(self.vnEz, order='F')
|
||||
Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:])
|
||||
|
||||
# determine number oE slices in x and y dimension
|
||||
@@ -395,7 +395,7 @@ class TensorView(object):
|
||||
mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
|
||||
tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
|
||||
|
||||
return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)])
|
||||
return animate(fig, animateFrame, frames=mesh.vnC['xyz'.index(normal)])
|
||||
|
||||
def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1):
|
||||
"""
|
||||
|
||||
@@ -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
@@ -150,7 +150,7 @@ class Vertical1DModel(BaseModel):
|
||||
|
||||
The number of cells in the
|
||||
last dimension of the mesh."""
|
||||
return self.mesh.nCv[self.mesh.dim-1]
|
||||
return self.mesh.vnC[self.mesh.dim-1]
|
||||
|
||||
def transform(self, m):
|
||||
"""
|
||||
@@ -158,7 +158,7 @@ class Vertical1DModel(BaseModel):
|
||||
:rtype: numpy.array
|
||||
:return: transformed model
|
||||
"""
|
||||
repNum = self.mesh.nCv[:self.mesh.dim-1].prod()
|
||||
repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
|
||||
return Utils.mkvc(m).repeat(repNum)
|
||||
|
||||
def transformDeriv(self, m):
|
||||
@@ -167,7 +167,7 @@ class Vertical1DModel(BaseModel):
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: derivative of transformed model
|
||||
"""
|
||||
repNum = self.mesh.nCv[:self.mesh.dim-1].prod()
|
||||
repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
|
||||
repVec = sp.csr_matrix(
|
||||
(np.ones(repNum),
|
||||
(range(repNum), np.zeros(repNum))
|
||||
|
||||
@@ -207,7 +207,7 @@ class Tikhonov(BaseRegularization):
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFv[0]].T*self.mesh.vol
|
||||
Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
|
||||
self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
|
||||
return self._Wx
|
||||
|
||||
@@ -215,7 +215,7 @@ class Tikhonov(BaseRegularization):
|
||||
def Wy(self):
|
||||
"""Regularization matrix Wy"""
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2])].T*self.mesh.vol
|
||||
Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
|
||||
self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
|
||||
return self._Wy
|
||||
|
||||
@@ -223,7 +223,7 @@ class Tikhonov(BaseRegularization):
|
||||
def Wz(self):
|
||||
"""Regularization matrix Wz"""
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.nFv[:2]):].T*self.mesh.vol
|
||||
Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
|
||||
self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
|
||||
return self._Wz
|
||||
|
||||
|
||||
@@ -44,43 +44,43 @@ class BasicLOMTests(unittest.TestCase):
|
||||
|
||||
def test_tangents(self):
|
||||
T = self.LOM2.tangents
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEv[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEv[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEv[1])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEv[1])))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEx)))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEx)))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEy)))
|
||||
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEy)))
|
||||
|
||||
T = self.LOM3.tangents
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEv[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEv[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEv[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEx)))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEx)))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEx)))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEv[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEv[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEv[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEy)))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEy)))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEy)))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEv[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEv[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEv[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEz)))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEz)))
|
||||
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEz)))
|
||||
|
||||
def test_normals(self):
|
||||
N = self.LOM2.normals
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFv[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFv[0])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFv[1])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFv[1])))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFx)))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFx)))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFy)))
|
||||
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFy)))
|
||||
|
||||
N = self.LOM3.normals
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFv[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFv[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFv[0])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFx)))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFx)))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFx)))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFv[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFv[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFv[1])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFy)))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFy)))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFy)))
|
||||
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFv[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFv[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFv[2])))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFz)))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFz)))
|
||||
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFz)))
|
||||
|
||||
def test_grid(self):
|
||||
self.assertTrue(np.all(self.LOM2.gridCC == self.TM2.gridCC))
|
||||
|
||||
@@ -1,54 +1,46 @@
|
||||
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)
|
||||
|
||||
def test_mesh_nc(self):
|
||||
self.assertTrue(np.all(self.mesh.nCv == [6, 2, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnC == [6, 2, 3]))
|
||||
|
||||
def test_mesh_nc_xyz(self):
|
||||
x = np.all(self.mesh.nCx == 6)
|
||||
y = np.all(self.mesh.nCy == 2)
|
||||
z = np.all(self.mesh.nCz == 3)
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
self.assertTrue(np.all(self.mesh.nCx == 6))
|
||||
self.assertTrue(np.all(self.mesh.nCy == 2))
|
||||
self.assertTrue(np.all(self.mesh.nCz == 3))
|
||||
|
||||
def test_mesh_nf(self):
|
||||
x = np.all(self.mesh.nFx == [7, 2, 3])
|
||||
y = np.all(self.mesh.nFy == [6, 3, 3])
|
||||
z = np.all(self.mesh.nFz == [6, 2, 4])
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [7, 2, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [6, 3, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnFz == [6, 2, 4]))
|
||||
|
||||
def test_mesh_ne(self):
|
||||
x = np.all(self.mesh.nEx == [6, 3, 4])
|
||||
y = np.all(self.mesh.nEy == [7, 2, 4])
|
||||
z = np.all(self.mesh.nEz == [7, 3, 3])
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [6, 3, 4]))
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [7, 2, 4]))
|
||||
self.assertTrue(np.all(self.mesh.vnEz == [7, 3, 3]))
|
||||
|
||||
def test_mesh_numbers(self):
|
||||
c = self.mesh.nC == 36
|
||||
fv = np.all(self.mesh.nFv == [42, 54, 48])
|
||||
ev = np.all(self.mesh.nEv == [72, 56, 63])
|
||||
f = np.all(self.mesh.nF == np.sum([42, 54, 48]))
|
||||
e = np.all(self.mesh.nE == np.sum([72, 56, 63]))
|
||||
|
||||
self.assertTrue(np.all([c, fv, ev, f, e]))
|
||||
self.assertTrue(self.mesh.nC == 36)
|
||||
self.assertTrue(np.all(self.mesh.vnF == [42, 54, 48]))
|
||||
self.assertTrue(np.all(self.mesh.vnE == [72, 56, 63]))
|
||||
self.assertTrue(np.all(self.mesh.nF == np.sum([42, 54, 48])))
|
||||
self.assertTrue(np.all(self.mesh.nE == np.sum([72, 56, 63])))
|
||||
|
||||
def test_mesh_r_E_V(self):
|
||||
ex = np.ones(self.mesh.nEv[0])
|
||||
ey = np.ones(self.mesh.nEv[1])*2
|
||||
ez = np.ones(self.mesh.nEv[2])*3
|
||||
ex = np.ones(self.mesh.nEx)
|
||||
ey = np.ones(self.mesh.nEy)*2
|
||||
ez = np.ones(self.mesh.nEz)*3
|
||||
e = np.r_[ex, ey, ez]
|
||||
tex = self.mesh.r(e, 'E', 'Ex', 'V')
|
||||
tey = self.mesh.r(e, 'E', 'Ey', 'V')
|
||||
@@ -62,9 +54,9 @@ class TestBaseMesh(unittest.TestCase):
|
||||
self.assertTrue(np.all(tez == ez))
|
||||
|
||||
def test_mesh_r_F_V(self):
|
||||
fx = np.ones(self.mesh.nFv[0])
|
||||
fy = np.ones(self.mesh.nFv[1])*2
|
||||
fz = np.ones(self.mesh.nFv[2])*3
|
||||
fx = np.ones(self.mesh.nFx)
|
||||
fy = np.ones(self.mesh.nFy)*2
|
||||
fz = np.ones(self.mesh.nFz)*3
|
||||
f = np.r_[fx, fy, fz]
|
||||
tfx = self.mesh.r(f, 'F', 'Fx', 'V')
|
||||
tfy = self.mesh.r(f, 'F', 'Fy', 'V')
|
||||
@@ -78,25 +70,25 @@ class TestBaseMesh(unittest.TestCase):
|
||||
self.assertTrue(np.all(tfz == fz))
|
||||
|
||||
def test_mesh_r_E_M(self):
|
||||
g = np.ones((np.prod(self.mesh.nEx), 3))
|
||||
g = np.ones((np.prod(self.mesh.vnEx), 3))
|
||||
g[:, 1] = 2
|
||||
g[:, 2] = 3
|
||||
Xex, Yex, Zex = self.mesh.r(g, 'Ex', 'Ex', 'M')
|
||||
self.assertTrue(np.all(Xex.shape == self.mesh.nEx))
|
||||
self.assertTrue(np.all(Yex.shape == self.mesh.nEx))
|
||||
self.assertTrue(np.all(Zex.shape == self.mesh.nEx))
|
||||
self.assertTrue(np.all(Xex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Yex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Zex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Xex == 1))
|
||||
self.assertTrue(np.all(Yex == 2))
|
||||
self.assertTrue(np.all(Zex == 3))
|
||||
|
||||
def test_mesh_r_F_M(self):
|
||||
g = np.ones((np.prod(self.mesh.nFx), 3))
|
||||
g = np.ones((np.prod(self.mesh.vnFx), 3))
|
||||
g[:, 1] = 2
|
||||
g[:, 2] = 3
|
||||
Xfx, Yfx, Zfx = self.mesh.r(g, 'Fx', 'Fx', 'M')
|
||||
self.assertTrue(np.all(Xfx.shape == self.mesh.nFx))
|
||||
self.assertTrue(np.all(Yfx.shape == self.mesh.nFx))
|
||||
self.assertTrue(np.all(Zfx.shape == self.mesh.nFx))
|
||||
self.assertTrue(np.all(Xfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Yfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Zfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Xfx == 1))
|
||||
self.assertTrue(np.all(Yfx == 2))
|
||||
self.assertTrue(np.all(Zfx == 3))
|
||||
@@ -106,9 +98,9 @@ class TestBaseMesh(unittest.TestCase):
|
||||
g[:, 1] = 2
|
||||
g[:, 2] = 3
|
||||
Xc, Yc, Zc = self.mesh.r(g, 'CC', 'CC', 'M')
|
||||
self.assertTrue(np.all(Xc.shape == self.mesh.nCv))
|
||||
self.assertTrue(np.all(Yc.shape == self.mesh.nCv))
|
||||
self.assertTrue(np.all(Zc.shape == self.mesh.nCv))
|
||||
self.assertTrue(np.all(Xc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Yc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Zc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Xc == 1))
|
||||
self.assertTrue(np.all(Yc == 2))
|
||||
self.assertTrue(np.all(Zc == 3))
|
||||
@@ -117,47 +109,44 @@ 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)
|
||||
|
||||
def test_mesh_nc(self):
|
||||
self.assertTrue(np.all(self.mesh.nCv == [6, 2]))
|
||||
self.assertTrue(np.all(self.mesh.vnC == [6, 2]))
|
||||
|
||||
def test_mesh_nc_xyz(self):
|
||||
x = np.all(self.mesh.nCx == 6)
|
||||
y = np.all(self.mesh.nCy == 2)
|
||||
z = self.mesh.nCz is None
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
self.assertTrue(np.all(self.mesh.nCx == 6))
|
||||
self.assertTrue(np.all(self.mesh.nCy == 2))
|
||||
self.assertTrue(self.mesh.nCz is None)
|
||||
|
||||
def test_mesh_nf(self):
|
||||
x = np.all(self.mesh.nFx == [7, 2])
|
||||
y = np.all(self.mesh.nFy == [6, 3])
|
||||
z = self.mesh.nFz is None
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [7, 2]))
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [6, 3]))
|
||||
self.assertTrue(self.mesh.vnFz is None)
|
||||
|
||||
def test_mesh_ne(self):
|
||||
x = np.all(self.mesh.nEx == [6, 3])
|
||||
y = np.all(self.mesh.nEy == [7, 2])
|
||||
z = self.mesh.nEz is None
|
||||
|
||||
self.assertTrue(np.all([x, y, z]))
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [6, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [7, 2]))
|
||||
self.assertTrue(self.mesh.vnEz is None)
|
||||
|
||||
def test_mesh_numbers(self):
|
||||
c = self.mesh.nC == 12
|
||||
fv = np.all(self.mesh.nFv == [14, 18])
|
||||
ev = np.all(self.mesh.nEv == [18, 14])
|
||||
f = np.all(self.mesh.nF == np.sum([14, 18]))
|
||||
e = np.all(self.mesh.nE == np.sum([18, 14]))
|
||||
|
||||
self.assertTrue(np.all([c, fv, ev, f, e]))
|
||||
self.assertTrue(np.all(self.mesh.vnF == [14, 18]))
|
||||
self.assertTrue(np.all(self.mesh.nFx == 14))
|
||||
self.assertTrue(np.all(self.mesh.nFy == 18))
|
||||
self.assertTrue(np.all(self.mesh.nEx == 18))
|
||||
self.assertTrue(np.all(self.mesh.nEy == 14))
|
||||
self.assertTrue(np.all(self.mesh.vnE == [18, 14]))
|
||||
self.assertTrue(np.all(self.mesh.vnE == [18, 14]))
|
||||
self.assertTrue(np.all(self.mesh.nF == np.sum([14, 18])))
|
||||
self.assertTrue(np.all(self.mesh.nE == np.sum([18, 14])))
|
||||
|
||||
def test_mesh_r_E_V(self):
|
||||
ex = np.ones(self.mesh.nEv[0])
|
||||
ey = np.ones(self.mesh.nEv[1])*2
|
||||
ex = np.ones(self.mesh.nEx)
|
||||
ey = np.ones(self.mesh.nEy)*2
|
||||
e = np.r_[ex, ey]
|
||||
tex = self.mesh.r(e, 'E', 'Ex', 'V')
|
||||
tey = self.mesh.r(e, 'E', 'Ey', 'V')
|
||||
@@ -169,8 +158,8 @@ class TestMeshNumbers2D(unittest.TestCase):
|
||||
self.assertRaises(AssertionError, self.mesh.r, e, 'E', 'Ez', 'V')
|
||||
|
||||
def test_mesh_r_F_V(self):
|
||||
fx = np.ones(self.mesh.nFv[0])
|
||||
fy = np.ones(self.mesh.nFv[1])*2
|
||||
fx = np.ones(self.mesh.nFx)
|
||||
fy = np.ones(self.mesh.nFy)*2
|
||||
f = np.r_[fx, fy]
|
||||
tfx = self.mesh.r(f, 'F', 'Fx', 'V')
|
||||
tfy = self.mesh.r(f, 'F', 'Fy', 'V')
|
||||
@@ -182,20 +171,20 @@ class TestMeshNumbers2D(unittest.TestCase):
|
||||
self.assertRaises(AssertionError, self.mesh.r, f, 'F', 'Fz', 'V')
|
||||
|
||||
def test_mesh_r_E_M(self):
|
||||
g = np.ones((np.prod(self.mesh.nEx), 2))
|
||||
g = np.ones((np.prod(self.mesh.vnEx), 2))
|
||||
g[:, 1] = 2
|
||||
Xex, Yex = self.mesh.r(g, 'Ex', 'Ex', 'M')
|
||||
self.assertTrue(np.all(Xex.shape == self.mesh.nEx))
|
||||
self.assertTrue(np.all(Yex.shape == self.mesh.nEx))
|
||||
self.assertTrue(np.all(Xex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Yex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Xex == 1))
|
||||
self.assertTrue(np.all(Yex == 2))
|
||||
|
||||
def test_mesh_r_F_M(self):
|
||||
g = np.ones((np.prod(self.mesh.nFx), 2))
|
||||
g = np.ones((np.prod(self.mesh.vnFx), 2))
|
||||
g[:, 1] = 2
|
||||
Xfx, Yfx = self.mesh.r(g, 'Fx', 'Fx', 'M')
|
||||
self.assertTrue(np.all(Xfx.shape == self.mesh.nFx))
|
||||
self.assertTrue(np.all(Yfx.shape == self.mesh.nFx))
|
||||
self.assertTrue(np.all(Xfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Yfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Xfx == 1))
|
||||
self.assertTrue(np.all(Yfx == 2))
|
||||
|
||||
@@ -203,8 +192,8 @@ class TestMeshNumbers2D(unittest.TestCase):
|
||||
g = np.ones((self.mesh.nC, 2))
|
||||
g[:, 1] = 2
|
||||
Xc, Yc = self.mesh.r(g, 'CC', 'CC', 'M')
|
||||
self.assertTrue(np.all(Xc.shape == self.mesh.nCv))
|
||||
self.assertTrue(np.all(Yc.shape == self.mesh.nCv))
|
||||
self.assertTrue(np.all(Xc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Yc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Xc == 1))
|
||||
self.assertTrue(np.all(Yc == 2))
|
||||
|
||||
|
||||
Reference in New Issue
Block a user