mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-04 04:45:06 +08:00
Initial clean up of the CylMesh code. Inherits from TensorMesh. Many more things to do!
This commit is contained in:
rowanc1
@@ -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
|
||||
|
||||
@@ -1,476 +0,0 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from scipy.constants import pi
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag
|
||||
|
||||
class Cyl1DMesh(object):
|
||||
"""
|
||||
Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems
|
||||
"""
|
||||
|
||||
_meshType = 'CYL1D'
|
||||
|
||||
def __init__(self, h, z0=None):
|
||||
assert len(h) == 2, "len(h) must equal 2"
|
||||
if z0 is not None:
|
||||
assert z0.size == 1, "z0.size must equal 1"
|
||||
|
||||
for i, h_i in enumerate(h):
|
||||
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)
|
||||
|
||||
# Ensure h contains 1D vectors
|
||||
self._h = [mkvc(x.astype(float)) for x in h]
|
||||
|
||||
if z0 is None:
|
||||
z0 = 0
|
||||
self._z0 = z0
|
||||
|
||||
####################################################
|
||||
# Mesh properties
|
||||
####################################################
|
||||
|
||||
def h():
|
||||
doc = "list containing the width of each cell"
|
||||
def fget(self):
|
||||
return self._h
|
||||
return locals()
|
||||
h = property(**h())
|
||||
|
||||
@property
|
||||
def dim(self): return 2
|
||||
|
||||
def z0():
|
||||
doc = "The z-origin"
|
||||
def fget(self):
|
||||
return self._z0
|
||||
return locals()
|
||||
z0 = property(**z0())
|
||||
|
||||
def hr():
|
||||
doc = "Width of the cells in the r direction"
|
||||
def fget(self):
|
||||
return self._h[0]
|
||||
return locals()
|
||||
hr = property(**hr())
|
||||
|
||||
def hz():
|
||||
doc = "Width of the cells in the z direction"
|
||||
def fget(self):
|
||||
return self._h[1]
|
||||
return locals()
|
||||
hz = property(**hz())
|
||||
|
||||
####################################################
|
||||
# Counting
|
||||
####################################################
|
||||
|
||||
def nCx():
|
||||
doc = "Number of cells in the radial direction"
|
||||
fget = lambda self: self.hr.size
|
||||
return locals()
|
||||
nCx = property(**nCx())
|
||||
|
||||
def nCz():
|
||||
doc = "Number of cells in the z direction"
|
||||
fget = lambda self: self.hz.size
|
||||
return locals()
|
||||
nCz = property(**nCz())
|
||||
|
||||
def nC():
|
||||
doc = "Total number of cells"
|
||||
fget = lambda self: self.nCx * self.nCz
|
||||
return locals()
|
||||
nC = property(**nC())
|
||||
|
||||
def vnC():
|
||||
doc = "Total number of cells in each direction"
|
||||
fget = lambda self: np.array([self.nCx, self.nCz])
|
||||
return locals()
|
||||
vnC = property(**vnC())
|
||||
|
||||
def nNr():
|
||||
doc = "Number of nodes in the radial direction"
|
||||
fget = lambda self: self.hr.size
|
||||
return locals()
|
||||
nNr = property(**nNr())
|
||||
|
||||
def nNz():
|
||||
doc = "Number of nodes in the radial direction"
|
||||
fget = lambda self: self.hz.size + 1
|
||||
return locals()
|
||||
nNz = property(**nNz())
|
||||
|
||||
def nN():
|
||||
doc = "Total number of nodes"
|
||||
fget = lambda self: self.nNr * self.nNz
|
||||
return locals()
|
||||
nN = property(**nN())
|
||||
|
||||
def nFr():
|
||||
doc = "Number of r faces"
|
||||
fget = lambda self: self.nNr * self.nCz
|
||||
return locals()
|
||||
nFr = property(**nFr())
|
||||
|
||||
def vnFz():
|
||||
doc = "Number of z faces"
|
||||
fget = lambda self: self.nNz * self.nCx
|
||||
return locals()
|
||||
vnFz = property(**vnFz())
|
||||
|
||||
def vnF():
|
||||
doc = "Total number of faces in each direction"
|
||||
fget = lambda self: np.array([self.nFr, self.vnFz])
|
||||
return locals()
|
||||
vnF = property(**vnF())
|
||||
|
||||
def nF():
|
||||
doc = "Total number of faces"
|
||||
fget = lambda self: self.nFr + self.vnFz
|
||||
return locals()
|
||||
nF = property(**nF())
|
||||
|
||||
def nE():
|
||||
doc = "Number of edges"
|
||||
fget = lambda self: self.nN
|
||||
return locals()
|
||||
nE = property(**nE())
|
||||
|
||||
####################################################
|
||||
# Vectors & Grids
|
||||
####################################################
|
||||
|
||||
def vectorNr():
|
||||
doc = "Nodal grid vector (1D) in the r direction"
|
||||
fget = lambda self: self.hr.cumsum()
|
||||
return locals()
|
||||
vectorNr = property(**vectorNr())
|
||||
|
||||
def vectorNz():
|
||||
doc = "Nodal grid vector (1D) in the z direction"
|
||||
fget = lambda self: np.r_[0, self.hz.cumsum()] + self._z0
|
||||
return locals()
|
||||
vectorNz = property(**vectorNz())
|
||||
|
||||
def vectorCCr():
|
||||
doc = "Cell centered grid vector (1D) in the r direction"
|
||||
fget = lambda self: np.r_[0, self.hr.cumsum()[1:] - self.hr[1:]/2]
|
||||
return locals()
|
||||
vectorCCr = property(**vectorCCr())
|
||||
|
||||
def vectorCCz():
|
||||
doc = "Cell centered grid vector (1D) in the z direction"
|
||||
fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0
|
||||
return locals()
|
||||
vectorCCz = property(**vectorCCz())
|
||||
|
||||
def gridCC():
|
||||
doc = "Cell-centered grid"
|
||||
def fget(self):
|
||||
if self._gridCC is None:
|
||||
self._gridCC = ndgrid([self.vectorCCr, self.vectorCCz])
|
||||
return self._gridCC
|
||||
return locals()
|
||||
_gridCC = None
|
||||
gridCC = property(**gridCC())
|
||||
|
||||
def gridN():
|
||||
doc = "Nodal grid"
|
||||
def fget(self):
|
||||
if self._gridN is None:
|
||||
self._gridN = ndgrid([self.vectorNr, self.vectorNz])
|
||||
return self._gridN
|
||||
return locals()
|
||||
_gridN = None
|
||||
gridN = property(**gridN())
|
||||
|
||||
def gridFr():
|
||||
doc = "r face grid"
|
||||
def fget(self):
|
||||
if self._gridFr is None:
|
||||
self._gridFr = ndgrid([self.vectorNr, self.vectorCCz])
|
||||
return self._gridFr
|
||||
return locals()
|
||||
_gridFr = None
|
||||
gridFr = property(**gridFr())
|
||||
|
||||
def gridFz():
|
||||
doc = "z face grid"
|
||||
def fget(self):
|
||||
if self._gridFz is None:
|
||||
self._gridFz = ndgrid([self.vectorCCr, self.vectorNz])
|
||||
return self._gridFz
|
||||
return locals()
|
||||
_gridFz = None
|
||||
gridFz = property(**gridFz())
|
||||
|
||||
####################################################
|
||||
# Geometries
|
||||
####################################################
|
||||
|
||||
def edge():
|
||||
doc = "Edge lengths"
|
||||
def fget(self):
|
||||
if self._edge is None:
|
||||
self._edge = 2*pi*self.gridN[:,0]
|
||||
return self._edge
|
||||
return locals()
|
||||
_edge = None
|
||||
edge = property(**edge())
|
||||
|
||||
def area():
|
||||
doc = "Face areas"
|
||||
def fget(self):
|
||||
if self._area is None:
|
||||
areaR = np.kron(self.hz, 2*pi*self.vectorNr)
|
||||
areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2))
|
||||
self._area = np.r_[areaR, areaZ]
|
||||
return self._area
|
||||
return locals()
|
||||
_area = None
|
||||
area = property(**area())
|
||||
|
||||
def vol():
|
||||
doc = "Volume of each cell"
|
||||
def fget(self):
|
||||
if self._vol is None:
|
||||
az = pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2)
|
||||
self._vol = np.kron(self.hz,az)
|
||||
return self._vol
|
||||
return locals()
|
||||
_vol = None
|
||||
vol = property(**vol())
|
||||
|
||||
####################################################
|
||||
# Operators
|
||||
####################################################
|
||||
|
||||
def edgeCurl():
|
||||
doc = "The edgeCurl property."
|
||||
def fget(self):
|
||||
if self._edgeCurl 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
|
||||
return locals()
|
||||
_edgeCurl = None
|
||||
edgeCurl = property(**edgeCurl())
|
||||
|
||||
def aveE2CC():
|
||||
doc = "Averaging operator from cell edges to cell centres"
|
||||
def fget(self):
|
||||
if self._aveE2CC is None:
|
||||
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 = sp.kron(az, ar).T
|
||||
return self._aveE2CC
|
||||
return locals()
|
||||
_aveE2CC = None
|
||||
aveE2CC = property(**aveE2CC())
|
||||
|
||||
def aveF2CC():
|
||||
doc = "Averaging operator from cell faces to cell centres"
|
||||
def fget(self):
|
||||
if self._aveF2CC is None:
|
||||
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
|
||||
Afr = sp.kron(sp.eye(self.nCz),ar)
|
||||
Afz = sp.kron(az,sp.eye(self.nCx))
|
||||
self._aveF2CC = sp.vstack((Afr,Afz)).T
|
||||
return self._aveF2CC
|
||||
return locals()
|
||||
_aveF2CC = None
|
||||
aveF2CC = property(**aveF2CC())
|
||||
|
||||
def getFaceMassDeriv(self):
|
||||
Av = self.aveF2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
def getEdgeMassDeriv(self):
|
||||
Av = self.aveE2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
|
||||
####################################################
|
||||
# Methods
|
||||
####################################################
|
||||
|
||||
|
||||
def getMass(self, materialProp=None, loc='e'):
|
||||
""" Produces mass matricies.
|
||||
|
||||
:param None,float,numpy.ndarray materialProp: property to be averaged (see below)
|
||||
:param str loc: Average to location: 'e'-edges, 'f'-faces
|
||||
: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)
|
||||
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 getInterpolationMat(self, loc, locType='fz'):
|
||||
""" Produces intrpolation matrix
|
||||
|
||||
:param numpy.ndarray loc: Location of points to interpolate to
|
||||
:param str locType: What to interpolate (see below)
|
||||
:rtype: scipy.sparse.csr.csr_matrix
|
||||
:return: M, the intrpolation matrix
|
||||
|
||||
locType can be::
|
||||
|
||||
'fz' -> z-component of field defined on faces
|
||||
'fr' -> r-component of field defined on faces
|
||||
'et' -> theta-component of field defined on edges
|
||||
"""
|
||||
|
||||
loc = np.atleast_2d(loc)
|
||||
|
||||
assert np.all(loc[:,0]<=self.vectorNr.max()) & \
|
||||
np.all(loc[:,1]>=self.vectorNz.min()) & \
|
||||
np.all(loc[:,1]<=self.vectorNz.max()), \
|
||||
"Points outside of mesh"
|
||||
|
||||
|
||||
if locType=='fz':
|
||||
Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float)
|
||||
|
||||
for i, iloc in enumerate(loc):
|
||||
# Point is on a z-interface
|
||||
if np.any(np.abs(self.vectorNz-iloc[1])<0.001):
|
||||
dFz = self.gridFz-iloc #Distance to z faces
|
||||
dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left...
|
||||
indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one
|
||||
if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation)
|
||||
zFL = self.gridFz[indL,:]
|
||||
zFLL = self.gridFz[indL-1,:]
|
||||
Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0])
|
||||
Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0])
|
||||
else:
|
||||
zFL = self.gridFz[indL,:]
|
||||
zFR = self.gridFz[indL+1,:]
|
||||
Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0])
|
||||
Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0])
|
||||
# Point is in a cell
|
||||
else:
|
||||
dFz = self.gridFz-iloc
|
||||
dFz[dFz>0] = np.inf
|
||||
dFz = np.sum(dFz**2, axis=1)
|
||||
|
||||
indBL = np.argmin(dFz) # Face below and to the left
|
||||
indAL = indBL + self.nCx # Face above and to the left
|
||||
|
||||
zF_BL = self.gridFz[indBL,:]
|
||||
zF_AL = self.gridFz[indAL,:]
|
||||
|
||||
dzB = iloc[1] - zF_BL[1] # z-distance to face below
|
||||
dzA = zF_AL[1] - iloc[1] # z-distance to face above
|
||||
|
||||
if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation)
|
||||
zF_BLL = self.gridFz[indBL-1,:]
|
||||
zF_ALL = self.gridFz[indAL-1,:]
|
||||
|
||||
DZ = zF_AL[1] - zF_BL[1]
|
||||
DR = zF_AL[0] - zF_ALL[0]
|
||||
|
||||
drL = iloc[0] - zF_AL[0]
|
||||
drLL = iloc[0] - zF_ALL[0]
|
||||
|
||||
Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR)
|
||||
Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR)
|
||||
Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR)
|
||||
Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR)
|
||||
else:
|
||||
indBR = indBL+1 # Face below and to the right
|
||||
indAR = indAL + 1 # Face above and to the right
|
||||
zF_BR = self.gridFz[indBR,:]
|
||||
|
||||
drL = iloc[0] - zF_BL[0] # r-distance to face on left
|
||||
drR = zF_BR[0] - iloc[0] # r-distance to face on right
|
||||
|
||||
drz = (drL + drR)*(dzB + dzA)
|
||||
Q[i,indBL+self.nFr] = drR*dzA/drz
|
||||
Q[i,indBR+self.nFr] = drL*dzA/drz
|
||||
Q[i,indAL+self.nFr] = drR*dzB/drz
|
||||
Q[i,indAR+self.nFr] = drL*dzB/drz
|
||||
|
||||
elif locType=='fr':
|
||||
raise NotImplementedError('locType==fr')
|
||||
elif locType=='et':
|
||||
raise NotImplementedError('locType==et')
|
||||
else:
|
||||
raise ValueError('Invalid locType')
|
||||
return Q.tocsr()
|
||||
|
||||
def getNearest(self, loc, locType):
|
||||
""" Returns the index of the closest face or edge to a given location
|
||||
|
||||
:param numpy.ndarray loc: Test point
|
||||
:param str locType: Type of location desired (see below)
|
||||
:rtype: int
|
||||
:return: ind:
|
||||
|
||||
locType can be::
|
||||
|
||||
'fz' -> location of nearest z-face
|
||||
'fr' -> location of nearest r-face
|
||||
'et' -> location of nearest edge
|
||||
"""
|
||||
|
||||
if locType=='et':
|
||||
dr = self.gridN[:,0] - loc[0]
|
||||
dz = self.gridN[:,1] - loc[1]
|
||||
elif locType=='fz':
|
||||
dr = self.gridFz[:,0] - loc[0]
|
||||
dz = self.gridFz[:,1] - loc[1]
|
||||
elif locType=='fr':
|
||||
dr = self.gridFr[:,0] - loc[0]
|
||||
dz = self.gridFr[:,1] - loc[1]
|
||||
else:
|
||||
raise ValueError('Invalid locType')
|
||||
R = np.sqrt(dr**2 + dz**2)
|
||||
ind = np.argmin(R)
|
||||
return ind
|
||||
@@ -0,0 +1,337 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from scipy.constants import pi
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag
|
||||
from TensorMesh import TensorMesh
|
||||
|
||||
class CylMesh(TensorMesh):
|
||||
"""
|
||||
CylMesh is a mesh class for cylindrically problems
|
||||
"""
|
||||
|
||||
_meshType = 'CYL'
|
||||
|
||||
def __init__(self, h, x0=None):
|
||||
assert len(h) == 3, "len(h) must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]"
|
||||
|
||||
if x0 is not None:
|
||||
assert x0.size == 3, "x0.size must equal 1"
|
||||
else:
|
||||
x0 = np.r_[0, 0, 0]
|
||||
|
||||
for i, h_i in enumerate(h):
|
||||
if type(h_i) in [int, long, float]:
|
||||
# This gives you something over the unit cylinder.
|
||||
h_i = (2*np.pi if i==1 else 1.)*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.
|
||||
|
||||
assert h[1].sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
|
||||
|
||||
TensorMesh.__init__(self, h, x0)
|
||||
|
||||
|
||||
@property
|
||||
def nNy(self):
|
||||
"""
|
||||
Number of nodes in the y-direction
|
||||
|
||||
:rtype: int
|
||||
:return: nNy
|
||||
"""
|
||||
return self.nCy
|
||||
|
||||
@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 (np.r_[self.nNx, self.nNy, self.nNz]).prod()
|
||||
|
||||
@property
|
||||
def nFy(self):
|
||||
"""
|
||||
Number of y-faces
|
||||
|
||||
:rtype: int
|
||||
:return: nFy
|
||||
"""
|
||||
return (self.vnC + np.r_[0,-1,0][:self.dim]).prod()
|
||||
|
||||
@property
|
||||
def vnF(self):
|
||||
"""Total number of faces in each direction"""
|
||||
return np.array([self.nFr, self.vnFz])
|
||||
|
||||
@property
|
||||
def nF(self):
|
||||
"""Total number of faces"""
|
||||
return self.nFr + self.vnFz
|
||||
|
||||
@property
|
||||
def nE(self):
|
||||
"""Number of edges"""
|
||||
return self.nN
|
||||
|
||||
@property
|
||||
def vectorNx(self):
|
||||
"""Nodal grid vector (1D) in the r direction"""
|
||||
return self.hr.cumsum()
|
||||
|
||||
@property
|
||||
def edge(self):
|
||||
"""Edge lengths"""
|
||||
if getattr(self, '_edge', None) is None:
|
||||
self._edge = 2*pi*self.gridN[:,0]
|
||||
return self._edge
|
||||
|
||||
@property
|
||||
def area(self):
|
||||
"""Face areas"""
|
||||
if getattr(self, '_area', None) is None:
|
||||
areaR = np.kron(self.hz, 2*pi*self.vectorNr)
|
||||
areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-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:
|
||||
az = pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2)
|
||||
self._vol = np.kron(self.hz,az)
|
||||
return self._vol
|
||||
|
||||
####################################################
|
||||
# Operators
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def edgeCurl(self):
|
||||
"""The edgeCurl property."""
|
||||
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 centres"""
|
||||
if getattr(self, '_aveE2CC', None) is None:
|
||||
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 = sp.kron(az, ar).T
|
||||
return self._aveE2CC
|
||||
|
||||
@property
|
||||
def aveF2CC(self):
|
||||
"""Averaging operator from cell faces to cell centres"""
|
||||
if getattr(self, '_aveF2CC', None) is None:
|
||||
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
|
||||
Afr = sp.kron(sp.eye(self.nCz),ar)
|
||||
Afz = sp.kron(az,sp.eye(self.nCx))
|
||||
self._aveF2CC = sp.vstack((Afr,Afz)).T
|
||||
return self._aveF2CC
|
||||
|
||||
def getFaceMassDeriv(self):
|
||||
Av = self.aveF2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
def getEdgeMassDeriv(self):
|
||||
Av = self.aveE2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
|
||||
####################################################
|
||||
# Methods
|
||||
####################################################
|
||||
|
||||
|
||||
def getMass(self, materialProp=None, loc='e'):
|
||||
""" Produces mass matricies.
|
||||
|
||||
:param None,float,numpy.ndarray materialProp: property to be averaged (see below)
|
||||
:param str loc: Average to location: 'e'-edges, 'f'-faces
|
||||
: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)
|
||||
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 getInterpolationMat(self, loc, locType='fz'):
|
||||
""" Produces intrpolation matrix
|
||||
|
||||
:param numpy.ndarray loc: Location of points to interpolate to
|
||||
:param str locType: What to interpolate (see below)
|
||||
:rtype: scipy.sparse.csr.csr_matrix
|
||||
:return: M, the intrpolation matrix
|
||||
|
||||
locType can be::
|
||||
|
||||
'fz' -> z-component of field defined on faces
|
||||
'fr' -> r-component of field defined on faces
|
||||
'et' -> theta-component of field defined on edges
|
||||
"""
|
||||
|
||||
loc = np.atleast_2d(loc)
|
||||
|
||||
assert np.all(loc[:,0]<=self.vectorNr.max()) & \
|
||||
np.all(loc[:,1]>=self.vectorNz.min()) & \
|
||||
np.all(loc[:,1]<=self.vectorNz.max()), \
|
||||
"Points outside of mesh"
|
||||
|
||||
|
||||
if locType=='fz':
|
||||
Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float)
|
||||
|
||||
for i, iloc in enumerate(loc):
|
||||
# Point is on a z-interface
|
||||
if np.any(np.abs(self.vectorNz-iloc[1])<0.001):
|
||||
dFz = self.gridFz-iloc #Distance to z faces
|
||||
dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left...
|
||||
indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one
|
||||
if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation)
|
||||
zFL = self.gridFz[indL,:]
|
||||
zFLL = self.gridFz[indL-1,:]
|
||||
Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0])
|
||||
Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0])
|
||||
else:
|
||||
zFL = self.gridFz[indL,:]
|
||||
zFR = self.gridFz[indL+1,:]
|
||||
Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0])
|
||||
Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0])
|
||||
# Point is in a cell
|
||||
else:
|
||||
dFz = self.gridFz-iloc
|
||||
dFz[dFz>0] = np.inf
|
||||
dFz = np.sum(dFz**2, axis=1)
|
||||
|
||||
indBL = np.argmin(dFz) # Face below and to the left
|
||||
indAL = indBL + self.nCx # Face above and to the left
|
||||
|
||||
zF_BL = self.gridFz[indBL,:]
|
||||
zF_AL = self.gridFz[indAL,:]
|
||||
|
||||
dzB = iloc[1] - zF_BL[1] # z-distance to face below
|
||||
dzA = zF_AL[1] - iloc[1] # z-distance to face above
|
||||
|
||||
if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation)
|
||||
zF_BLL = self.gridFz[indBL-1,:]
|
||||
zF_ALL = self.gridFz[indAL-1,:]
|
||||
|
||||
DZ = zF_AL[1] - zF_BL[1]
|
||||
DR = zF_AL[0] - zF_ALL[0]
|
||||
|
||||
drL = iloc[0] - zF_AL[0]
|
||||
drLL = iloc[0] - zF_ALL[0]
|
||||
|
||||
Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR)
|
||||
Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR)
|
||||
Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR)
|
||||
Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR)
|
||||
else:
|
||||
indBR = indBL+1 # Face below and to the right
|
||||
indAR = indAL + 1 # Face above and to the right
|
||||
zF_BR = self.gridFz[indBR,:]
|
||||
|
||||
drL = iloc[0] - zF_BL[0] # r-distance to face on left
|
||||
drR = zF_BR[0] - iloc[0] # r-distance to face on right
|
||||
|
||||
drz = (drL + drR)*(dzB + dzA)
|
||||
Q[i,indBL+self.nFr] = drR*dzA/drz
|
||||
Q[i,indBR+self.nFr] = drL*dzA/drz
|
||||
Q[i,indAL+self.nFr] = drR*dzB/drz
|
||||
Q[i,indAR+self.nFr] = drL*dzB/drz
|
||||
|
||||
elif locType=='fr':
|
||||
raise NotImplementedError('locType==fr')
|
||||
elif locType=='et':
|
||||
raise NotImplementedError('locType==et')
|
||||
else:
|
||||
raise ValueError('Invalid locType')
|
||||
return Q.tocsr()
|
||||
|
||||
def getNearest(self, loc, locType):
|
||||
""" Returns the index of the closest face or edge to a given location
|
||||
|
||||
:param numpy.ndarray loc: Test point
|
||||
:param str locType: Type of location desired (see below)
|
||||
:rtype: int
|
||||
:return: ind:
|
||||
|
||||
locType can be::
|
||||
|
||||
'fz' -> location of nearest z-face
|
||||
'fr' -> location of nearest r-face
|
||||
'et' -> location of nearest edge
|
||||
"""
|
||||
|
||||
if locType=='et':
|
||||
dr = self.gridN[:,0] - loc[0]
|
||||
dz = self.gridN[:,1] - loc[1]
|
||||
elif locType=='fz':
|
||||
dr = self.gridFz[:,0] - loc[0]
|
||||
dz = self.gridFz[:,1] - loc[1]
|
||||
elif locType=='fr':
|
||||
dr = self.gridFr[:,0] - loc[0]
|
||||
dz = self.gridFr[:,1] - loc[1]
|
||||
else:
|
||||
raise ValueError('Invalid locType')
|
||||
R = np.sqrt(dr**2 + dz**2)
|
||||
ind = np.argmin(R)
|
||||
return ind
|
||||
+145
-199
@@ -105,226 +105,172 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
|
||||
return outStr
|
||||
|
||||
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 h(self):
|
||||
"""h is a list containing the cell widths of the tensor mesh in each dimension."""
|
||||
return self._h
|
||||
|
||||
def hx():
|
||||
doc = "Width of cells in the x direction"
|
||||
fget = lambda self: self._h[0]
|
||||
return locals()
|
||||
hx = property(**hx())
|
||||
@property
|
||||
def hx(self):
|
||||
"Width of cells in the x direction"
|
||||
return self._h[0]
|
||||
|
||||
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 hy(self):
|
||||
"Width of cells in the y direction"
|
||||
return None if self.dim < 2 else self._h[1]
|
||||
|
||||
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 hz(self):
|
||||
"Width of cells in the z direction"
|
||||
return None if self.dim < 3 else self._h[2]
|
||||
|
||||
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 vectorNx(self):
|
||||
"""Nodal grid vector (1D) in the x direction."""
|
||||
return np.r_[0., self.hx.cumsum()] + self.x0[0]
|
||||
|
||||
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 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 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 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 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 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 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 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 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 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 gridCC():
|
||||
doc = "Cell-centered grid."
|
||||
@property
|
||||
def gridCC(self):
|
||||
"""Cell-centered grid."""
|
||||
if getattr(self, '_gridCC', None) is None:
|
||||
self._gridCC = Utils.ndgrid(self.getTensor('CC'))
|
||||
return self._gridCC
|
||||
|
||||
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 gridN(self):
|
||||
"""Nodal grid."""
|
||||
if getattr(self, '_gridN', None) is None:
|
||||
self._gridN = Utils.ndgrid(self.getTensor('N'))
|
||||
return self._gridN
|
||||
|
||||
def gridN():
|
||||
doc = "Nodal grid."
|
||||
@property
|
||||
def gridFx(self):
|
||||
"""Face staggered grid in the x direction."""
|
||||
if getattr(self, '_gridFx', None) is None:
|
||||
self._gridFx = Utils.ndgrid(self.getTensor('Fx'))
|
||||
return self._gridFx
|
||||
|
||||
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())
|
||||
@property
|
||||
def gridFy(self):
|
||||
"""Face staggered grid in the y direction."""
|
||||
if getattr(self, '_gridFy', None) is None and self.dim > 1:
|
||||
self._gridFy = Utils.ndgrid(self.getTensor('Fy'))
|
||||
return self._gridFy
|
||||
|
||||
def gridFx():
|
||||
doc = "Face staggered grid in the x direction."
|
||||
@property
|
||||
def gridFz(self):
|
||||
"""Face staggered grid in the z direction."""
|
||||
if getattr(self, '_gridFz', None) is None and self.dim > 2:
|
||||
self._gridFz = Utils.ndgrid(self.getTensor('Fz'))
|
||||
return self._gridFz
|
||||
|
||||
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())
|
||||
@property
|
||||
def gridEx(self):
|
||||
"""Edge staggered grid in the x direction."""
|
||||
if getattr(self, '_gridEx', None) is None:
|
||||
self._gridEx = Utils.ndgrid(self.getTensor('Ex'))
|
||||
return self._gridEx
|
||||
|
||||
def gridFy():
|
||||
doc = "Face staggered grid in the y direction."
|
||||
@property
|
||||
def gridEy(self):
|
||||
"""Edge staggered grid in the y direction."""
|
||||
if getattr(self, '_gridEy', None) is None and self.dim > 1:
|
||||
self._gridEy = Utils.ndgrid(self.getTensor('Ey'))
|
||||
return self._gridEy
|
||||
|
||||
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())
|
||||
@property
|
||||
def gridEz(self):
|
||||
"""Edge staggered grid in the z direction."""
|
||||
if getattr(self, '_gridEz', None) is None and self.dim > 2:
|
||||
self._gridEz = Utils.ndgrid(self.getTensor('Ez'))
|
||||
return self._gridEz
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
def vol():
|
||||
doc = "Construct cell volumes of the 3D model as 1d array."
|
||||
@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
|
||||
|
||||
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())
|
||||
@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
|
||||
|
||||
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())
|
||||
@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
|
||||
|
||||
# --------------- Methods ---------------------
|
||||
|
||||
|
||||
@@ -1,8 +1,8 @@
|
||||
from Cyl1DMesh import Cyl1DMesh
|
||||
from TensorMesh import TensorMesh
|
||||
from TreeMesh import TreeMesh
|
||||
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
|
||||
from BaseMesh import BaseMesh, BaseRectangularMesh
|
||||
from TensorMesh import TensorMesh
|
||||
from CylMesh import CylMesh
|
||||
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
|
||||
from TreeMesh import TreeMesh
|
||||
from TensorView import TensorView
|
||||
from LomView import LomView
|
||||
from InnerProducts import InnerProducts
|
||||
|
||||
@@ -0,0 +1,29 @@
|
||||
import unittest
|
||||
import sys
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class TestCylMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
hx = np.ones(10)
|
||||
hz = np.ones(5)
|
||||
|
||||
self.mesh = Mesh.CylMesh([hx, 1,hz])
|
||||
|
||||
def test_cylMeshInheritance(self):
|
||||
self.assertTrue(isinstance(self.mesh, Mesh.BaseMesh))
|
||||
|
||||
def test_cylMeshDimensions(self):
|
||||
self.assertTrue(self.mesh.dim == 3)
|
||||
|
||||
def test_cylMesh_nc(self):
|
||||
self.assertTrue(np.all(self.mesh.vnC == [10, 1, 5]))
|
||||
|
||||
def test_cylMesh_nc_xyz(self):
|
||||
self.assertTrue(self.mesh.nCx == 10)
|
||||
self.assertTrue(self.mesh.nCy == 1)
|
||||
self.assertTrue(self.mesh.nCz == 5)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user