mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 18:52:05 +08:00
Rowan Cockett
+93
-9
@@ -98,8 +98,7 @@ class BaseMesh(object):
|
||||
:rtype: int
|
||||
:return: nEy
|
||||
"""
|
||||
if self.dim < 2: return None
|
||||
return (self._n + np.r_[1,0,1][:self.dim]).prod()
|
||||
return None if self.dim < 2 else (self._n + np.r_[1,0,1][:self.dim]).prod()
|
||||
|
||||
@property
|
||||
def nEz(self):
|
||||
@@ -109,8 +108,7 @@ class BaseMesh(object):
|
||||
:rtype: int
|
||||
:return: nEz
|
||||
"""
|
||||
if self.dim < 3: return None
|
||||
return (self._n + np.r_[1,1,0][:self.dim]).prod()
|
||||
return None if self.dim < 3 else (self._n + np.r_[1,1,0][:self.dim]).prod()
|
||||
|
||||
@property
|
||||
def vnE(self):
|
||||
@@ -157,8 +155,7 @@ class BaseMesh(object):
|
||||
:rtype: int
|
||||
:return: nFy
|
||||
"""
|
||||
if self.dim < 2: return None
|
||||
return (self._n + np.r_[0,1,0][:self.dim]).prod()
|
||||
return None if self.dim < 2 else (self._n + np.r_[0,1,0][:self.dim]).prod()
|
||||
|
||||
@property
|
||||
def nFz(self):
|
||||
@@ -168,8 +165,7 @@ class BaseMesh(object):
|
||||
:rtype: int
|
||||
:return: nFz
|
||||
"""
|
||||
if self.dim < 3: return None
|
||||
return (self._n + np.r_[0,0,1][:self.dim]).prod()
|
||||
return None if self.dim < 3 else (self._n + np.r_[0,0,1][:self.dim]).prod()
|
||||
|
||||
@property
|
||||
def vnF(self):
|
||||
@@ -316,7 +312,7 @@ class BaseRectangularMesh(BaseMesh):
|
||||
@property
|
||||
def nNy(self):
|
||||
"""
|
||||
Number of noes in the y-direction
|
||||
Number of nodes in the y-direction
|
||||
|
||||
:rtype: int
|
||||
:return: nNy or None if dim < 2
|
||||
@@ -403,6 +399,94 @@ class BaseRectangularMesh(BaseMesh):
|
||||
"""
|
||||
return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
|
||||
|
||||
##################################
|
||||
# Redo the numbering so they are dependent of the vector numbers
|
||||
##################################
|
||||
|
||||
@property
|
||||
def nC(self):
|
||||
"""
|
||||
Total number of cells
|
||||
|
||||
:rtype: int
|
||||
:return: nC
|
||||
"""
|
||||
return self.vnC.prod()
|
||||
|
||||
@property
|
||||
def nN(self):
|
||||
"""
|
||||
Total number of nodes
|
||||
|
||||
:rtype: int
|
||||
:return: nN
|
||||
"""
|
||||
return self.vnN.prod()
|
||||
|
||||
@property
|
||||
def nEx(self):
|
||||
"""
|
||||
Number of x-edges
|
||||
|
||||
:rtype: int
|
||||
:return: nEx
|
||||
"""
|
||||
return self.vnEx.prod()
|
||||
|
||||
@property
|
||||
def nEy(self):
|
||||
"""
|
||||
Number of y-edges
|
||||
|
||||
:rtype: int
|
||||
:return: nEy
|
||||
"""
|
||||
if self.dim < 2: return
|
||||
return self.vnEy.prod()
|
||||
|
||||
@property
|
||||
def nEz(self):
|
||||
"""
|
||||
Number of z-edges
|
||||
|
||||
:rtype: int
|
||||
:return: nEz
|
||||
"""
|
||||
if self.dim < 3: return
|
||||
return self.vnEz.prod()
|
||||
|
||||
@property
|
||||
def nFx(self):
|
||||
"""
|
||||
Number of x-faces
|
||||
|
||||
:rtype: int
|
||||
:return: nFx
|
||||
"""
|
||||
return self.vnFx.prod()
|
||||
|
||||
@property
|
||||
def nFy(self):
|
||||
"""
|
||||
Number of y-faces
|
||||
|
||||
:rtype: int
|
||||
:return: nFy
|
||||
"""
|
||||
if self.dim < 2: return
|
||||
return self.vnFy.prod()
|
||||
|
||||
@property
|
||||
def nFz(self):
|
||||
"""
|
||||
Number of z-faces
|
||||
|
||||
:rtype: int
|
||||
:return: nFz
|
||||
"""
|
||||
if self.dim < 3: return
|
||||
return self.vnFz.prod()
|
||||
|
||||
def r(self, x, xType='CC', outType='CC', format='V'):
|
||||
"""
|
||||
Mesh.r is a quick reshape command that will do the best it can at giving you what you want.
|
||||
|
||||
@@ -12,6 +12,7 @@ class Cyl1DMesh(object):
|
||||
|
||||
def __init__(self, h, z0=None):
|
||||
assert len(h) == 2, "len(h) must equal 2"
|
||||
print 'PendingDeprecationWarning: Cyl1DMesh has been replaced by CylMesh. hy == theta == 2pi, use one cell to make this cylindrically symmetric.'
|
||||
if z0 is not None:
|
||||
assert z0.size == 1, "z0.size must equal 1"
|
||||
|
||||
|
||||
@@ -0,0 +1,293 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from scipy.constants import pi
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap
|
||||
from TensorMesh import BaseTensorMesh
|
||||
from InnerProducts import InnerProducts
|
||||
|
||||
|
||||
class CylMesh(BaseTensorMesh, InnerProducts):
|
||||
"""
|
||||
CylMesh is a mesh class for cylindrical problems
|
||||
|
||||
::
|
||||
|
||||
cs, nc, npad = 20., 30, 8
|
||||
hx = Utils.meshTensors(((npad+10,cs,0.7), (nc,cs), (npad,cs)))
|
||||
hz = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs)))
|
||||
mesh = Mesh.CylMesh([hx,1,hz], [0.,0,-hz.sum()/2.])
|
||||
"""
|
||||
|
||||
_meshType = 'CYL'
|
||||
|
||||
_unitDimensions = [1, 2*np.pi, 1]
|
||||
|
||||
def __init__(self, h, x0=None):
|
||||
BaseTensorMesh.__init__(self, h, x0)
|
||||
assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]"
|
||||
assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
|
||||
|
||||
@property
|
||||
def isSymmetric(self):
|
||||
return self.nCy == 1
|
||||
|
||||
@property
|
||||
def nNx(self):
|
||||
"""
|
||||
Number of nodes in the x-direction
|
||||
|
||||
:rtype: int
|
||||
:return: nNx
|
||||
"""
|
||||
if self.isSymmetric: return self.nCx
|
||||
return self.nCx + 1
|
||||
|
||||
@property
|
||||
def nNy(self):
|
||||
"""
|
||||
Number of nodes in the y-direction
|
||||
|
||||
:rtype: int
|
||||
:return: nNy
|
||||
"""
|
||||
if self.isSymmetric: return 0
|
||||
return self.nCy
|
||||
|
||||
@property
|
||||
def vnFx(self):
|
||||
"""
|
||||
Number of x-faces in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: vnFx
|
||||
"""
|
||||
return self.vnC
|
||||
|
||||
@property
|
||||
def vnEy(self):
|
||||
"""
|
||||
Number of y-edges in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: vnEy or None if dim < 2
|
||||
"""
|
||||
nNx = self.nNx if self.isSymmetric else self.nNx - 1
|
||||
return np.r_[nNx, self.nCy, self.nNz]
|
||||
|
||||
@property
|
||||
def vnEz(self):
|
||||
"""
|
||||
Number of z-edges in each direction
|
||||
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: vnEz or None if nCy > 1
|
||||
"""
|
||||
if self.isSymmetric:
|
||||
return np.r_[self.nNx, self.nNy, self.nCz]
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def nEz(self):
|
||||
"""
|
||||
Number of z-edges
|
||||
|
||||
:rtype: int
|
||||
:return: nEz
|
||||
"""
|
||||
if self.isSymmetric:
|
||||
return self.vnEz.prod()
|
||||
return (np.r_[self.nNx-1, self.nNy, self.nCz]).prod() + self.nCz
|
||||
|
||||
@property
|
||||
def vectorCCx(self):
|
||||
"""Cell-centered grid vector (1D) in the x direction."""
|
||||
return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5
|
||||
|
||||
@property
|
||||
def vectorCCy(self):
|
||||
"""Cell-centered grid vector (1D) in the y direction."""
|
||||
return np.r_[0, self.hy[:-1]]
|
||||
|
||||
@property
|
||||
def vectorNx(self):
|
||||
"""Nodal grid vector (1D) in the x direction."""
|
||||
if self.isSymmetric:
|
||||
return self.hx.cumsum()
|
||||
return np.r_[0, self.hx].cumsum()
|
||||
|
||||
@property
|
||||
def vectorNy(self):
|
||||
"""Nodal grid vector (1D) in the y direction."""
|
||||
if self.isSymmetric:
|
||||
# There aren't really any nodes, but all the grids need
|
||||
# somewhere to live, why not zero?!
|
||||
return np.r_[0]
|
||||
return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5
|
||||
|
||||
@property
|
||||
def edge(self):
|
||||
"""Edge lengths"""
|
||||
if getattr(self, '_edge', None) is None:
|
||||
if self.isSymmetric:
|
||||
self._edge = 2*pi*self.gridN[:,0]
|
||||
else:
|
||||
raise NotImplementedError('edges not yet implemented for 3D cyl mesh')
|
||||
return self._edge
|
||||
|
||||
@property
|
||||
def area(self):
|
||||
"""Face areas"""
|
||||
if getattr(self, '_area', None) is None:
|
||||
if self.nCy > 1:
|
||||
raise NotImplementedError('area not yet implemented for 3D cyl mesh')
|
||||
areaR = np.kron(self.hz, 2*pi*self.vectorNx)
|
||||
areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2))
|
||||
self._area = np.r_[areaR, areaZ]
|
||||
return self._area
|
||||
|
||||
@property
|
||||
def vol(self):
|
||||
"""Volume of each cell"""
|
||||
if getattr(self, '_vol', None) is None:
|
||||
if self.nCy > 1:
|
||||
raise NotImplementedError('vol not yet implemented for 3D cyl mesh')
|
||||
az = pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2)
|
||||
self._vol = np.kron(self.hz, az)
|
||||
return self._vol
|
||||
|
||||
####################################################
|
||||
# Operators
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def faceDiv(self):
|
||||
"""Construct divergence operator (face-stg to cell-centres)."""
|
||||
if getattr(self, '_faceDiv', None) is None:
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
D1 = self.faceDivx
|
||||
D3 = self.faceDivz
|
||||
if self.isSymmetric:
|
||||
D = sp.hstack((D1, D3), format="csr")
|
||||
elif self.nCy > 1:
|
||||
D2 = self.faceDivy
|
||||
D = sp.hstack((D1, D2, D3), format="csr")
|
||||
self._faceDiv = D
|
||||
return self._faceDiv
|
||||
|
||||
@property
|
||||
def faceDivx(self):
|
||||
"""Construct divergence operator in the x component (face-stg to cell-centres)."""
|
||||
if getattr(self, '_faceDivx', None) is None:
|
||||
D1 = kron3(speye(self.nCz), speye(self.nCy), ddx(self.nCx)[:,1:])
|
||||
S = self.r(self.area, 'F', 'Fx', 'V')
|
||||
V = self.vol
|
||||
self._faceDivx = sdiag(1/V)*D1*sdiag(S)
|
||||
return self._faceDivx
|
||||
|
||||
@property
|
||||
def faceDivy(self):
|
||||
"""Construct divergence operator in the y component (face-stg to cell-centres)."""
|
||||
raise NotImplementedError('Wrapping the ddx is not yet implemented.')
|
||||
if getattr(self, '_faceDivy', None) is None:
|
||||
# TODO: this needs to wrap to join up faces which are connected in the cylinder
|
||||
D2 = kron3(speye(self.nCz), ddx(self.nCy), speye(self.nCx))
|
||||
S = self.r(self.area, 'F', 'Fy', 'V')
|
||||
V = self.vol
|
||||
self._faceDivy = sdiag(1/V)*D2*sdiag(S)
|
||||
return self._faceDivy
|
||||
|
||||
@property
|
||||
def faceDivz(self):
|
||||
"""Construct divergence operator in the z component (face-stg to cell-centres)."""
|
||||
if getattr(self, '_faceDivz', None) is None:
|
||||
D3 = kron3(ddx(self.nCz), speye(self.nCy), speye(self.nCx))
|
||||
S = self.r(self.area, 'F', 'Fz', 'V')
|
||||
V = self.vol
|
||||
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
|
||||
return self._faceDivz
|
||||
|
||||
|
||||
@property
|
||||
def cellGrad(self):
|
||||
"""The cell centered Gradient, takes you to cell faces."""
|
||||
raise NotImplementedError('Cell Grad is not yet implemented.')
|
||||
|
||||
|
||||
@property
|
||||
def nodalGrad(self):
|
||||
"""Construct gradient operator (nodes to edges)."""
|
||||
# Nodal grad does not make sense for cylindrically symmetric mesh.
|
||||
if self.isSymmetric: return None
|
||||
raise NotImplementedError('nodalGrad not yet implemented')
|
||||
|
||||
@property
|
||||
def nodalLaplacian(self):
|
||||
"""Construct laplacian operator (nodes to edges)."""
|
||||
raise NotImplementedError('nodalLaplacian not yet implemented')
|
||||
|
||||
@property
|
||||
def edgeCurl(self):
|
||||
"""The edgeCurl property."""
|
||||
if self.nCy > 1:
|
||||
raise NotImplementedError('Edge curl not yet implemented for nCy > 1')
|
||||
if getattr(self, '_edgeCurl', None) is None:
|
||||
#1D Difference matricies
|
||||
dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1,0], self.nCx, self.nCx, format="csr")
|
||||
dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0,1], self.nCz, self.nCz+1, format="csr")
|
||||
|
||||
#2D Difference matricies
|
||||
Dr = sp.kron(sp.eye(self.nNz), dr)
|
||||
Dz = -sp.kron(dz, sp.eye(self.nCx)) #Not sure about this negative
|
||||
|
||||
#Edge curl operator
|
||||
self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0)
|
||||
return self._edgeCurl
|
||||
|
||||
# @property
|
||||
# def aveE2CC(self):
|
||||
# """Averaging operator from cell edges to cell centers"""
|
||||
# if getattr(self, '_aveE2CC', None) is None:
|
||||
# if self.isSymmetric:
|
||||
# az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr')
|
||||
# ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr')
|
||||
# ar[0,0] = 1
|
||||
# self._aveE2CC = (0.5)*sp.kron(az, ar).T
|
||||
# else:
|
||||
# raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
# return self._aveE2CC
|
||||
|
||||
@property
|
||||
def aveE2CC(self):
|
||||
"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.vnC
|
||||
if self.isSymmetric:
|
||||
avR = av(n[0])[:,1:]
|
||||
avR[0,0] = 1.
|
||||
self._aveE2CC = (0.5)*sp.kron(av(n[2]), avR, format="csr")
|
||||
else:
|
||||
raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
# self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
# kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
# kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveE2CC
|
||||
|
||||
|
||||
@property
|
||||
def aveF2CC(self):
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
if getattr(self, '_aveF2CC', None) is None:
|
||||
n = self.vnC
|
||||
if self.isSymmetric:
|
||||
avR = av(n[0])[:,1:]
|
||||
avR[0,0] = 1.
|
||||
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[2]), avR), sp.kron(av(n[2]), speye(n[0]))), format="csr")
|
||||
else:
|
||||
raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
# self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
|
||||
# kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
# kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._aveF2CC
|
||||
@@ -167,7 +167,7 @@ class DiffOperators(object):
|
||||
elif(self.dim == 3):
|
||||
D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.r(self.area, 'F','Fx', 'V')
|
||||
S = self.r(self.area, 'F', 'Fx', 'V')
|
||||
V = self.vol
|
||||
self._faceDivx = sdiag(1/V)*D1*sdiag(S)
|
||||
|
||||
@@ -190,7 +190,7 @@ class DiffOperators(object):
|
||||
elif(self.dim == 3):
|
||||
D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.r(self.area, 'F','Fy', 'V')
|
||||
S = self.r(self.area, 'F', 'Fy', 'V')
|
||||
V = self.vol
|
||||
self._faceDivy = sdiag(1/V)*D2*sdiag(S)
|
||||
|
||||
@@ -210,7 +210,7 @@ class DiffOperators(object):
|
||||
# Compute faceDivergence operator on faces
|
||||
D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.r(self.area, 'F','Fz', 'V')
|
||||
S = self.r(self.area, 'F', 'Fz', 'V')
|
||||
V = self.vol
|
||||
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
|
||||
|
||||
@@ -623,11 +623,11 @@ class DiffOperators(object):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
elif(self.dim == 2):
|
||||
self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveE2CC
|
||||
|
||||
@property
|
||||
@@ -696,57 +696,3 @@ class DiffOperators(object):
|
||||
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveN2F
|
||||
|
||||
# --------------- Methods ---------------------
|
||||
|
||||
def getMass(self, materialProp=None, loc='e'):
|
||||
""" Produces mass matricies.
|
||||
|
||||
:param str loc: Average to location: 'e'-edges, 'f'-faces
|
||||
:param None,float,numpy.ndarray materialProp: property to be averaged (see below)
|
||||
:rtype: scipy.sparse.csr.csr_matrix
|
||||
:return: M, the mass matrix
|
||||
|
||||
materialProp can be::
|
||||
|
||||
None -> takes materialProp = 1 (default)
|
||||
float -> a constant value for entire domain
|
||||
numpy.ndarray -> if materialProp.size == self.nC
|
||||
3D property model
|
||||
if materialProp.size = self.nCz
|
||||
1D (layered eath) property model
|
||||
"""
|
||||
if materialProp is None:
|
||||
materialProp = np.ones(self.nC)
|
||||
elif type(materialProp) is float:
|
||||
materialProp = np.ones(self.nC)*materialProp
|
||||
elif materialProp.shape == (self.nCz,):
|
||||
materialProp = materialProp.repeat(self.nCx*self.nCy)
|
||||
materialProp = mkvc(materialProp)
|
||||
assert materialProp.shape == (self.nC,), "materialProp incorrect shape"
|
||||
|
||||
if loc=='e':
|
||||
Av = self.aveE2CC
|
||||
elif loc=='f':
|
||||
Av = self.aveF2CC
|
||||
else:
|
||||
raise ValueError('Invalid loc')
|
||||
|
||||
diag = Av.T * (self.vol * mkvc(materialProp))
|
||||
|
||||
return sdiag(diag)
|
||||
|
||||
def getEdgeMass(self, materialProp=None):
|
||||
"""mass matrix for products of edge functions w'*M(materialProp)*e"""
|
||||
return self.getMass(loc='e', materialProp=materialProp)
|
||||
|
||||
def getFaceMass(self, materialProp=None):
|
||||
"""mass matrix for products of face functions w'*M(materialProp)*f"""
|
||||
return self.getMass(loc='f', materialProp=materialProp)
|
||||
|
||||
def getFaceMassDeriv(self):
|
||||
Av = self.aveF2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
def getEdgeMassDeriv(self):
|
||||
Av = self.aveE2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
@@ -284,7 +284,6 @@ class InnerProducts(object):
|
||||
# | |/
|
||||
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
|
||||
def _getFacePx(M):
|
||||
assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
|
||||
return _getFacePx_Rectangular(M)
|
||||
|
||||
+302
-337
@@ -1,50 +1,31 @@
|
||||
from SimPEG import Utils, np, sp
|
||||
from BaseMesh import BaseRectangularMesh
|
||||
from TensorView import TensorView
|
||||
from View import TensorView
|
||||
from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
|
||||
class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
"""
|
||||
TensorMesh is a mesh class that deals with tensor product meshes.
|
||||
def _getTensorGrid(self, key):
|
||||
if getattr(self, '_grid' + key, None) is None:
|
||||
setattr(self, '_grid' + key, Utils.ndgrid(self.getTensor(key)))
|
||||
return getattr(self, '_grid' + key)
|
||||
|
||||
Any Mesh that has a constant width along the entire axis
|
||||
such that it can defined by a single width vector, called 'h'.
|
||||
|
||||
::
|
||||
|
||||
hx = np.array([1,1,1])
|
||||
hy = np.array([1,2])
|
||||
hz = np.array([1,1,1,1])
|
||||
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz])
|
||||
|
||||
Example of a padded tensor mesh:
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Mesh, Utils
|
||||
M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
|
||||
M.plotGrid()
|
||||
|
||||
For a quick tensor mesh on a (10x12x15) unit cube::
|
||||
|
||||
mesh = Mesh.TensorMesh([10, 12, 15])
|
||||
|
||||
"""
|
||||
class BaseTensorMesh(BaseRectangularMesh):
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
_meshType = 'TENSOR'
|
||||
_meshType = 'BASETENSOR'
|
||||
|
||||
_unitDimensions = [1, 1, 1]
|
||||
|
||||
def __init__(self, h_in, x0=None):
|
||||
assert type(h_in) is list, 'h_in must be a list'
|
||||
assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3'
|
||||
h = range(len(h_in))
|
||||
for i, h_i in enumerate(h_in):
|
||||
if type(h_i) in [int, long, float]:
|
||||
if type(h_i) in [int, long, float, np.int_]:
|
||||
# This gives you something over the unit cube.
|
||||
h_i = np.ones(int(h_i))/int(h_i)
|
||||
h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i)
|
||||
assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
|
||||
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
|
||||
h[i] = h_i[:] # make a copy.
|
||||
@@ -55,279 +36,101 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
# Ensure h contains 1D vectors
|
||||
self._h = [Utils.mkvc(x.astype(float)) for x in h]
|
||||
|
||||
def __str__(self):
|
||||
outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim)
|
||||
def printH(hx, outStr=''):
|
||||
i = -1
|
||||
while True:
|
||||
i = i + 1
|
||||
if i > hx.size:
|
||||
break
|
||||
elif i == hx.size:
|
||||
break
|
||||
h = hx[i]
|
||||
n = 1
|
||||
for j in range(i+1, hx.size):
|
||||
if hx[j] == h:
|
||||
n = n + 1
|
||||
i = i + 1
|
||||
else:
|
||||
break
|
||||
@property
|
||||
def h(self):
|
||||
"""h is a list containing the cell widths of the tensor mesh in each dimension."""
|
||||
return self._h
|
||||
|
||||
if n == 1:
|
||||
outStr = outStr + ' {0:.2f},'.format(h)
|
||||
else:
|
||||
outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h)
|
||||
@property
|
||||
def hx(self):
|
||||
"Width of cells in the x direction"
|
||||
return self._h[0]
|
||||
|
||||
return outStr[:-1]
|
||||
@property
|
||||
def hy(self):
|
||||
"Width of cells in the y direction"
|
||||
return None if self.dim < 2 else self._h[1]
|
||||
|
||||
if self.dim == 1:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
pass
|
||||
elif self.dim == 2:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr = outStr + printH(self.hy, outStr='\n hy:')
|
||||
elif self.dim == 3:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr = outStr + '\n nCz: {0:d}'.format(self.nCz)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr = outStr + printH(self.hy, outStr='\n hy:')
|
||||
outStr = outStr + printH(self.hz, outStr='\n hz:')
|
||||
@property
|
||||
def hz(self):
|
||||
"Width of cells in the z direction"
|
||||
return None if self.dim < 3 else self._h[2]
|
||||
|
||||
return outStr
|
||||
@property
|
||||
def vectorNx(self):
|
||||
"""Nodal grid vector (1D) in the x direction."""
|
||||
return np.r_[0., self.hx.cumsum()] + self.x0[0]
|
||||
|
||||
def h():
|
||||
doc = "h is a list containing the cell widths of the tensor mesh in each dimension."
|
||||
fget = lambda self: self._h
|
||||
return locals()
|
||||
h = property(**h())
|
||||
@property
|
||||
def vectorNy(self):
|
||||
"""Nodal grid vector (1D) in the y direction."""
|
||||
return None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1]
|
||||
|
||||
def hx():
|
||||
doc = "Width of cells in the x direction"
|
||||
fget = lambda self: self._h[0]
|
||||
return locals()
|
||||
hx = property(**hx())
|
||||
@property
|
||||
def vectorNz(self):
|
||||
"""Nodal grid vector (1D) in the z direction."""
|
||||
return None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2]
|
||||
|
||||
def hy():
|
||||
doc = "Width of cells in the y direction"
|
||||
fget = lambda self: None if self.dim < 2 else self._h[1]
|
||||
return locals()
|
||||
hy = property(**hy())
|
||||
@property
|
||||
def vectorCCx(self):
|
||||
"""Cell-centered grid vector (1D) in the x direction."""
|
||||
return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0]
|
||||
|
||||
def hz():
|
||||
doc = "Width of cells in the z direction"
|
||||
fget = lambda self: None if self.dim < 3 else self._h[2]
|
||||
return locals()
|
||||
hz = property(**hz())
|
||||
@property
|
||||
def vectorCCy(self):
|
||||
"""Cell-centered grid vector (1D) in the y direction."""
|
||||
return None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1]
|
||||
|
||||
def vectorNx():
|
||||
doc = "Nodal grid vector (1D) in the x direction."
|
||||
fget = lambda self: np.r_[0., self.hx.cumsum()] + self.x0[0]
|
||||
return locals()
|
||||
vectorNx = property(**vectorNx())
|
||||
@property
|
||||
def vectorCCz(self):
|
||||
"""Cell-centered grid vector (1D) in the z direction."""
|
||||
return None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2]
|
||||
|
||||
def vectorNy():
|
||||
doc = "Nodal grid vector (1D) in the y direction."
|
||||
fget = lambda self: None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1]
|
||||
return locals()
|
||||
vectorNy = property(**vectorNy())
|
||||
@property
|
||||
def gridCC(self):
|
||||
"""Cell-centered grid."""
|
||||
return _getTensorGrid(self, 'CC')
|
||||
|
||||
def vectorNz():
|
||||
doc = "Nodal grid vector (1D) in the z direction."
|
||||
fget = lambda self: None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2]
|
||||
return locals()
|
||||
vectorNz = property(**vectorNz())
|
||||
@property
|
||||
def gridN(self):
|
||||
"""Nodal grid."""
|
||||
return _getTensorGrid(self, 'N')
|
||||
|
||||
def vectorCCx():
|
||||
doc = "Cell-centered grid vector (1D) in the x direction."
|
||||
fget = lambda self: np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0]
|
||||
return locals()
|
||||
vectorCCx = property(**vectorCCx())
|
||||
@property
|
||||
def gridFx(self):
|
||||
"""Face staggered grid in the x direction."""
|
||||
if self.nFx == 0: return
|
||||
return _getTensorGrid(self, 'Fx')
|
||||
|
||||
def vectorCCy():
|
||||
doc = "Cell-centered grid vector (1D) in the y direction."
|
||||
fget = lambda self: None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1]
|
||||
return locals()
|
||||
vectorCCy = property(**vectorCCy())
|
||||
@property
|
||||
def gridFy(self):
|
||||
"""Face staggered grid in the y direction."""
|
||||
if self.nFy == 0 or self.dim < 2: return
|
||||
return _getTensorGrid(self, 'Fy')
|
||||
|
||||
def vectorCCz():
|
||||
doc = "Cell-centered grid vector (1D) in the z direction."
|
||||
fget = lambda self: None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2]
|
||||
return locals()
|
||||
vectorCCz = property(**vectorCCz())
|
||||
@property
|
||||
def gridFz(self):
|
||||
"""Face staggered grid in the z direction."""
|
||||
if self.nFz == 0 or self.dim < 3: return
|
||||
return _getTensorGrid(self, 'Fz')
|
||||
|
||||
def gridCC():
|
||||
doc = "Cell-centered grid."
|
||||
@property
|
||||
def gridEx(self):
|
||||
"""Edge staggered grid in the x direction."""
|
||||
if self.nEx == 0: return
|
||||
return _getTensorGrid(self, 'Ex')
|
||||
|
||||
def fget(self):
|
||||
if self._gridCC is None:
|
||||
self._gridCC = Utils.ndgrid(self.getTensor('CC'))
|
||||
return self._gridCC
|
||||
return locals()
|
||||
_gridCC = None # Store grid by default
|
||||
gridCC = property(**gridCC())
|
||||
@property
|
||||
def gridEy(self):
|
||||
"""Edge staggered grid in the y direction."""
|
||||
if self.nEy == 0 or self.dim < 2: return
|
||||
return _getTensorGrid(self, 'Ey')
|
||||
|
||||
def gridN():
|
||||
doc = "Nodal grid."
|
||||
|
||||
def fget(self):
|
||||
if self._gridN is None:
|
||||
self._gridN = Utils.ndgrid(self.getTensor('N'))
|
||||
return self._gridN
|
||||
return locals()
|
||||
_gridN = None # Store grid by default
|
||||
gridN = property(**gridN())
|
||||
|
||||
def gridFx():
|
||||
doc = "Face staggered grid in the x direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFx is None:
|
||||
self._gridFx = Utils.ndgrid(self.getTensor('Fx'))
|
||||
return self._gridFx
|
||||
return locals()
|
||||
_gridFx = None # Store grid by default
|
||||
gridFx = property(**gridFx())
|
||||
|
||||
def gridFy():
|
||||
doc = "Face staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFy is None and self.dim > 1:
|
||||
self._gridFy = Utils.ndgrid(self.getTensor('Fy'))
|
||||
return self._gridFy
|
||||
return locals()
|
||||
_gridFy = None # Store grid by default
|
||||
gridFy = property(**gridFy())
|
||||
|
||||
def gridFz():
|
||||
doc = "Face staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridFz is None and self.dim > 2:
|
||||
self._gridFz = Utils.ndgrid(self.getTensor('Fz'))
|
||||
return self._gridFz
|
||||
return locals()
|
||||
_gridFz = None # Store grid by default
|
||||
gridFz = property(**gridFz())
|
||||
|
||||
def gridEx():
|
||||
doc = "Edge staggered grid in the x direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEx is None:
|
||||
self._gridEx = Utils.ndgrid(self.getTensor('Ex'))
|
||||
return self._gridEx
|
||||
return locals()
|
||||
_gridEx = None # Store grid by default
|
||||
gridEx = property(**gridEx())
|
||||
|
||||
def gridEy():
|
||||
doc = "Edge staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEy is None and self.dim > 1:
|
||||
self._gridEy = Utils.ndgrid(self.getTensor('Ey'))
|
||||
return self._gridEy
|
||||
return locals()
|
||||
_gridEy = None # Store grid by default
|
||||
gridEy = property(**gridEy())
|
||||
|
||||
def gridEz():
|
||||
doc = "Edge staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEz is None and self.dim > 2:
|
||||
self._gridEz = Utils.ndgrid(self.getTensor('Ez'))
|
||||
return self._gridEz
|
||||
return locals()
|
||||
_gridEz = None # Store grid by default
|
||||
gridEz = property(**gridEz())
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
def vol():
|
||||
doc = "Construct cell volumes of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._vol is None):
|
||||
vh = self.h
|
||||
# Compute cell volumes
|
||||
if(self.dim == 1):
|
||||
self._vol = Utils.mkvc(vh[0])
|
||||
elif(self.dim == 2):
|
||||
# Cell sizes in each direction
|
||||
self._vol = Utils.mkvc(np.outer(vh[0], vh[1]))
|
||||
elif(self.dim == 3):
|
||||
# Cell sizes in each direction
|
||||
self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2]))
|
||||
return self._vol
|
||||
return locals()
|
||||
_vol = None
|
||||
vol = property(**vol())
|
||||
|
||||
def area():
|
||||
doc = "Construct face areas of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._area is None):
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 1):
|
||||
self._area = np.ones(n[0]+1)
|
||||
elif(self.dim == 2):
|
||||
area1 = np.outer(np.ones(n[0]+1), vh[1])
|
||||
area2 = np.outer(vh[0], np.ones(n[1]+1))
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
|
||||
elif(self.dim == 3):
|
||||
area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2])))
|
||||
area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
|
||||
area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
|
||||
return self._area
|
||||
return locals()
|
||||
_area = None
|
||||
area = property(**area())
|
||||
|
||||
def edge():
|
||||
doc = "Construct edge legnths of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._edge is None):
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute edge lengths
|
||||
if(self.dim == 1):
|
||||
self._edge = Utils.mkvc(vh[0])
|
||||
elif(self.dim == 2):
|
||||
l1 = np.outer(vh[0], np.ones(n[1]+1))
|
||||
l2 = np.outer(np.ones(n[0]+1), vh[1])
|
||||
self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)]
|
||||
elif(self.dim == 3):
|
||||
l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
|
||||
l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
|
||||
l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
|
||||
self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)]
|
||||
return self._edge
|
||||
return locals()
|
||||
_edge = None
|
||||
edge = property(**edge())
|
||||
|
||||
# --------------- Methods ---------------------
|
||||
@property
|
||||
def gridEz(self):
|
||||
"""Edge staggered grid in the z direction."""
|
||||
if self.nEz == 0 or self.dim < 3: return
|
||||
return _getTensorGrid(self, 'Ez')
|
||||
|
||||
def getTensor(self, locType):
|
||||
""" Returns a tensor list.
|
||||
@@ -367,6 +170,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
|
||||
return [t for t in ten if t is not None]
|
||||
|
||||
# --------------- Methods ---------------------
|
||||
|
||||
def isInside(self, pts, locType='N'):
|
||||
"""
|
||||
@@ -429,69 +233,26 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
indZeros = np.logical_not(self.isInside(loc))
|
||||
loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')])
|
||||
|
||||
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:
|
||||
if locType in ['Fx','Fy','Fz','Ex','Ey','Ez']:
|
||||
ind = {'x':0, 'y':1, 'z':2}[locType[1]]
|
||||
assert self.dim >= ind, 'mesh is not high enough dimension.'
|
||||
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))
|
||||
# remove any zero blocks (hstack complains)
|
||||
components = [comp for comp in components if comp.shape[1] > 0]
|
||||
Q = sp.hstack(components)
|
||||
elif locType in ['CC', 'N']:
|
||||
Q = Utils.interpmat(loc, *self.getTensor(locType))
|
||||
else:
|
||||
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
|
||||
|
||||
if zerosOutside:
|
||||
Q[indZeros, :] = 0
|
||||
|
||||
return Q.tocsr()
|
||||
|
||||
|
||||
@property
|
||||
def faceBoundaryInd(self):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
"""
|
||||
if self.dim==1:
|
||||
indxd = (self.gridFx==min(self.gridFx))
|
||||
indxu = (self.gridFx==max(self.gridFx))
|
||||
return indxd, indxu
|
||||
elif self.dim==2:
|
||||
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
|
||||
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
|
||||
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
|
||||
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif self.dim==3:
|
||||
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
|
||||
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
|
||||
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
|
||||
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
|
||||
indzd = (self.gridFz[:,2]==min(self.gridFz[:,2]))
|
||||
indzu = (self.gridFz[:,2]==max(self.gridFz[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
@property
|
||||
def cellBoundaryInd(self):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
"""
|
||||
if self.dim==1:
|
||||
indxd = (self.gridCC==min(self.gridCC))
|
||||
indxu = (self.gridCC==max(self.gridCC))
|
||||
return indxd, indxu
|
||||
elif self.dim==2:
|
||||
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
|
||||
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
|
||||
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
|
||||
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif self.dim==3:
|
||||
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
|
||||
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
|
||||
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
|
||||
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
|
||||
indzd = (self.gridCC[:,2]==min(self.gridCC[:,2]))
|
||||
indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
@@ -577,7 +338,8 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
if materialProperty is None:
|
||||
return None
|
||||
return 0.0
|
||||
|
||||
if Utils.isScalar(materialProperty):
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
V = Utils.sdiag(self.vol)
|
||||
@@ -585,12 +347,14 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
if v is None:
|
||||
return self.dim * Av.T * V * ones
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V * ones
|
||||
|
||||
if materialProperty.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
V = Utils.sdiag(self.vol)
|
||||
if v is None:
|
||||
return self.dim * Av.T * V
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V
|
||||
|
||||
if materialProperty.size == self.nC*self.dim: # anisotropic
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
|
||||
@@ -599,6 +363,207 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
return Utils.sdiag(v) * Av.T * V
|
||||
|
||||
|
||||
|
||||
class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts):
|
||||
"""
|
||||
TensorMesh is a mesh class that deals with tensor product meshes.
|
||||
|
||||
Any Mesh that has a constant width along the entire axis
|
||||
such that it can defined by a single width vector, called 'h'.
|
||||
|
||||
::
|
||||
|
||||
hx = np.array([1,1,1])
|
||||
hy = np.array([1,2])
|
||||
hz = np.array([1,1,1,1])
|
||||
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz])
|
||||
|
||||
Example of a padded tensor mesh:
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Mesh, Utils
|
||||
M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
|
||||
M.plotGrid()
|
||||
|
||||
For a quick tensor mesh on a (10x12x15) unit cube::
|
||||
|
||||
mesh = Mesh.TensorMesh([10, 12, 15])
|
||||
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
_meshType = 'TENSOR'
|
||||
|
||||
def __init__(self, h_in, x0=None):
|
||||
BaseTensorMesh.__init__(self, h_in, x0)
|
||||
|
||||
def __str__(self):
|
||||
outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim)
|
||||
def printH(hx, outStr=''):
|
||||
i = -1
|
||||
while True:
|
||||
i = i + 1
|
||||
if i > hx.size:
|
||||
break
|
||||
elif i == hx.size:
|
||||
break
|
||||
h = hx[i]
|
||||
n = 1
|
||||
for j in range(i+1, hx.size):
|
||||
if hx[j] == h:
|
||||
n = n + 1
|
||||
i = i + 1
|
||||
else:
|
||||
break
|
||||
|
||||
if n == 1:
|
||||
outStr = outStr + ' {0:.2f},'.format(h)
|
||||
else:
|
||||
outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h)
|
||||
|
||||
return outStr[:-1]
|
||||
|
||||
if self.dim == 1:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
pass
|
||||
elif self.dim == 2:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr = outStr + printH(self.hy, outStr='\n hy:')
|
||||
elif self.dim == 3:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr = outStr + '\n nCz: {0:d}'.format(self.nCz)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr = outStr + printH(self.hy, outStr='\n hy:')
|
||||
outStr = outStr + printH(self.hz, outStr='\n hz:')
|
||||
|
||||
return outStr
|
||||
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
@property
|
||||
def vol(self):
|
||||
"""Construct cell volumes of the 3D model as 1d array."""
|
||||
if getattr(self, '_vol', None) is None:
|
||||
vh = self.h
|
||||
# Compute cell volumes
|
||||
if self.dim == 1:
|
||||
self._vol = Utils.mkvc(vh[0])
|
||||
elif self.dim == 2:
|
||||
# Cell sizes in each direction
|
||||
self._vol = Utils.mkvc(np.outer(vh[0], vh[1]))
|
||||
elif self.dim == 3:
|
||||
# Cell sizes in each direction
|
||||
self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2]))
|
||||
return self._vol
|
||||
|
||||
@property
|
||||
def area(self):
|
||||
"""Construct face areas of the 3D model as 1d array."""
|
||||
if getattr(self, '_area', None) is None:
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 1):
|
||||
self._area = np.ones(n[0]+1)
|
||||
elif(self.dim == 2):
|
||||
area1 = np.outer(np.ones(n[0]+1), vh[1])
|
||||
area2 = np.outer(vh[0], np.ones(n[1]+1))
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
|
||||
elif(self.dim == 3):
|
||||
area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2])))
|
||||
area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
|
||||
area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
|
||||
return self._area
|
||||
|
||||
@property
|
||||
def edge(self):
|
||||
"""Construct edge legnths of the 3D model as 1d array."""
|
||||
if getattr(self, '_edge', None) is None:
|
||||
# Ensure that we are working with column vectors
|
||||
vh = self.h
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute edge lengths
|
||||
if(self.dim == 1):
|
||||
self._edge = Utils.mkvc(vh[0])
|
||||
elif(self.dim == 2):
|
||||
l1 = np.outer(vh[0], np.ones(n[1]+1))
|
||||
l2 = np.outer(np.ones(n[0]+1), vh[1])
|
||||
self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)]
|
||||
elif(self.dim == 3):
|
||||
l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
|
||||
l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
|
||||
l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
|
||||
self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)]
|
||||
return self._edge
|
||||
|
||||
@property
|
||||
def faceBoundaryInd(self):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
"""
|
||||
if self.dim==1:
|
||||
indxd = (self.gridFx==min(self.gridFx))
|
||||
indxu = (self.gridFx==max(self.gridFx))
|
||||
return indxd, indxu
|
||||
elif self.dim==2:
|
||||
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
|
||||
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
|
||||
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
|
||||
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif self.dim==3:
|
||||
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
|
||||
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
|
||||
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
|
||||
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
|
||||
indzd = (self.gridFz[:,2]==min(self.gridFz[:,2]))
|
||||
indzu = (self.gridFz[:,2]==max(self.gridFz[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
@property
|
||||
def cellBoundaryInd(self):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
"""
|
||||
if self.dim==1:
|
||||
indxd = (self.gridCC==min(self.gridCC))
|
||||
indxu = (self.gridCC==max(self.gridCC))
|
||||
return indxd, indxu
|
||||
elif self.dim==2:
|
||||
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
|
||||
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
|
||||
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
|
||||
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif self.dim==3:
|
||||
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
|
||||
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
|
||||
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
|
||||
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
|
||||
indzd = (self.gridCC[:,2]==min(self.gridCC[:,2]))
|
||||
indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
print('Welcome to tensor mesh!')
|
||||
|
||||
|
||||
@@ -516,6 +516,106 @@ class TensorView(object):
|
||||
return animate(fig, animateFrame, frames=len(frames))
|
||||
|
||||
|
||||
|
||||
class LomView(object):
|
||||
"""
|
||||
Provides viewing functions for LogicallyOrthogonalMesh
|
||||
|
||||
This class is inherited by LogicallyOrthogonalMesh
|
||||
|
||||
"""
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def plotGrid(self, length=0.05, showIt=False):
|
||||
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
|
||||
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh, Utils
|
||||
X, Y = Utils.exampleLomGird([3,3],'rotate')
|
||||
M = Mesh.LogicallyOrthogonalMesh([X, Y])
|
||||
M.plotGrid(showIt=True)
|
||||
|
||||
"""
|
||||
NN = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
fig = plt.figure(2)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2]
|
||||
Y = np.r_[Y1, Y2]
|
||||
|
||||
plt.plot(X, Y)
|
||||
|
||||
plt.hold(True)
|
||||
Nx = self.r(self.normals, 'F', 'Fx', 'V')
|
||||
Ny = self.r(self.normals, 'F', 'Fy', 'V')
|
||||
Tx = self.r(self.tangents, 'E', 'Ex', 'V')
|
||||
Ty = self.r(self.tangents, 'E', 'Ey', 'V')
|
||||
|
||||
plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
|
||||
|
||||
nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
|
||||
plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
|
||||
plt.plot(nX, nY, 'r-')
|
||||
|
||||
nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
|
||||
#plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
|
||||
plt.plot(nX, nY, 'g-')
|
||||
|
||||
tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
|
||||
tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
|
||||
plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
|
||||
plt.plot(tX, tY, 'r-')
|
||||
|
||||
nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
|
||||
#plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
|
||||
plt.plot(nX, nY, 'g-')
|
||||
plt.axis('equal')
|
||||
|
||||
elif self.dim == 3:
|
||||
fig = plt.figure(3)
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
|
||||
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
|
||||
Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
|
||||
|
||||
X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
|
||||
Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
|
||||
Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2, X3]
|
||||
Y = np.r_[Y1, Y2, Y3]
|
||||
Z = np.r_[Z1, Z2, Z3]
|
||||
|
||||
plt.plot(X, Y, 'b', zs=Z)
|
||||
ax.set_zlabel('x3')
|
||||
|
||||
ax.grid(True)
|
||||
ax.hold(False)
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
|
||||
if showIt: plt.show()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import *
|
||||
mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
|
||||
@@ -525,6 +625,7 @@ if __name__ == '__main__':
|
||||
q = Utils.mkvc(q)
|
||||
A = M.faceDiv*M.cellGrad
|
||||
b = Solver(A).solve(q)
|
||||
|
||||
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, pcolorOpts={'alpha':0.8})
|
||||
M2 = Mesh.TensorMesh([10,20],x0=[10,5])
|
||||
f = np.r_[np.sin(M2.gridFx[:,0]*2*np.pi), np.sin(M2.gridFy[:,1]*2*np.pi)]
|
||||
@@ -1,8 +1,5 @@
|
||||
from Cyl1DMesh import Cyl1DMesh
|
||||
from TensorMesh import TensorMesh
|
||||
from TreeMesh import TreeMesh
|
||||
from CylMesh import CylMesh
|
||||
from Cyl1DMesh import Cyl1DMesh
|
||||
from LogicallyRectMesh import LogicallyRectMesh
|
||||
from BaseMesh import BaseMesh, BaseRectangularMesh
|
||||
from TensorView import TensorView
|
||||
from InnerProducts import InnerProducts
|
||||
from DiffOperators import DiffOperators
|
||||
from TreeMesh import TreeMesh
|
||||
|
||||
+1
-1
@@ -337,7 +337,7 @@ if __name__ == '__main__':
|
||||
|
||||
D = M.faceDiv
|
||||
G = M.cellGrad
|
||||
Msig = M.getFaceMass()
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
|
||||
@@ -209,6 +209,7 @@ class BaseRx(object):
|
||||
def nD(self):
|
||||
return self.locs.shape[0]
|
||||
|
||||
|
||||
class BaseTx(object):
|
||||
"""SimPEG Transmitter Object"""
|
||||
|
||||
|
||||
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
|
||||
from numpy.linalg import norm
|
||||
from SimPEG.Utils import mkvc, sdiag
|
||||
from SimPEG import Utils
|
||||
from SimPEG.Mesh import TensorMesh, LogicallyRectMesh
|
||||
from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
@@ -88,10 +88,7 @@ class OrderTest(unittest.TestCase):
|
||||
"""
|
||||
if 'TensorMesh' in self._meshType:
|
||||
if 'uniform' in self._meshType:
|
||||
h1 = np.ones(nc)/nc
|
||||
h2 = np.ones(nc)/nc
|
||||
h3 = np.ones(nc)/nc
|
||||
h = [h1, h2, h3]
|
||||
h = [nc, nc, nc]
|
||||
elif 'random' in self._meshType:
|
||||
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
@@ -104,6 +101,20 @@ class OrderTest(unittest.TestCase):
|
||||
max_h = max([np.max(hi) for hi in self.M.h])
|
||||
return max_h
|
||||
|
||||
elif 'CylMesh' in self._meshType:
|
||||
if 'uniform' in self._meshType:
|
||||
h = [nc, nc, nc]
|
||||
else:
|
||||
raise Exception('Unexpected meshType')
|
||||
|
||||
if self.meshDimension == 2:
|
||||
self.M = CylMesh([h[0], 1, h[2]])
|
||||
max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]])
|
||||
elif self.meshDimension == 3:
|
||||
self.M = CylMesh(h)
|
||||
max_h = max([np.max(hi) for hi in self.M.h])
|
||||
return max_h
|
||||
|
||||
elif 'LRM' in self._meshType:
|
||||
if 'uniform' in self._meshType:
|
||||
kwrd = 'rect'
|
||||
|
||||
@@ -22,7 +22,7 @@ class TestSolver(unittest.TestCase):
|
||||
|
||||
D = M.faceDiv
|
||||
G = M.cellGrad
|
||||
Msig = M.getFaceMass()
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
import unittest
|
||||
import sys
|
||||
from SimPEG.Mesh import BaseRectangularMesh
|
||||
from SimPEG.Mesh.BaseMesh import BaseRectangularMesh
|
||||
import numpy as np
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,327 @@
|
||||
import unittest
|
||||
import sys
|
||||
from SimPEG import *
|
||||
from TestUtils import OrderTest
|
||||
|
||||
|
||||
class TestCyl2DMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
hx = np.r_[1,1,0.5]
|
||||
hz = np.r_[2,1]
|
||||
self.mesh = Mesh.CylMesh([hx, 1,hz])
|
||||
|
||||
def test_dim(self):
|
||||
self.assertTrue(self.mesh.dim == 3)
|
||||
|
||||
def test_nC(self):
|
||||
self.assertTrue(self.mesh.nC == 6)
|
||||
self.assertTrue(self.mesh.nCx == 3)
|
||||
self.assertTrue(self.mesh.nCy == 1)
|
||||
self.assertTrue(self.mesh.nCz == 2)
|
||||
self.assertTrue(np.all(self.mesh.vnC == [3, 1, 2]))
|
||||
|
||||
def test_nN(self):
|
||||
self.assertTrue(self.mesh.nN == 0)
|
||||
self.assertTrue(self.mesh.nNx == 3)
|
||||
self.assertTrue(self.mesh.nNy == 0)
|
||||
self.assertTrue(self.mesh.nNz == 3)
|
||||
self.assertTrue(np.all(self.mesh.vnN == [3, 0, 3]))
|
||||
|
||||
def test_nF(self):
|
||||
self.assertTrue(self.mesh.nFx == 6)
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [3, 1, 2]))
|
||||
self.assertTrue(self.mesh.nFy == 0)
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [3, 0, 2]))
|
||||
self.assertTrue(self.mesh.nFz == 9)
|
||||
self.assertTrue(np.all(self.mesh.vnFz == [3, 1, 3]))
|
||||
self.assertTrue(self.mesh.nF == 15)
|
||||
self.assertTrue(np.all(self.mesh.vnF == [6, 0, 9]))
|
||||
|
||||
def test_nE(self):
|
||||
self.assertTrue(self.mesh.nEx == 0)
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [3, 0, 3]))
|
||||
self.assertTrue(self.mesh.nEy == 9)
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [3, 1, 3]))
|
||||
self.assertTrue(self.mesh.nEz == 0)
|
||||
self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2]))
|
||||
self.assertTrue(self.mesh.nE == 9)
|
||||
self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0]))
|
||||
|
||||
def test_vectorsCC(self):
|
||||
v = np.r_[0.5, 1.5, 2.25]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0)
|
||||
v = np.r_[0]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0)
|
||||
v = np.r_[1, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0)
|
||||
|
||||
def test_vectorsN(self):
|
||||
v = np.r_[1, 2, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0)
|
||||
v = np.r_[0]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0)
|
||||
v = np.r_[0, 2, 3.]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0)
|
||||
|
||||
def test_edge(self):
|
||||
edge = np.r_[1, 2, 2.5, 1, 2, 2.5, 1, 2, 2.5] * 2 * np.pi
|
||||
self.assertTrue(np.linalg.norm((edge-self.mesh.edge)) == 0)
|
||||
|
||||
def test_area(self):
|
||||
r = np.r_[0, 1, 2, 2.5]
|
||||
a = r[1:]*2*np.pi
|
||||
areaX = np.r_[2*a,a]
|
||||
a = (r[1:]**2 - r[:-1]**2)*np.pi
|
||||
areaZ = np.r_[a,a,a]
|
||||
area = np.r_[areaX, areaZ]
|
||||
self.assertTrue(np.linalg.norm((area-self.mesh.area)) == 0)
|
||||
|
||||
def test_vol(self):
|
||||
r = np.r_[0, 1, 2, 2.5]
|
||||
a = (r[1:]**2 - r[:-1]**2)*np.pi
|
||||
vol = np.r_[2*a,a]
|
||||
self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0)
|
||||
|
||||
def test_gridSizes(self):
|
||||
self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3))
|
||||
self.assertTrue(self.mesh.gridN.shape == (9, 3))
|
||||
|
||||
self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3))
|
||||
self.assertTrue(self.mesh.gridFy is None)
|
||||
self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3))
|
||||
|
||||
self.assertTrue(self.mesh.gridEx is None)
|
||||
self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3))
|
||||
self.assertTrue(self.mesh.gridEz is None)
|
||||
|
||||
def test_gridCC(self):
|
||||
x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25]
|
||||
y = np.zeros(6)
|
||||
z = np.r_[1,1,1,2.5,2.5,2.5]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridCC).ravel()) == 0)
|
||||
|
||||
def test_gridN(self):
|
||||
x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5]
|
||||
y = np.zeros(9)
|
||||
z = np.r_[0,0,0,2,2,2,3,3,3.]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridN).ravel()) == 0)
|
||||
|
||||
def test_gridFx(self):
|
||||
x = np.r_[1,2,2.5,1,2,2.5]
|
||||
y = np.zeros(6)
|
||||
z = np.r_[1,1,1,2.5,2.5,2.5]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridFx).ravel()) == 0)
|
||||
|
||||
def test_gridFz(self):
|
||||
x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25,0.5,1.5,2.25]
|
||||
y = np.zeros(9)
|
||||
z = np.r_[0,0,0,2,2,2,3,3,3.]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridFz).ravel()) == 0)
|
||||
|
||||
def test_gridEy(self):
|
||||
x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5]
|
||||
y = np.zeros(9)
|
||||
z = np.r_[0,0,0,2,2,2,3,3,3.]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridEy).ravel()) == 0)
|
||||
|
||||
def test_lightOperators(self):
|
||||
self.assertTrue(self.mesh.nodalGrad is None)
|
||||
|
||||
|
||||
|
||||
MESHTYPES = ['uniformCylMesh']
|
||||
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 2])
|
||||
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
cyl_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
|
||||
cyl_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
|
||||
cylF2 = lambda M, fx, fy: np.vstack((cyl_row2(M.gridFx, fx, fy), cyl_row2(M.gridFz, fx, fy)))
|
||||
|
||||
|
||||
class TestFaceDiv2D(OrderTest):
|
||||
name = "FaceDiv"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
|
||||
funR = lambda r, z: np.sin(2.*np.pi*r)
|
||||
funZ = lambda r, z: np.sin(2.*np.pi*z)
|
||||
|
||||
sol = lambda r, t, z: (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r + 2*np.pi*np.cos(2*np.pi*z)
|
||||
|
||||
Fc = cylF2(self.M, funR, funZ)
|
||||
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
divF = self.M.faceDiv.dot(F)
|
||||
divF_anal = call3(sol, self.M.gridCC)
|
||||
|
||||
err = np.linalg.norm((divF-divF_anal), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestEdgeCurl2D(OrderTest):
|
||||
name = "EdgeCurl"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
# To Recreate or change the functions:
|
||||
|
||||
# import sympy
|
||||
# r,t,z = sympy.symbols('r,t,z')
|
||||
|
||||
# fR = 0
|
||||
# fZ = 0
|
||||
# fT = sympy.sin(2.*sympy.pi*z)
|
||||
|
||||
# print 1/r*sympy.diff(fZ,t) - sympy.diff(fT,z)
|
||||
# print sympy.diff(fR,z) - sympy.diff(fZ,r)
|
||||
# print 1/r*(sympy.diff(r*fT,r) - sympy.diff(fR,t))
|
||||
|
||||
funT = lambda r, t, z: np.sin(2.*np.pi*z)
|
||||
|
||||
solR = lambda r, z: -2.0*np.pi*np.cos(2.0*np.pi*z)
|
||||
solZ = lambda r, z: np.sin(2.0*np.pi*z)/r
|
||||
|
||||
E = call3(funT, self.M.gridEy)
|
||||
|
||||
curlE = self.M.edgeCurl.dot(E)
|
||||
|
||||
Fc = cylF2(self.M, solR, solZ)
|
||||
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
|
||||
curlE_anal = self.M.projectFaceVector(Fc)
|
||||
|
||||
err = np.linalg.norm((curlE-curlE_anal), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
# class TestInnerProducts2D(OrderTest):
|
||||
# """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
# meshTypes = MESHTYPES
|
||||
# meshDimension = 2
|
||||
# meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
# def getError(self):
|
||||
|
||||
# funR = lambda r, t, z: np.cos(2.0*np.pi*z)
|
||||
# funT = lambda r, t, z: 0*t
|
||||
# funZ = lambda r, t, z: np.sin(2.0*np.pi*r)
|
||||
|
||||
# call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
|
||||
# sigma1 = lambda r, t, z: z+1
|
||||
# sigma2 = lambda r, t, z: r*z+50
|
||||
# sigma3 = lambda r, t, z: 3+t*r
|
||||
# sigma4 = lambda r, t, z: 0.1*r*t*z
|
||||
# sigma5 = lambda r, t, z: 0.2*z*r*t
|
||||
# sigma6 = lambda r, t, z: 0.1*t
|
||||
|
||||
# Gc = self.M.gridCC
|
||||
# if self.sigmaTest == 1:
|
||||
# sigma = np.c_[call(sigma1, Gc)]
|
||||
# analytic = 144877./360 # Found using sympy. z=5
|
||||
# elif self.sigmaTest == 2:
|
||||
# sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
||||
# analytic = 189959./120 # Found using sympy. z=5
|
||||
# elif self.sigmaTest == 3:
|
||||
# sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
# analytic = 781427./360 # Found using sympy. z=5
|
||||
|
||||
# if self.location == 'edges':
|
||||
# E = call(funT, self.M.gridEy)
|
||||
# A = self.M.getEdgeInnerProduct(sigma)
|
||||
# numeric = E.T.dot(A.dot(E))
|
||||
# elif self.location == 'faces':
|
||||
# Fr = call(funR, self.M.gridFx)
|
||||
# Fz = call(funZ, self.M.gridFz)
|
||||
# A = self.M.getFaceInnerProduct(sigma)
|
||||
# F = np.r_[Fr,Fz]
|
||||
# numeric = F.T.dot(A.dot(F))
|
||||
|
||||
# print numeric
|
||||
# err = np.abs(numeric - analytic)
|
||||
# return err
|
||||
|
||||
# def test_order1_faces(self):
|
||||
# self.name = "2D Face Inner Product - Isotropic"
|
||||
# self.location = 'faces'
|
||||
# self.sigmaTest = 1
|
||||
# self.orderTest()
|
||||
|
||||
|
||||
class TestCyl3DMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
hx = np.r_[1,1,0.5]
|
||||
hy = np.r_[np.pi, np.pi]
|
||||
hz = np.r_[2,1]
|
||||
self.mesh = Mesh.CylMesh([hx, hy,hz])
|
||||
|
||||
def test_dim(self):
|
||||
self.assertTrue(self.mesh.dim == 3)
|
||||
|
||||
def test_nC(self):
|
||||
self.assertTrue(self.mesh.nCx == 3)
|
||||
self.assertTrue(self.mesh.nCy == 2)
|
||||
self.assertTrue(self.mesh.nCz == 2)
|
||||
self.assertTrue(np.all(self.mesh.vnC == [3, 2, 2]))
|
||||
|
||||
def test_nN(self):
|
||||
self.assertTrue(self.mesh.nN == 24)
|
||||
self.assertTrue(self.mesh.nNx == 4)
|
||||
self.assertTrue(self.mesh.nNy == 2)
|
||||
self.assertTrue(self.mesh.nNz == 3)
|
||||
self.assertTrue(np.all(self.mesh.vnN == [4, 2, 3]))
|
||||
|
||||
def test_nF(self):
|
||||
self.assertTrue(self.mesh.nFx == 12)
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [3, 2, 2]))
|
||||
self.assertTrue(self.mesh.nFy == 12)
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [3, 2, 2]))
|
||||
self.assertTrue(self.mesh.nFz == 18)
|
||||
self.assertTrue(np.all(self.mesh.vnFz == [3, 2, 3]))
|
||||
self.assertTrue(self.mesh.nF == 42)
|
||||
self.assertTrue(np.all(self.mesh.vnF == [12, 12, 18]))
|
||||
|
||||
def test_nE(self):
|
||||
self.assertTrue(self.mesh.nEx == 18)
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [3, 2, 3]))
|
||||
self.assertTrue(self.mesh.nEy == 18)
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [3, 2, 3]))
|
||||
self.assertTrue(self.mesh.nEz == 12 + 2)
|
||||
self.assertTrue(self.mesh.vnEz is None)
|
||||
self.assertTrue(self.mesh.nE == 50)
|
||||
self.assertTrue(np.all(self.mesh.vnE == [18, 18, 14]))
|
||||
|
||||
def test_vectorsCC(self):
|
||||
v = np.r_[0.5, 1.5, 2.25]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0)
|
||||
v = np.r_[0, np.pi]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0)
|
||||
v = np.r_[1, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0)
|
||||
|
||||
def test_vectorsN(self):
|
||||
v = np.r_[0, 1, 2, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0)
|
||||
v = np.r_[np.pi/2, 1.5*np.pi]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0)
|
||||
v = np.r_[0, 2, 3]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -3,6 +3,7 @@ import unittest
|
||||
from TestUtils import OrderTest
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
#TODO: 'randomTensorMesh'
|
||||
MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
|
||||
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
|
||||
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
|
||||
Reference in New Issue
Block a user