wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-30 08:10:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of https://github.com/simpeg/simpeg into Map

Browse Source
This commit is contained in:
rowanc1
2014-04-15 11:13:21 -07:00
parent 987b759088 ca7c74e835
commit a14a0dbaa1
18 changed files with 1200 additions and 425 deletions
Download Patch File Download Diff File Expand all files Collapse all files
AUTHORS.rst
+8
View File
@@ -0,0 +1,8 @@
Credits
=======
- Rowan Cockett
- Eldad Haber
- Lindsey Heagy
- Seogi Kang
- Dave Marchant
SimPEG/Mesh/BaseMesh.py
+94 -10
View File
@@ -26,7 +26,7 @@ class BaseMesh(object):
# Ensure x0 & n are 1D vectors
self._n = np.array(n, dtype=int).ravel()
self._x0 = np.array(x0).ravel()
self._x0 = np.array(x0, dtype=float).ravel()
@property
def x0(self):
@@ -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.
SimPEG/Mesh/Cyl1DMesh.py
+1
View File
@@ -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"
SimPEG/Mesh/CylMesh.py
+293
View File
@@ -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
SimPEG/Mesh/DiffOperators.py
+6 -60
View File
@@ -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)
SimPEG/Mesh/InnerProducts.py
-1
View File
@@ -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)
SimPEG/Mesh/TensorMesh.py
+302 -337
View File
@@ -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!')
SimPEG/Mesh/TensorView.py → SimPEG/Mesh/View.py
+101
View File
@@ -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)]
SimPEG/Mesh/__init__.py
+3 -6
View File
@@ -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
SimPEG/Solver.py
+1 -1
View File
@@ -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
SimPEG/Survey.py
+1
View File
@@ -209,6 +209,7 @@ class BaseRx(object):
def nD(self):
return self.locs.shape[0]
class BaseTx(object):
"""SimPEG Transmitter Object"""
SimPEG/Tests/TestUtils.py
+16 -5
View File
@@ -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'
SimPEG/Tests/test_Solver.py
+1 -1
View File
@@ -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
SimPEG/Tests/test_basemesh.py
+1 -1
View File
@@ -1,6 +1,6 @@
import unittest
import sys
from SimPEG.Mesh import BaseRectangularMesh
from SimPEG.Mesh.BaseMesh import BaseRectangularMesh
import numpy as np
SimPEG/Tests/test_cylMesh.py
+327
View File
@@ -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()
SimPEG/Tests/test_operators.py
+1
View File
@@ -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])
SimPEG/Utils/__init__.py
+1 -1
View File
@@ -1,5 +1,5 @@
from matutils import *
from meshutils import exampleLrmGrid, meshTensors, points2nodes
from meshutils import exampleLrmGrid, meshTensors, points2nodes, writeUBCTensorMesh, writeUBCTensorModel
from lrmutils import volTetra, faceInfo, indexCube
from interputils import interpmat
from ipythonutils import easyAnimate as animate
SimPEG/Utils/meshutils.py
+43 -2
View File
@@ -24,7 +24,6 @@ def exampleLrmGrid(nC, exType):
amt[amt < 0] = 0
return [X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt]
def meshTensors(*args):
"""
**meshTensors** takes any number of tuples that have the form::
@@ -58,7 +57,7 @@ def points2nodes(mesh, pts):
Move a list of the nearest nodes to a set of points
:param simpeg.Mesh.TensorMesh mesh: The mesh
:param numpy.ndarray pts: Points to move}
:param numpy.ndarray pts: Points to move
:rtype: numpy.ndarray
:return: nodeInds
"""
@@ -76,6 +75,48 @@ def points2nodes(mesh, pts):
return nodeInds
def writeUBCTensorMesh(mesh, fileName):
"""
Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
:param simpeg.Mesh.TensorMesh mesh: The mesh
:param str fileName: File to write to
"""
assert mesh.dim == 3
s = ''
s += '%i %i %i\n' %tuple(mesh.vnC)
origin = mesh.x0
origin.dtype = float
origin[2] = origin[2]+mesh.hz.sum()
s += '%.2f %.2f %.2f\n' %tuple(origin)
s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz)
f = open(fileName, 'w')
f.write(s)
f.close()
def writeUBCTensorModel(mesh, model, fileName):
"""
Writes a model associated with a SimPEG TensorMesh
to a UBC-GIF format model file.
:param simpeg.Mesh.TensorMesh mesh: The mesh
:param numpy.ndarray model: The model
:param str fileName: File to write to
"""
# Reshape to [z,y,x]
model3D = np.reshape(model, mesh.vnC)
# Permute to [z,x,y]
model3D = np.swapaxes(model3D, 1, 2)
# Flip z to positive down
model3D = model3D[::-1,:,:]
np.savetxt(fileName, mkvc(model3D))
if __name__ == '__main__':
from SimPEG import Mesh