wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 19:48:52 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
f2e13182bf75deea943feb672d983cf4c7d46186
simpeg/SimPEG/Mesh/CylMesh.py
T

402 lines
15 KiB
Python
Raw Blame History

import numpy as np
import scipy.sparse as sp
from scipy.constants import pi
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
from TensorMesh import BaseTensorMesh, BaseRectangularMesh
from InnerProducts import InnerProducts
from View import CylView
class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
CylMesh is a mesh class for cylindrical problems
.. note::
for a cylindrically symmetric mesh use [hx, 1, hz]
::
cs, nc, npad = 20., 30, 8
hx = Utils.meshTensor([(cs,npad+10,-0.7), (cs,nc), (cs,npad,1.3)])
hz = Utils.meshTensor([(cs,npad ,-1.3), (cs,nc), (cs,npad,1.3)])
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, cartesianOrigin=None):
BaseTensorMesh.__init__(self, h, x0)
assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
if self.dim == 2:
print 'Warning, a disk mesh has not been tested thoroughly.'
cartesianOrigin = np.zeros(self.dim) if cartesianOrigin is None else cartesianOrigin
assert len(cartesianOrigin) == self.dim, "cartesianOrigin must be the same length as the dimension of the mesh."
self.cartesianOrigin = np.array(cartesianOrigin, dtype=float)
@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.identity(self.nNz), dr)
Dz = -sp.kron(dz, sp.identity(self.nCx))
A = self.area
E = self.edge
#Edge curl operator
self._edgeCurl = sdiag(1/A)*sp.vstack((Dz, Dr))*sdiag(E)
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 aveE2CCV(self):
"Construct the averaging operator on cell edges to cell centers."
if getattr(self, '_aveE2CCV', None) is None:
# The number of cell centers in each direction
n = self.vnC
if self.isSymmetric:
return self.aveE2CC
else:
raise NotImplementedError('wrapping in the averaging is not yet implemented')
return self._aveE2CCV
@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
@property
def aveF2CCV(self):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CCV', None) is None:
n = self.vnC
if self.isSymmetric:
avR = av(n[0])[:,1:]
avR[0,0] = 1.
self._aveF2CCV = sp.block_diag((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')
return self._aveF2CCV
def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None):
"""
Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
"""
assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
if locTypeTo is None:
locTypeTo = locType
if locType == 'F':
# do this three times for each component
X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x')
Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y')
Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z')
return sp.vstack((X,Y,Z))
if locType == 'E':
X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x')
Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y')
Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE)
return sp.vstack((X,Y,Z))
grid = getattr(Mrect, 'grid' + locTypeTo)
# This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
theta[theta < 0] += np.pi*2.0
r = ((grid[:,0] - self.cartesianOrigin[0])**2 + (grid[:,1] - self.cartesianOrigin[1])**2)**0.5
if locType in ['CC', 'N', 'Fz', 'Ez']:
G, proj = np.c_[r, theta, grid[:,2]], np.ones(r.size)
else:
dotMe = {
'Fx': Mrect.normals[:Mrect.nFx,:],
'Fy': Mrect.normals[Mrect.nFx:(Mrect.nFx+Mrect.nFy),:],
'Fz': Mrect.normals[-Mrect.nFz:,:],
'Ex': Mrect.tangents[:Mrect.nEx,:],
'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
'Ez': Mrect.tangents[-Mrect.nEz:,:],
}[locTypeTo]
if 'F' in locType:
normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
proj = ( normals * dotMe ).sum(axis=1)
if 'E' in locType:
tangents = np.c_[-np.sin(theta), np.cos(theta), np.zeros(theta.size)]
proj = ( tangents * dotMe ).sum(axis=1)
G = np.c_[r, theta, grid[:,2]]
interpType = locType
if interpType == 'Fy':
interpType = 'Fx'
elif interpType == 'Ex':
interpType = 'Ey'
Pc2r = self.getInterpolationMat(G, interpType)
Proj = sdiag(proj)
return Proj * Pc2r
if __name__ == '__main__':
from SimPEG import *
hx = np.r_[1,1,0.5]
hz = np.r_[2,1]
M = Mesh.CylMesh([hx, 1,hz], x0='00N')
M.plotImage(np.random.rand(M.nC), showIt=False)
M.plotGrid(centers=True, showIt=True)
Reference in New Issue View Git Blame Copy Permalink