wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-15 11:26:09 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #53 from simpeg/develop

Browse Source 
Merge Develop into Master
This commit is contained in:
Rowan Cockett
2014-02-13 17:28:10 -08:00
parent 6c269b54dc 7f9a5075d9
commit 26c3d8d616
39 changed files with 2469 additions and 2190 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.travis.yml
+13 -2
View File
@@ -4,9 +4,20 @@ python:
virtualenv:
system_site_packages: true
before_install:
- sudo apt-get install -qq python-numpy python-scipy python-matplotlib
- python SimPEG/setup.py
- sudo apt-get install -qq gcc gfortran libblas-dev liblapack-dev python-numpy python-scipy python-matplotlib python-pip
- sudo pip install scipy --upgrade
- sudo pip install numpy --upgrade
- cd SimPEG
- python setup.py
- cd ../
# command to install dependencies
install: "pip install -r requirements.txt --use-mirrors"
# command to run tests
script: nosetests -v
notifications:
email:
- rowanc1@gmail.com
- sgkang09@gmail.com
- dwfmarchant@gmail.com
- lindseyheagy@gmail.com
SimPEG/Data.py
+23 -5
View File
@@ -4,7 +4,7 @@ import Utils, numpy as np
class BaseData(object):
"""Data holds the observed data, and the standard deviations."""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
std = None #: Estimated Standard Deviations
dobs = None #: Observed data
@@ -56,21 +56,39 @@ class BaseData(object):
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.field(m)
return Utils.mkvc(self.projectField(u))
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.projectFields(u))
@Utils.count
def projectField(self, u):
def projectFields(self, u):
"""
This function projects the fields onto the data space.
.. math::
d_\\text{pred} = P(u(m))
d_\\text{pred} = \mathbf{P} u(m)
"""
return u
@Utils.count
def projectFieldsAdjoint(self, d):
"""
This function is the adjoint of the projection.
**projectFieldsAdjoint** is used in the
calculation of the sensitivities.
.. math::
u = \mathbf{P}^\\top d
:param numpy.array d: data
:param numpy.array u: fields (ish)
:rtype: fields like object
:return: data
"""
return d
#TODO: def projectFieldDeriv(self, u): Does this need to be made??!
@Utils.count
SimPEG/Examples/DC.py
-249
View File
@@ -1,249 +0,0 @@
from SimPEG import *
class DCData(Data.BaseData):
"""
**DCData**
Geophysical DC resistivity data.
"""
P = None #: projection
def __init__(self, **kwargs):
Data.BaseData.__init__(self, **kwargs)
Utils.setKwargs(self, **kwargs)
def reshapeFields(self, u):
if len(u.shape) == 1:
u = u.reshape([-1, self.RHS.shape[1]], order='F')
return u
def projectField(self, u):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
"""
u = self.reshapeFields(u)
return Utils.mkvc(self.P*u)
class DCProblem(Problem.BaseProblem):
"""
**DCProblem**
Geophysical DC resistivity problem.
"""
dataPair = DCData
def __init__(self, mesh, model, **kwargs):
Problem.BaseProblem.__init__(self, mesh, model)
self.mesh.setCellGradBC('neumann')
Utils.setKwargs(self, **kwargs)
def createMatrix(self, m):
"""
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.csc_matrix
:return: A(m)
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
Where M() is the mass matrix and mT is the model transform.
"""
D = self.mesh.faceDiv
G = self.mesh.cellGrad
sigma = self.model.transform(m)
Msig = self.mesh.getFaceMass(sigma)
A = D*Msig*G
return A.tocsc()
def field(self, m):
A = self.createMatrix(m)
solve = Solver(A)
phi = solve.solve(self.data.RHS)
return Utils.mkvc(phi)
def J(self, m, v, u=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: Jv
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
\\nabla_u (A(m)u - q) = A(m)
\\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
Where M() is the mass matrix and mT is the model transform.
.. math::
J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
"""
if u is None:
u = self.field(m)
u = self.data.reshapeFields(u)
P = self.data.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
Av_dm = self.mesh.getFaceMassDeriv()
mT_dm = self.model.transformDeriv(m)
dCdu = A
dCdm = np.empty_like(u)
for i, ui in enumerate(u.T): # loop over each column
dCdm[:, i] = D * ( Utils.sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) )
solve = Solver(dCdu)
Jv = - P * solve.solve(dCdm)
return Utils.mkvc(Jv)
def Jt(self, m, v, u=None):
"""Takes data, turns it into a model..ish"""
if u is None:
u = self.field(m)
u = self.data.reshapeFields(u)
v = self.data.reshapeFields(v)
P = self.data.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
Av_dm = self.mesh.getFaceMassDeriv()
mT_dm = self.model.transformDeriv(m)
dCdu = A.T
solve = Solver(dCdu)
w = solve.solve(P.T*v)
Jtv = 0
for i, ui in enumerate(u.T): # loop over each column
Jtv += Utils.sdiag( G * ui ) * ( D.T * w[:,i] )
Jtv = - mT_dm.T * ( Av_dm.T * Jtv )
return Jtv
def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
""" Generate projection matrix (Q) and """
elecend = 0.5+spacelec*(nelec-1)
elecLocR = np.linspace(elecini, elecend, nelec)
elecLocT = elecLocR+1
nrx = nelec-1
ntx = nelec-1
q = np.zeros((mesh.nC, ntx))
Q = np.zeros((mesh.nC, nrx))
for i in range(nrx):
rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
q[txind1,i] = 1
q[txind2,i] = -1
Q[rxind1,i] = 1
Q[rxind2,i] = -1
Q = sp.csr_matrix(Q)
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
return q, Q, rxmidLoc
if __name__ == '__main__':
import matplotlib.pyplot as plt
# Create the mesh
h1 = np.ones(20)
h2 = np.ones(100)
M = Mesh.TensorMesh([h1,h2])
# Create some parameters for the model
sig1 = np.log(1)
sig2 = np.log(0.01)
# Create a synthetic model from a block in a half-space
p0 = [5, 10]
p1 = [15, 50]
condVals = [sig1, sig2]
mSynth = Utils.ModelBuilder.defineBlockConductivity(M.gridCC,p0,p1,condVals)
plt.colorbar(M.plotImage(mSynth))
# plt.show()
# Set up the projection
nelec = 50
spacelec = 2
surfloc = 0.5
elecini = 0.5
elecend = 0.5+spacelec*(nelec-1)
elecLocR = np.linspace(elecini, elecend, nelec)
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, M)
P = Q.T
model = Model.LogModel(M)
prob = DCProblem(M, model)
# Create some data
data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=q)
u = prob.field(mSynth)
u = data.reshapeFields(u)
M.plotImage(u[:,10])
plt.show()
# Now set up the prob to do some minimization
# prob.dobs = dobs
# prob.std = dobs*0 + 0.05
m0 = M.gridCC[:,0]*0+sig2
reg = Regularization.Tikhonov(model)
objFunc = ObjFunction.BaseObjFunction(data, reg)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
inv = Inversion.BaseInversion(objFunc, opt)
# Check Derivative
derChk = lambda m: [objFunc.dataObj(m), objFunc.dataObjDeriv(m)]
# Tests.checkDerivative(derChk, mSynth)
print objFunc.dataObj(m0)
print objFunc.dataObj(mSynth)
m = inv.run(m0)
plt.colorbar(M.plotImage(m))
print m
plt.show()
SimPEG/Examples/__init__.py
-2
View File
@@ -1,2 +0,0 @@
import DC
import Linear
SimPEG/Inversion.py
+1 -1
View File
@@ -7,7 +7,7 @@ class BaseInversion(object):
"""BaseInversion(objFunc, opt, **kwargs)
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
name = 'BaseInversion'
SimPEG/Mesh/BaseMesh.py
+120 -188
View File
@@ -27,18 +27,16 @@ class BaseMesh(object):
# Ensure x0 & n are 1D vectors
self._n = np.array(n, dtype=int).ravel()
self._x0 = np.array(x0).ravel()
self._dim = len(n)
def x0():
doc = """
@property
def x0(self):
"""
Origin of the mesh
:rtype: numpy.array (dim, )
:return: x0
"""
fget = lambda self: self._x0
return locals()
x0 = property(**x0())
return self._x0
def r(self, x, xType='CC', outType='CC', format='V'):
"""
@@ -147,61 +145,57 @@ class BaseMesh(object):
else:
return switchKernal(x)
def dim():
doc = """
@property
def dim(self):
"""
The dimension of the mesh (1, 2, or 3).
:rtype: int
:return: dim
"""
fget = lambda self: self._dim
return locals()
dim = property(**dim())
return len(self._n)
def nCx():
doc = """
@property
def nCx(self):
"""
Number of cells in the x direction
:rtype: int
:return: nCx
"""
fget = lambda self: self._n[0]
return locals()
nCx = property(**nCx())
return self._n[0]
def nCy():
doc = """
@property
def nCy(self):
"""
Number of cells in the y direction
:rtype: int
:return: nCy or None if dim < 2
"""
return None if self.dim < 2 else self._n[1]
def fget(self):
if self.dim > 1:
return self._n[1]
else:
return None
return locals()
nCy = property(**nCy())
def nCz():
doc = """Number of cells in the z direction
@property
def nCz(self):
"""Number of cells in the z direction
:rtype: int
:return: nCz or None if dim < 3
"""
return None if self.dim < 3 else self._n[2]
def fget(self):
if self.dim > 2:
return self._n[2]
else:
return None
return locals()
nCz = property(**nCz())
@property
def nCv(self):
"""
Total number of cells in each direction
def nC():
:rtype: numpy.array (dim, )
:return: [nCx, nCy, nCz]
"""
return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
@property
def nC(self):
doc = """
Total number of cells in the model.
@@ -214,65 +208,50 @@ class BaseMesh(object):
from SimPEG import Mesh, np
Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True)
"""
fget = lambda self: np.prod(self._n)
return locals()
nC = property(**nC())
return self.nCv.prod()
def nCv():
doc = """
Total number of cells in each direction
:rtype: numpy.array (dim, )
:return: [nCx, nCy, nCz]
@property
def nNx(self):
"""
fget = lambda self: np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
return locals()
nCv = property(**nCv())
def nNx():
doc = """
Number of nodes in the x-direction
:rtype: int
:return: nNx
"""
fget = lambda self: self.nCx + 1
return locals()
nNx = property(**nNx())
return self.nCx + 1
def nNy():
doc = """
@property
def nNy(self):
"""
Number of noes in the y-direction
:rtype: int
:return: nNy or None if dim < 2
"""
return None if self.dim < 2 else self.nCy + 1
def fget(self):
if self.dim > 1:
return self._n[1] + 1
else:
return None
return locals()
nNy = property(**nNy())
def nNz():
doc = """
@property
def nNz(self):
"""
Number of nodes in the z-direction
:rtype: int
:return: nNz or None if dim < 3
"""
return None if self.dim < 3 else self.nCz + 1
def fget(self):
if self.dim > 2:
return self._n[2] + 1
else:
return None
return locals()
nNz = property(**nNz())
@property
def nNv(self):
"""
Total number of nodes in each direction
def nN():
:rtype: numpy.array (dim, )
:return: [nNx, nNy, nNz]
"""
return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
@property
def nN(self):
doc = """
Total number of nodes
@@ -285,66 +264,41 @@ class BaseMesh(object):
from SimPEG import Mesh, np
Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True)
"""
fget = lambda self: np.prod(self.nCv + 1)
return locals()
nN = property(**nN())
return self.nNv.prod()
def nNv():
doc = """
Total number of nodes in each direction
:rtype: numpy.array (dim, )
:return: [nNx, nNy, nNz]
@property
def nEx(self):
"""
fget = lambda self: np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
return locals()
nNv = property(**nNv())
def nEx():
doc = """
Number of x-edges in each direction
:rtype: numpy.array (dim, )
:return: nEx
"""
fget = lambda self: np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
return locals()
nEx = property(**nEx())
return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
def nEy():
doc = """
@property
def nEy(self):
"""
Number of y-edges in each direction
:rtype: numpy.array (dim, )
:return: nEy or None if dim < 2
"""
return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
def fget(self):
if self.dim > 1:
return np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
else:
return None
return locals()
nEy = property(**nEy())
def nEz():
doc = """
@property
def nEz(self):
"""
Number of z-edges in each direction
:rtype: numpy.array (dim, )
:return: nEz or None if dim < 3
"""
return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
def fget(self):
if self.dim > 2:
return np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
else:
return None
return locals()
nEz = property(**nEz())
def nEv():
doc = """
@property
def nEv(self):
"""
Total number of edges in each direction
:rtype: numpy.array (dim, )
@@ -356,67 +310,53 @@ class BaseMesh(object):
from SimPEG import Mesh, np
Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True)
"""
fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
return locals()
nEv = property(**nEv())
return np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
def nE():
doc = """
@property
def nE(self):
"""
Total number of edges.
:rtype: int
:return: sum([prod(nEx), prod(nEy), prod(nEz)])
"""
fget = lambda self: np.sum(self.nEv)
return locals()
nE = property(**nE())
return self.nEv.sum()
def nFx():
doc = """
@property
def nFx(self):
"""
Number of x-faces in each direction
:rtype: numpy.array (dim, )
:return: nFx
"""
fget = lambda self: np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
return locals()
nFx = property(**nFx())
return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
def nFy():
doc = """
@property
def nFy(self):
"""
Number of y-faces in each direction
:rtype: numpy.array (dim, )
:return: nFy or None if dim < 2
"""
return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
def fget(self):
if self.dim > 1:
return np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
else:
return None
return locals()
nFy = property(**nFy())
def nFz():
doc = """
@property
def nFz(self):
"""
Number of z-faces in each direction
:rtype: numpy.array (dim, )
:return: nFz or None if dim < 3
"""
return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
def fget(self):
if self.dim > 2:
return np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
else:
return None
return locals()
nFz = property(**nFz())
def nFv():
doc = """
@property
def nFv(self):
"""
Total number of faces in each direction
:rtype: numpy.array (dim, )
@@ -428,64 +368,56 @@ class BaseMesh(object):
from SimPEG import Mesh, np
Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True)
"""
fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
return locals()
nFv = property(**nFv())
return np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
def nF():
doc = """
@property
def nF(self):
"""
Total number of faces.
:rtype: int
:return: sum([prod(nFx), prod(nFy), prod(nFz)])
:return: sum([nFx, nFy, nFz])
"""
fget = lambda self: np.sum(self.nFv)
return locals()
nF = property(**nF())
return self.nFv.sum()
def normals():
doc = """
@property
def normals(self):
"""
Face Normals
:rtype: numpy.array (sum(nF), dim)
:return: normals
"""
if self.dim == 2:
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])]
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])]
return np.r_[nX, nY]
elif self.dim == 3:
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])]
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])]
nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])]
return np.r_[nX, nY, nZ]
def fget(self):
if self.dim == 2:
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])]
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])]
return np.r_[nX, nY]
elif self.dim == 3:
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])]
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])]
nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])]
return np.r_[nX, nY, nZ]
return locals()
normals = property(**normals())
def tangents():
doc = """
@property
def tangents(self):
"""
Edge Tangents
:rtype: numpy.array (sum(nE), dim)
:return: normals
"""
if self.dim == 2:
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])]
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])]
return np.r_[tX, tY]
elif self.dim == 3:
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])]
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])]
tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])]
return np.r_[tX, tY, tZ]
def fget(self):
if self.dim == 2:
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])]
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])]
return np.r_[tX, tY]
elif self.dim == 3:
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])]
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])]
tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])]
return np.r_[tX, tY, tZ]
return locals()
tangents = property(**tangents())
def projectFaceVector(self, fV):
"""
SimPEG/Mesh/Cyl1DMesh.py
+12
View File
@@ -37,6 +37,9 @@ class Cyl1DMesh(object):
return locals()
h = property(**h())
@property
def dim(self): return 2
def z0():
doc = "The z-origin"
def fget(self):
@@ -290,6 +293,15 @@ class Cyl1DMesh(object):
_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
####################################################
SimPEG/Mesh/DiffOperators.py
+126 -111
View File
@@ -462,126 +462,137 @@ class DiffOperators(object):
# --------------- Averaging ---------------------
def aveF2CC():
doc = "Construct the averaging operator on cell faces to cell centers."
@property
def aveF2CC(self):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CC', None) is None:
n = self.nCv
if(self.dim == 1):
self._aveF2CC = av(n[0])
elif(self.dim == 2):
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
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
def fget(self):
if(self._aveF2CC is None):
n = self.nCv
if(self.dim == 1):
self._aveF2CC = av(n[0])
elif(self.dim == 2):
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
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
return locals()
_aveF2CC = None
aveF2CC = property(**aveF2CC())
def aveCC2F():
doc = "Construct the averaging operator on cell cell centers to faces."
@property
def aveF2CCV(self):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CCV', None) is None:
n = self.nCv
if(self.dim == 1):
self._aveF2CCV = av(n[0])
elif(self.dim == 2):
self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
self._aveF2CCV = sp.block_diag((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._aveF2CCV
def fget(self):
if(self._aveCC2F is None):
n = self.nCv
if(self.dim == 1):
self._aveCC2F = avExtrap(n[0])
elif(self.dim == 2):
self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
return self._aveCC2F
return locals()
_aveCC2F = None
aveCC2F = property(**aveCC2F())
@property
def aveCC2F(self):
"Construct the averaging operator on cell cell centers to faces."
if getattr(self, '_aveCC2F', None) is None:
n = self.nCv
if(self.dim == 1):
self._aveCC2F = avExtrap(n[0])
elif(self.dim == 2):
self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
return self._aveCC2F
def aveE2CC():
doc = "Construct the averaging operator on cell edges to cell centers."
@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.nCv
if(self.dim == 1):
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")
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")
return self._aveE2CC
def fget(self):
if(self._aveE2CC is None):
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
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")
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")
return self._aveE2CC
return locals()
_aveE2CC = None
aveE2CC = property(**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.nCv
if(self.dim == 1):
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
elif(self.dim == 2):
self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])),
sp.kron(speye(n[1]), av(n[0]))), format="csr")
elif(self.dim == 3):
self._aveE2CCV = sp.block_diag((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._aveE2CCV
def aveN2CC():
doc = "Construct the averaging operator on cell nodes to cell centers."
@property
def aveN2CC(self):
"Construct the averaging operator on cell nodes to cell centers."
if getattr(self, '_aveN2CC', None) is None:
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
self._aveN2CC = av(n[0])
elif(self.dim == 2):
self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr()
elif(self.dim == 3):
self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr()
return self._aveN2CC
def fget(self):
if(self._aveN2CC is None):
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
self._aveN2CC = av(n[0])
elif(self.dim == 2):
self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr()
elif(self.dim == 3):
self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr()
return self._aveN2CC
return locals()
_aveN2CC = None
aveN2CC = property(**aveN2CC())
@property
def aveN2E(self):
"Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."
def aveN2E():
doc = "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."
if getattr(self, '_aveN2E', None) is None:
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
self._aveN2E = av(n[0])
elif(self.dim == 2):
self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
elif(self.dim == 3):
self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
return self._aveN2E
def fget(self):
if(self._aveN2E is None):
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
self._aveN2E = av(n[0])
elif(self.dim == 2):
self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
elif(self.dim == 3):
self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
return self._aveN2E
return locals()
_aveN2E = None
aveN2E = property(**aveN2E())
def aveN2F():
doc = "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
def fget(self):
if(self._aveN2F is None):
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
self._aveN2F = av(n[0])
elif(self.dim == 2):
self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
elif(self.dim == 3):
self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
return self._aveN2F
return locals()
_aveN2F = None
aveN2F = property(**aveN2F())
@property
def aveN2F(self):
"Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
if getattr(self, '_aveN2F', None) is None:
# The number of cell centers in each direction
n = self.nCv
if(self.dim == 1):
self._aveN2F = av(n[0])
elif(self.dim == 2):
self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
elif(self.dim == 3):
self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
return self._aveN2F
# --------------- Methods ---------------------
@@ -633,3 +644,7 @@ class DiffOperators(object):
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
+388 -379
View File
@@ -78,29 +78,208 @@ class InnerProducts(object):
def __init__(self):
raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.')
def getFaceInnerProduct(self, mu=None, returnP=False):
"""Wrapper function,
:py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct`
:py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D`
def getFaceInnerProduct(M, mu=None, returnP=False):
"""
if self.dim == 2:
return getFaceInnerProduct2D(self, mu, returnP)
elif self.dim == 3:
return getFaceInnerProduct(self, mu, returnP)
:param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (sum(nF), sum(nF))
def getEdgeInnerProduct(self, sigma=None, returnP=False):
"""Wrapper function,
Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
:py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct`
.. math::
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right]
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right]
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right]
\mathbf{M}(\\vec{\mu}) = {1\over 8}
\left(\sum_{i=1}^8
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P000, P100, P010, P110, P001, P101, P011, P111]
Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
**For 2D:**
Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows:
.. math::
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right]
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right]
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right]
.. math::
\mathbf{M}(\\vec{\mu}) = {1\over 4}
\left(\sum_{i=1}^4
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P00, P10, P01, P11]
Here each P (2*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
:py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D`
"""
if self.dim == 2:
return getEdgeInnerProduct2D(self, sigma, returnP)
elif self.dim == 3:
return getEdgeInnerProduct(self, sigma, returnP)
if M.dim == 2:
# Square root of cell volume multiplied by 1/4
v = np.sqrt(0.25*M.vol)
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
Pxx = _getFacePxx(M)
P000 = V2*Pxx('fXm', 'fYm')
P100 = V2*Pxx('fXp', 'fYm')
P010 = V2*Pxx('fXm', 'fYp')
P110 = V2*Pxx('fXp', 'fYp')
elif M.dim == 3:
# Square root of cell volume multiplied by 1/8
v = np.sqrt(0.125*M.vol)
V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry
Pxxx = _getFacePxxx(M)
P000 = V3*Pxxx('fXm', 'fYm', 'fZm')
P100 = V3*Pxxx('fXp', 'fYm', 'fZm')
P010 = V3*Pxxx('fXm', 'fYp', 'fZm')
P110 = V3*Pxxx('fXp', 'fYp', 'fZm')
P001 = V3*Pxxx('fXm', 'fYm', 'fZp')
P101 = V3*Pxxx('fXp', 'fYm', 'fZp')
P011 = V3*Pxxx('fXm', 'fYp', 'fZp')
P111 = V3*Pxxx('fXp', 'fYp', 'fZp')
Mu = _makeTensor(M, mu)
A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
P = [P000, P100, P010, P110]
if M.dim == 3:
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
P += [P001, P101, P011, P111]
if returnP:
return A, P
else:
return A
def getEdgeInnerProduct(M, sigma=None, returnP=False):
"""
:param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (sum(nE), sum(nE))
Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows:
.. math::
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right]
What is returned:
.. math::
\mathbf{M}(\Sigma) = {1\over 8}
\left(\sum_{i=1}^8
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P000, P100, P010, P110, P001, P101, P011, P111]
Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
**For 2D:**
Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows:
.. math::
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right]
.. math::
\mathbf{M}(\Sigma) = {1\over 4}
\left(\sum_{i=1}^4
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P00, P10, P01, P11]
Here each P (2*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
"""
# We will multiply by V on each side to keep symmetry
if M.dim == 2:
# Square root of cell volume multiplied by 1/4
v = np.sqrt(0.25*M.vol)
V = sdiag(np.r_[v, v])
eP = _getEdgePxx(M)
P000 = V*eP('eX0', 'eY0')
P100 = V*eP('eX0', 'eY1')
P010 = V*eP('eX1', 'eY0')
P110 = V*eP('eX1', 'eY1')
elif M.dim == 3:
# Square root of cell volume multiplied by 1/8
v = np.sqrt(0.125*M.vol)
V = sdiag(np.r_[v, v, v])
eP = _getEdgePxxx(M)
P000 = V*eP('eX0', 'eY0', 'eZ0')
P100 = V*eP('eX0', 'eY1', 'eZ1')
P010 = V*eP('eX1', 'eY0', 'eZ2')
P110 = V*eP('eX1', 'eY1', 'eZ3')
P001 = V*eP('eX2', 'eY2', 'eZ0')
P101 = V*eP('eX2', 'eY3', 'eZ1')
P011 = V*eP('eX3', 'eY2', 'eZ2')
P111 = V*eP('eX3', 'eY3', 'eZ3')
Sigma = _makeTensor(M, sigma)
A = P000.T*Sigma*P000 + P100.T*Sigma*P100 + P010.T*Sigma*P010 + P110.T*Sigma*P110
P = [P000, P100, P010, P110]
if M.dim == 3:
A = A + P001.T*Sigma*P001 + P101.T*Sigma*P101 + P011.T*Sigma*P011 + P111.T*Sigma*P111
P += [P001, P101, P011, P111]
if returnP:
return A, P
else:
return A
# ------------------------ Geometries ------------------------------
#
@@ -121,434 +300,264 @@ class InnerProducts(object):
# | |/
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
def _makeTensor(M, sigma):
if sigma is None: # default is ones
sigma = np.ones((M.nC, 1))
def getFaceInnerProduct(mesh, mu=None, returnP=False):
"""
:param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (sum(nF), sum(nF))
if M.dim == 2:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma])
elif sigma.shape[1] == 2: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
elif sigma.shape[1] == 3: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
Sigma = sp.vstack((row1, row2))
elif M.dim == 3:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma, sigma])
elif sigma.shape[1] == 3: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
elif sigma.shape[1] == 6: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
Sigma = sp.vstack((row1, row2, row3))
return Sigma
Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
def _getFacePxx(M):
if M._meshType == 'TREE':
return M._getFacePxx
.. math::
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right]
return _getFacePxx_Rectangular(M)
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right]
def _getFacePxxx(M):
if M._meshType == 'TREE':
return M._getFacePxxx
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right]
return _getFacePxxx_Rectangular(M)
\mathbf{M}(\\vec{\mu}) = {1\over 8}
\left(\sum_{i=1}^8
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
def _getEdgePxx(M):
if M._meshType == 'TREE':
return M._getEdgePxx
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
return _getEdgePxx_Rectangular(M)
P = [P000, P001, P010, P011, P100, P101, P110, P111]
def _getEdgePxxx(M):
if M._meshType == 'TREE':
return M._getEdgePxxx
Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
return _getEdgePxxx_Rectangular(M)
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
def _getFacePxx_Rectangular(M):
"""returns a function for creating projection matrices
Note that this is completed for each cell in the mesh at the same time.
Mats takes you from faces a subset of all faces on only the
faces that you ask for.
"""
These are centered around a single nodes.
if mu is None: # default is ones
mu = np.ones((mesh.nC, 1))
For example, if this was your entire mesh:
m = np.array([mesh.nCx, mesh.nCy, mesh.nCz])
nc = mesh.nC
f3(Yp)
2_______________3
| |
| |
| |
f0(Xm) | x | f1(Xp)
| |
| |
|_______________|
0 1
f2(Ym)
i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2]))
Pxx('m','m') = | 1, 0, 0, 0 |
| 0, 0, 1, 0 |
iijjkk = ndgrid(i, j, k)
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
Pxx('p','m') = | 0, 1, 0, 0 |
| 0, 0, 1, 0 |
if mesh._meshType == 'LOM':
fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M')
fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M')
fN3 = mesh.r(mesh.normals, 'F', 'Fz', 'M')
def Pxxx(pos):
ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nFv[0]
ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nFv[0] + mesh.nFv[1]
IND = np.r_[ind1, ind2, ind3].flatten()
PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr()
if mesh._meshType == 'LOM':
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]),
getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]),
getSubArray(fN3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]))
PXXX = I3x3 * PXXX
return PXXX
# no | node | f1 | f2 | f3
# 000 | i ,j ,k | i , j, k | i, j , k | i, j, k
# 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k
# 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k
# 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k
# 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1
# 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1
# 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1
# 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1
# Square root of cell volume multiplied by 1/8
v = np.sqrt(0.125*mesh.vol)
V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry
P000 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
P100 = V3*Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]])
P010 = V3*Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
P110 = V3*Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
P001 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]])
P101 = V3*Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]])
P011 = V3*Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]])
P111 = V3*Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
if mu.size == mesh.nC: # Isotropic!
mu = mkvc(mu) # ensure it is a vector.
Mu = sdiag(np.r_[mu, mu, mu])
elif mu.shape[1] == 3: # Diagonal tensor
Mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]])
elif mu.shape[1] == 6: # Fully anisotropic
row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4])))
row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5])))
row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2])))
Mu = sp.vstack((row1, row2, row3))
A = P000.T*Mu*P000 + P001.T*Mu*P001 + P010.T*Mu*P010 + P011.T*Mu*P011 + P100.T*Mu*P100 + P101.T*Mu*P101 + P110.T*Mu*P110 + P111.T*Mu*P111
P = [P000, P001, P010, P011, P100, P101, P110, P111]
if returnP:
return A, P
else:
return A
def getFaceInnerProduct2D(mesh, mu=None, returnP=False):
"""
:param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3))
:param bool returnP: returns the projection matrices
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (sum(nF), sum(nF))
Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows:
.. math::
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right]
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right]
\\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right]
.. math::
\mathbf{M}(\\vec{\mu}) = {1\over 4}
\left(\sum_{i=1}^4
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P00, P10, P01, P11]
Here each P (2*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
"""
if mu is None: # default is ones
mu = np.ones((mesh.nC, 1))
m = np.array([mesh.nCx, mesh.nCy])
nc = mesh.nC
i, j = np.int64(range(m[0])), np.int64(range(m[1]))
"""
i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
iijj = ndgrid(i, j)
ii, jj = iijj[:, 0], iijj[:, 1]
if mesh._meshType == 'LOM':
fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M')
fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M')
if M._meshType == 'LOM':
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
def Pxx(pos):
ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]])
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nFv[0]
def Pxx(xFace, yFace):
"""
xFace is 'fXp' or 'fXm'
yFace is 'fYp' or 'fYm'
"""
# no | node | f1 | f2
# 00 | i ,j | i , j | i, j
# 10 | i+1,j | i+1, j | i, j
# 01 | i ,j+1 | i , j | i, j+1
# 11 | i+1,j+1 | i+1, j | i, j+1
posFx = 0 if xFace == 'fXm' else 1
posFy = 0 if yFace == 'fYm' else 1
ind1 = sub2ind(M.nFx, np.c_[ii + posFx, jj])
ind2 = sub2ind(M.nFy, np.c_[ii, jj + posFy]) + M.nFv[0]
IND = np.r_[ind1, ind2].flatten()
PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nF))).tocsr()
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nF)))
if mesh._meshType == 'LOM':
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1]]),
getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1]]))
if M._meshType == 'LOM':
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
PXX = I2x2 * PXX
return PXX
# no | node | f1 | f2
# 00 | i ,j | i , j | i, j
# 10 | i+1,j | i+1, j | i, j
# 01 | i ,j+1 | i , j | i, j+1
# 11 | i+1,j+1 | i+1, j | i, j+1
return Pxx
# Square root of cell volume multiplied by 1/4
v = np.sqrt(0.25*mesh.vol)
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
def _getFacePxxx_Rectangular(M):
"""returns a function for creating projection matrices
P00 = V2*Pxx([[0, 0], [0, 0]])
P10 = V2*Pxx([[1, 0], [0, 0]])
P01 = V2*Pxx([[0, 0], [0, 1]])
P11 = V2*Pxx([[1, 0], [0, 1]])
Mats takes you from faces a subset of all faces on only the
faces that you ask for.
if mu.size == mesh.nC: # Isotropic!
mu = mkvc(mu) # ensure it is a vector.
Mu = sdiag(np.r_[mu, mu])
elif mu.shape[1] == 2: # Diagonal tensor
Mu = sdiag(np.r_[mu[:, 0], mu[:, 1]])
elif mu.shape[1] == 3: # Fully anisotropic
row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 2])))
row2 = sp.hstack((sdiag(mu[:, 2]), sdiag(mu[:, 1])))
Mu = sp.vstack((row1, row2))
A = P00.T*Mu*P00 + P10.T*Mu*P10 + P01.T*Mu*P01 + P11.T*Mu*P11
P = [P00, P10, P01, P11]
if returnP:
return A, P
else:
return A
def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
"""
:param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (sum(nE), sum(nE))
Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows:
.. math::
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right]
What is returned:
.. math::
\mathbf{M}(\Sigma) = {1\over 8}
\left(\sum_{i=1}^8
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P000, P001, P010, P011, P100, P101, P110, P111]
Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
These are centered around a single nodes.
"""
if sigma is None: # default is ones
sigma = np.ones((mesh.nC, 1))
m = np.array([mesh.nCx, mesh.nCy, mesh.nCz])
nc = mesh.nC
i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2]))
i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz))
iijjkk = ndgrid(i, j, k)
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
if mesh._meshType == 'LOM':
eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M')
eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M')
eT3 = mesh.r(mesh.tangents, 'E', 'Ez', 'M')
if M._meshType == 'LOM':
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
fN3 = M.r(M.normals, 'F', 'Fz', 'M')
def Pxxx(pos):
ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nEv[0]
ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nEv[0] + mesh.nEv[1]
def Pxxx(xFace, yFace, zFace):
"""
xFace is 'fXp' or 'fXm'
yFace is 'fYp' or 'fYm'
zFace is 'fZp' or 'fZm'
"""
# no | node | f1 | f2 | f3
# 000 | i ,j ,k | i , j, k | i, j , k | i, j, k
# 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k
# 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k
# 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k
# 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1
# 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1
# 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1
# 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1
posX = 0 if xFace == 'fXm' else 1
posY = 0 if yFace == 'fYm' else 1
posZ = 0 if zFace == 'fZm' else 1
ind1 = sub2ind(M.nFx, np.c_[ii + posX, jj, kk])
ind2 = sub2ind(M.nFy, np.c_[ii, jj + posY, kk]) + M.nFv[0]
ind3 = sub2ind(M.nFz, np.c_[ii, jj, kk + posZ]) + M.nFv[0] + M.nFv[1]
IND = np.r_[ind1, ind2, ind3].flatten()
PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr()
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nF))).tocsr()
if mesh._meshType == 'LOM':
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]),
getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]),
getSubArray(eT3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]))
if M._meshType == 'LOM':
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
PXXX = I3x3 * PXXX
return PXXX
return Pxxx
# no | node | e1 | e2 | e3
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
# Square root of cell volume multiplied by 1/8
v = np.sqrt(0.125*mesh.vol)
V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry
P000 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
P100 = V3*Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]])
P010 = V3*Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]])
P110 = V3*Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]])
P001 = V3*Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]])
P101 = V3*Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]])
P011 = V3*Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]])
P111 = V3*Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
if sigma.size == mesh.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma, sigma])
elif sigma.shape[1] == 3: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
elif sigma.shape[1] == 6: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
Sigma = sp.vstack((row1, row2, row3))
A = P000.T*Sigma*P000 + P001.T*Sigma*P001 + P010.T*Sigma*P010 + P011.T*Sigma*P011 + P100.T*Sigma*P100 + P101.T*Sigma*P101 + P110.T*Sigma*P110 + P111.T*Sigma*P111
P = [P000, P001, P010, P011, P100, P101, P110, P111]
if returnP:
return A, P
else:
return A
def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False):
"""
:param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3))
:param bool returnP: returns the projection matrices
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (sum(nE), sum(nE))
Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows:
.. math::
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right]
\Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right]
.. math::
\mathbf{M}(\Sigma) = {1\over 4}
\left(\sum_{i=1}^4
\mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c
\\right)
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P00, P10, P01, P11]
Here each P (2*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
.. math::
\mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
Note that this is completed for each cell in the mesh at the same time.
"""
if sigma is None: # default is ones
sigma = np.ones((mesh.nC, 1))
m = np.array([mesh.nCx, mesh.nCy])
nc = mesh.nC
i, j = np.int64(range(m[0])), np.int64(range(m[1]))
def _getEdgePxx_Rectangular(M):
i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
iijj = ndgrid(i, j)
ii, jj = iijj[:, 0], iijj[:, 1]
if mesh._meshType == 'LOM':
eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M')
eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M')
if M._meshType == 'LOM':
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
def Pxx(pos):
ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]])
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nEv[0]
def Pxx(xEdge, yEdge):
# no | node | e1 | e2
# 00 | i ,j | i ,j | i ,j
# 10 | i+1,j | i ,j | i+1,j
# 01 | i ,j+1 | i ,j+1 | i ,j
# 11 | i+1,j+1 | i ,j+1 | i+1,j
posX = 0 if xEdge == 'eX0' else 1
posY = 0 if yEdge == 'eY0' else 1
ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX])
ind2 = sub2ind(M.nEy, np.c_[ii + posY, jj]) + M.nEv[0]
IND = np.r_[ind1, ind2].flatten()
PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nE))).tocsr()
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nE))).tocsr()
if mesh._meshType == 'LOM':
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1]]),
getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1]]))
if M._meshType == 'LOM':
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
PXX = I2x2 * PXX
return PXX
return Pxx
# no | node | e1 | e2
# 00 | i ,j | i ,j | i ,j
# 10 | i+1,j | i ,j | i+1,j
# 01 | i ,j+1 | i ,j+1 | i ,j
# 11 | i+1,j+1 | i ,j+1 | i+1,j
def _getEdgePxxx_Rectangular(M):
i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz))
# Square root of cell volume multiplied by 1/4
v = np.sqrt(0.25*mesh.vol)
V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry
iijjkk = ndgrid(i, j, k)
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
P00 = V2*Pxx([[0, 0], [0, 0]])
P10 = V2*Pxx([[0, 0], [1, 0]])
P01 = V2*Pxx([[0, 1], [0, 0]])
P11 = V2*Pxx([[0, 1], [1, 0]])
if M._meshType == 'LOM':
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
if sigma.size == mesh.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma])
elif sigma.shape[1] == 2: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
elif sigma.shape[1] == 3: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
Sigma = sp.vstack((row1, row2))
def Pxxx(xEdge, yEdge, zEdge):
A = P00.T*Sigma*P00 + P10.T*Sigma*P10 + P01.T*Sigma*P01 + P11.T*Sigma*P11
P = [P00, P10, P01, P11]
if returnP:
return A, P
else:
return A
# no | node | e1 | e2 | e3
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1]
posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1]
posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1]
ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX[0], kk + posX[1]])
ind2 = sub2ind(M.nEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEv[0]
ind3 = sub2ind(M.nEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEv[0] + M.nEv[1]
IND = np.r_[ind1, ind2, ind3].flatten()
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nE))).tocsr()
if M._meshType == 'LOM':
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
PXXX = I3x3 * PXXX
return PXXX
return Pxxx
if __name__ == '__main__':
from TensorMesh import TensorMesh
h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])]
mesh = TensorMesh(h)
mu = np.ones((mesh.nC, 6))
A, P = mesh.getFaceInnerProduct(mu, returnP=True)
B, P = mesh.getEdgeInnerProduct(mu, returnP=True)
M = TensorMesh(h)
mu = np.ones((M.nC, 6))
A, P = M.getFaceInnerProduct(mu, returnP=True)
B, P = M.getEdgeInnerProduct(mu, returnP=True)
SimPEG/Mesh/LogicallyOrthogonalMesh.py
+1 -1
View File
@@ -26,7 +26,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
M.plotGrid(showIt=True)
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
_meshType = 'LOM'
SimPEG/Mesh/TensorMesh.py
+1 -1
View File
@@ -33,7 +33,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
_meshType = 'TENSOR'
SimPEG/Mesh/TensorView.py
+1 -1
View File
@@ -404,7 +404,7 @@ class TensorView(object):
::
def function(var, ax, clim, tlt, i):
tlt.set_text('%%d'%%i)
tlt.set_text('%d'%i)
return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim)
mesh.video([model1, model2, ..., modeln],function)
SimPEG/Mesh/TreeMesh.py
+1116
View File
File diff suppressed because it is too large Load Diff
SimPEG/Mesh/__init__.py
+1
View File
@@ -1,5 +1,6 @@
from Cyl1DMesh import Cyl1DMesh
from TensorMesh import TensorMesh
from TreeMesh import TreeMesh
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
from BaseMesh import BaseMesh
from TensorView import TensorView
SimPEG/Model.py
+86 -4
View File
@@ -1,5 +1,5 @@
import Utils, Parameters, numpy as np, scipy.sparse as sp
from Tests import checkDerivative
class BaseModel(object):
"""
@@ -7,7 +7,7 @@ class BaseModel(object):
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
mesh = None #: A SimPEG Mesh
@@ -55,9 +55,14 @@ class BaseModel(object):
"""Number of parameters in the model."""
return self.mesh.nC
def example(self, modelType=None):
return np.random.rand(self.mesh.nC)
def example(self):
return np.random.rand(self.nP)
def test(self, m=None):
print 'Testing the %s Class!' % self.__class__.__name__
if m is None:
m = self.example()
return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, plotIt=False)
class LogModel(BaseModel):
@@ -127,3 +132,80 @@ class LogModel(BaseModel):
\\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m})
"""
return Utils.sdiag(np.exp(Utils.mkvc(m)))
class Vertical1DModel(BaseModel):
"""Vertical1DModel
Given a 1D vector through the last dimension
of the mesh, this will extend to the full
model space.
"""
def __init__(self, mesh, **kwargs):
BaseModel.__init__(self, mesh, **kwargs)
@property
def nP(self):
"""Number of model properties.
The number of cells in the
last dimension of the mesh."""
return self.mesh.nCv[self.mesh.dim-1]
def transform(self, m):
"""
:param numpy.array m: model
:rtype: numpy.array
:return: transformed model
"""
repNum = self.mesh.nCv[:self.mesh.dim-1].prod()
return Utils.mkvc(m).repeat(repNum)
def transformDeriv(self, m):
"""
:param numpy.array m: model
:rtype: scipy.csr_matrix
:return: derivative of transformed model
"""
repNum = self.mesh.nCv[:self.mesh.dim-1].prod()
repVec = sp.csr_matrix(
(np.ones(repNum),
(range(repNum), np.zeros(repNum))
), shape=(repNum, 1))
return sp.kron(sp.identity(self.nP), repVec)
class ComboModel(BaseModel):
"""Combination of various models."""
def __init__(self, mesh, models, **kwargs):
BaseModel.__init__(self, mesh, **kwargs)
self.models = [m(mesh, **kwargs) for m in models]
@property
def nP(self):
"""Number of model properties.
The number of cells in the
last dimension of the mesh."""
return self.models[-1].nP
def transform(self, m):
for model in reversed(self.models):
m = model.transform(m)
return m
def transformDeriv(self, m):
deriv = 1
mi = m
for model in reversed(self.models):
deriv = model.transformDeriv(mi) * deriv
mi = model.transform(mi)
return deriv
if __name__ == '__main__':
from SimPEG import *
mesh = Mesh.TensorMesh([10,8])
combo = ComboModel(mesh, [LogModel, Vertical1DModel])
m = combo.example()
print m.shape
print combo.test(np.arange(8))
SimPEG/ObjFunction.py
+4 -4
View File
@@ -3,7 +3,7 @@ import Utils, Parameters, numpy as np, scipy.sparse as sp
class BaseObjFunction(object):
"""BaseObjFunction(data, reg, **kwargs)"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
beta = Parameters.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter')
@@ -73,7 +73,7 @@ class BaseObjFunction(object):
self.u_current = None
self.m_current = m
u = self.data.prob.field(m)
u = self.data.prob.fields(m)
self.u_current = u
phi_d = self.dataObj(m, u=u)
@@ -160,7 +160,7 @@ class BaseObjFunction(object):
\\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R}
"""
if u is None: u = self.data.prob.field(m)
if u is None: u = self.data.prob.fields(m)
R = self.data.residualWeighted(m, u=u)
@@ -204,7 +204,7 @@ class BaseObjFunction(object):
\\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J}
"""
if u is None: u = self.data.prob.field(m)
if u is None: u = self.data.prob.fields(m)
R = self.data.residualWeighted(m, u=u)
SimPEG/Optimization.py
+1 -1
View File
@@ -82,7 +82,7 @@ class Minimize(object):
Minimize is a general class for derivative based optimization.
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
name = "General Optimization Algorithm" #: The name of the optimization algorithm
SimPEG/Parameters.py
+1 -1
View File
@@ -137,7 +137,7 @@ class BetaEstimate(Parameter):
u = objFunc.u_current
if u is None:
u = data.prob.field(m)
u = data.prob.fields(m)
x0 = np.random.rand(*m.shape)
t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u))
SimPEG/Problem.py
+2 -2
View File
@@ -34,7 +34,7 @@ class BaseProblem(object):
to (locally) find how model parameters change the data, and optimize!
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
@@ -142,7 +142,7 @@ class BaseProblem(object):
"""
return self.Jt(m, v, u)
def field(self, m):
def fields(self, m):
"""
The field given the model.
SimPEG/Regularization.py
+1 -1
View File
@@ -10,7 +10,7 @@ class BaseRegularization(object):
"""
__metaclass__ = Utils.Save.Savable
__metaclass__ = Utils.SimPEGMetaClass
modelPair = Model.BaseModel #: Some regularizations only work on specific models
SimPEG/Solver.py
+1 -1
View File
@@ -156,7 +156,7 @@ class Solver(object):
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
if factorize:
return self.dsolve(b)
return self.dsolve(b.flatten())
else:
return linalg.dsolve.spsolve(self.A, b)
SimPEG/Tests/test_TreeMesh.py
+503
View File
@@ -0,0 +1,503 @@
from SimPEG.Mesh import TensorMesh
from SimPEG.Mesh.TreeMesh import TreeMesh, TreeFace, TreeCell
import numpy as np
import unittest
import matplotlib.pyplot as plt
class TestOcTreeObjects(unittest.TestCase):
def setUp(self):
self.M = TreeMesh([2,1,1])
self.M.number()
self.Mr = TreeMesh([2,1,1])
self.Mr.children[0,0,0].refine()
self.Mr.number()
def q(s):
if s[0] == 'M':
m = self.M
s = s[1:]
else:
m = self.Mr
c = m.sortedCells[int(s[1])]
if len(s) == 2: return c
if s[2] == 'f' and len(s) == 5: return c.faceDict[s[2:]]
if s[2] == 'f': return getattr(c.faceDict[s[2:5]], 'edg' +s[5:])
if s[2] == 'e': return getattr(c,s[2:])
if s[2] == 'n': return getattr(c,'node'+s[3:])
self.q = q
def test_counts(self):
self.assertTrue(self.M.nC == 2)
self.assertTrue(self.M.nFx == 3)
self.assertTrue(self.M.nFy == 4)
self.assertTrue(self.M.nFz == 4)
self.assertTrue(self.M.nF == 11)
self.assertTrue(self.M.nEx == 8)
self.assertTrue(self.M.nEy == 6)
self.assertTrue(self.M.nEz == 6)
self.assertTrue(self.M.nE == 20)
self.assertTrue(self.M.nN == 12)
self.assertTrue(self.Mr.nC == 9)
self.assertTrue(self.Mr.nFx == 13)
self.assertTrue(self.Mr.nFy == 14)
self.assertTrue(self.Mr.nFz == 14)
self.assertTrue(self.Mr.nF == 41)
for cell in self.Mr.sortedCells:
for e in cell.edgeDict:
self.assertTrue(cell.edgeDict[e].edgeType==e[1].lower())
self.assertTrue(self.Mr.nN == 31)
self.assertTrue(self.Mr.nEx == 22)
self.assertTrue(self.Mr.nEy == 20)
self.assertTrue(self.Mr.nEz == 20)
def test_sizes(self):
q = self.q
for key in ['Mc0','Mc1']:
self.assertTrue(q(key).vol == 0.5)
self.assertTrue(q(key+'fXm').area == 1.)
self.assertTrue(q(key+'fXp').area == 1.)
self.assertTrue(q(key+'fYm').area == 0.5)
self.assertTrue(q(key+'fYp').area == 0.5)
self.assertTrue(q(key+'fZm').area == 0.5)
self.assertTrue(q(key+'fZp').area == 0.5)
def test_pointersM(self):
q = self.q
self.assertTrue(q('Mc0fXp') is q('Mc1fXm'))
self.assertTrue(q('Mc0fXpe0') is q('Mc1fXme0'))
self.assertTrue(q('Mc0fXpe1') is q('Mc1fXme1'))
self.assertTrue(q('Mc0fXpe2') is q('Mc1fXme2'))
self.assertTrue(q('Mc0fXpe3') is q('Mc1fXme3'))
self.assertTrue(q('Mc0fYp') is not q('c1fYm'))
self.assertTrue(q('Mc0fXm') is not q('c1fXm'))
# Test connectivity of shared edges
self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe0'))
self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe1'))
self.assertTrue(q('Mc0fZpe3') is q('Mc1fZpe2'))
self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe3'))
self.assertTrue(q('Mc0fZme3') is not q('c1fZme0'))
self.assertTrue(q('Mc0fZme3') is not q('c1fZme1'))
self.assertTrue(q('Mc0fZme3') is q('Mc1fZme2'))
self.assertTrue(q('Mc0fZme3') is not q('c1fZme3'))
self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe0'))
self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe1'))
self.assertTrue(q('Mc0fYpe3') is q('Mc1fYpe2'))
self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe3'))
self.assertTrue(q('Mc0fYme3') is not q('c1fYme0'))
self.assertTrue(q('Mc0fYme3') is not q('c1fYme1'))
self.assertTrue(q('Mc0fYme3') is q('Mc1fYme2'))
self.assertTrue(q('Mc0fYme3') is not q('c1fYme3'))
self.assertTrue(q('Mc0fZme3') is q('Mc1fXme0'))
self.assertTrue(q('Mc0fZpe3') is q('Mc1fXme1'))
self.assertTrue(q('Mc0fYme3') is q('Mc1fXme2'))
self.assertTrue(q('Mc0fYpe3') is q('Mc1fXme3'))
self.assertTrue(q('Mc0fZme3') is q('Mc0fXpe0'))
self.assertTrue(q('Mc0fZpe3') is q('Mc0fXpe1'))
self.assertTrue(q('Mc0fYme3') is q('Mc0fXpe2'))
self.assertTrue(q('Mc0fYpe3') is q('Mc0fXpe3'))
self.assertTrue(q('Mc1fZme2') is q('Mc1fXme0'))
self.assertTrue(q('Mc1fZpe2') is q('Mc1fXme1'))
self.assertTrue(q('Mc1fYme2') is q('Mc1fXme2'))
self.assertTrue(q('Mc1fYpe2') is q('Mc1fXme3'))
self.assertTrue(q('Mc1fZme2') is q('Mc0fXpe0'))
self.assertTrue(q('Mc1fZpe2') is q('Mc0fXpe1'))
self.assertTrue(q('Mc1fYme2') is q('Mc0fXpe2'))
self.assertTrue(q('Mc1fYpe2') is q('Mc0fXpe3'))
def test_nodePointers(self):
q = self.q
c0 = self.Mr.sortedCells[0]
c0n0 = c0.node0
self.assertTrue(c0n0 is q('c0n0'))
self.assertTrue(np.all(q('c0n0').center == np.r_[0,0,0.]))
self.assertTrue(q('c0n0').num == 0)
self.assertTrue(q('c0n1').num == 1)
self.assertTrue(q('c0n2').num == 4)
self.assertTrue(q('c0n3').num == 5)
self.assertTrue(q('c0n4').num == 11)
self.assertTrue(q('c0n5').num == 12)
self.assertTrue(q('c0n6').num == 14)
self.assertTrue(q('c0n7').num == 15)
def test_pointersMr(self):
q = self.q
c0 = self.Mr.sortedCells[0]
c0fXm = c0.fXm
c0eX0 = c0.eX0
c0fYme0 = c0.fYm.edge0
self.assertTrue(c0 is q('c0'))
self.assertTrue(c0fXm is q('c0fXm'))
self.assertTrue(c0eX0 is q('c0eX0'))
self.assertTrue(c0fYme0 is q('c0fYme0'))
self.assertTrue(q('c0').depth == 1)
self.assertTrue(q('c1').depth == 1)
self.assertTrue(q('c2').depth == 0)
# Make sure we know where the center of the cells are.
self.assertTrue(np.all(q('c0').center == np.r_[0.125,0.25,0.25]))
self.assertTrue(np.all(q('c1').center == np.r_[0.375,0.25,0.25]))
self.assertTrue(np.all(q('c2').center == np.r_[0.75,0.5,0.5]))
self.assertTrue(np.all(q('c3').center == np.r_[0.125,0.75,0.25]))
self.assertTrue(np.all(q('c4').center == np.r_[0.375,0.75,0.25]))
self.assertTrue(np.all(q('c5').center == np.r_[0.125,0.25,0.75]))
self.assertTrue(np.all(q('c6').center == np.r_[0.375,0.25,0.75]))
self.assertTrue(np.all(q('c7').center == np.r_[0.125,0.75,0.75]))
self.assertTrue(np.all(q('c8').center == np.r_[0.375,0.75,0.75]))
# Test X face connectivity and locations and stuff...
self.assertTrue(np.all(q('c0fXm').center == np.r_[0,0.25,0.25]))
self.assertTrue(np.all(q('c0fXp').center == np.r_[0.25,0.25,0.25]))
self.assertTrue(q('c0fXp') is q('c1fXm'))
self.assertTrue(np.all(q('c1fXp').center == np.r_[0.5,0.25,0.25]))
self.assertTrue(np.all(q('c2fXm').center == np.r_[0.5,0.5,0.5]))
self.assertTrue(q('c2fXm').branchdepth == 1)
self.assertTrue(q('c2fXm').children[0,0] is q('c1fXp'))
self.assertTrue(np.all(q('c3fXm').center == np.r_[0,0.75,0.25]))
self.assertTrue(np.all(q('c3fXp').center == np.r_[0.25,0.75,0.25]))
self.assertTrue(q('c4fXm') is q('c3fXp'))
self.assertTrue(q('c2fXm').children[1,0] is q('c4fXp'))
#Test some internal stuff (edges held by cell should be same as inside)
for key in ['Mc0', 'Mc1'] + ['c%d'%i for i in range(9)]:
self.assertTrue(q(key+'eX0') is q(key+'fZme0'))
self.assertTrue(q(key+'eX1') is q(key+'fZme1'))
self.assertTrue(q(key+'eX2') is q(key+'fZpe0'))
self.assertTrue(q(key+'eX3') is q(key+'fZpe1'))
self.assertTrue(q(key+'eX0') is q(key+'fYme0'))
self.assertTrue(q(key+'eX1') is q(key+'fYpe0'))
self.assertTrue(q(key+'eX2') is q(key+'fYme1'))
self.assertTrue(q(key+'eX3') is q(key+'fYpe1'))
self.assertTrue(q(key+'eY0') is q(key+'fXme0'))
self.assertTrue(q(key+'eY1') is q(key+'fXpe0'))
self.assertTrue(q(key+'eY2') is q(key+'fXme1'))
self.assertTrue(q(key+'eY3') is q(key+'fXpe1'))
self.assertTrue(q(key+'eY0') is q(key+'fZme2'))
self.assertTrue(q(key+'eY1') is q(key+'fZme3'))
self.assertTrue(q(key+'eY2') is q(key+'fZpe2'))
self.assertTrue(q(key+'eY3') is q(key+'fZpe3'))
self.assertTrue(q(key+'eZ0') is q(key+'fXme2'))
self.assertTrue(q(key+'eZ1') is q(key+'fXpe2'))
self.assertTrue(q(key+'eZ2') is q(key+'fXme3'))
self.assertTrue(q(key+'eZ3') is q(key+'fXpe3'))
self.assertTrue(q(key+'eZ0') is q(key+'fYme2'))
self.assertTrue(q(key+'eZ1') is q(key+'fYme3'))
self.assertTrue(q(key+'eZ2') is q(key+'fYpe2'))
self.assertTrue(q(key+'eZ3') is q(key+'fYpe3'))
#Test some edge stuff
self.assertTrue(np.all(q('c0eX0').center == np.r_[0.125,0,0]))
self.assertTrue(np.all(q('c0eX1').center == np.r_[0.125,0.5,0]))
self.assertTrue(np.all(q('c0eX2').center == np.r_[0.125,0,0.5]))
self.assertTrue(np.all(q('c0eX3').center == np.r_[0.125,0.5,0.5]))
self.assertTrue(np.all(q('c5eX0').center == np.r_[0.125,0,0.5]))
self.assertTrue(np.all(q('c5eX1').center == np.r_[0.125,0.5,0.5]))
self.assertTrue(q('c5eX0') is q('c0eX2'))
self.assertTrue(q('c5eX1') is q('c0eX3'))
self.assertTrue(np.all(q('c0eY0').center == np.r_[0,0.25,0]))
self.assertTrue(np.all(q('c0eY1').center == np.r_[0.25,0.25,0]))
self.assertTrue(np.all(q('c0eY2').center == np.r_[0,0.25,0.5]))
self.assertTrue(np.all(q('c0eY3').center == np.r_[0.25,0.25,0.5]))
self.assertTrue(np.all(q('c1eY0').center == np.r_[0.25,0.25,0]))
self.assertTrue(np.all(q('c1eY2').center == np.r_[0.25,0.25,0.5]))
self.assertTrue(q('c1eY0') is q('c0eY1'))
self.assertTrue(q('c1eY2') is q('c0eY3'))
self.assertTrue(np.all(q('c0eZ0').center == np.r_[0,0,0.25]))
self.assertTrue(np.all(q('c0eZ1').center == np.r_[0.25,0,0.25]))
self.assertTrue(np.all(q('c0eZ2').center == np.r_[0,0.5,0.25]))
self.assertTrue(np.all(q('c0eZ3').center == np.r_[0.25,0.5,0.25]))
self.assertTrue(np.all(q('c3eZ0').center == np.r_[0,0.5,0.25]))
self.assertTrue(np.all(q('c3eZ1').center == np.r_[0.25,0.5,0.25]))
self.assertTrue(q('c3eZ0') is q('c0eZ2'))
self.assertTrue(q('c3eZ1') is q('c0eZ3'))
self.assertTrue(q('c0fXp') is q('c1fXm'))
self.assertTrue(q('c0fYp') is not q('c1fYm'))
self.assertTrue(q('c0fXm') is not q('c1fXm'))
self.assertTrue(q('c1fXp') is q('c2fXm').children[0,0])
self.assertTrue(q('c1fYp') is q('c4fYm'))
self.assertTrue(q('c1fZp') is q('c6fZm'))
self.assertTrue(q('c6fXp') is q('c2fXm').children[0,1])
self.assertTrue(q('c4fXp') is q('c2fXm').children[1,0])
def test_gridCC(self):
x = np.r_[0.25,0.75]
y = np.r_[0.5,0.5]
z = np.r_[0.5,0.5]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridCC).flatten()) == 0)
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375]
y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75]
z = np.r_[0.25,0.25,0.5,0.25,0.25,0.75,0.75,0.75,0.75]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridCC).flatten()) == 0)
def test_gridN(self):
x = np.r_[0,0.5,1,0,0.5,1,0,0.5,1,0,0.5,1]
y = np.r_[0,0,0,1,1,1,0,0,0,1,1,1.]
z = np.r_[0,0,0,0,0,0,1,1,1,1,1,1.]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridN).flatten()) == 0)
x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1]
y = np.r_[0,0,0,0,0.5,0.5,0.5,1,1,1,1,0,0,0,0.5,0.5,0.5,1,1,1,0,0,0,0,0.5,0.5,0.5,1,1,1,1]
z = np.r_[0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1,1,1,1]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridN).flatten()) == 0)
def test_gridFx(self):
x = np.r_[0.0,0.5,1.0]
y = np.r_[0.5,0.5,0.5]
z = np.r_[0.5,0.5,0.5]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFx).flatten()) == 0)
x = np.r_[0.0,0.25,0.5,1.0,0.0,0.25,0.5,0.0,0.25,0.5,0.0,0.25,0.5]
y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75]
z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.75,0.75,0.75,0.75,0.75,0.75]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFx).flatten()) == 0)
def test_gridFy(self):
x = np.r_[0.25,0.75,0.25,0.75]
y = np.r_[0,0,1.,1.]
z = np.r_[0.5,0.5,0.5,0.5]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFy).flatten()) == 0)
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375]
y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1]
z = np.r_[0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFy).flatten()) == 0)
def test_gridFz(self):
x = np.r_[0.25,0.75,0.25,0.75]
y = np.r_[0.5,0.5,0.5,0.5]
z = np.r_[0,0,1.,1.]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFz).flatten()) == 0)
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375]
y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75,0.25,0.25,0.5,0.75,0.75]
z = np.r_[0,0,0,0,0,0.5,0.5,0.5,0.5,1,1,1,1,1]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFz).flatten()) == 0)
def test_gridEx(self):
x = np.r_[0.25,0.75,0.25,0.75,0.25,0.75,0.25,0.75]
y = np.r_[0,0,1.,1.,0,0,1.,1.]
z = np.r_[0,0,0,0,1.,1.,1.,1.]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEx).flatten()) == 0)
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75]
y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1,0,0,0,0.5,0.5,1,1,1]
z = np.r_[0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEx).flatten()) == 0)
def test_gridEy(self):
x = np.r_[0,0.5,1,0,0.5,1]
y = np.r_[0.5,0.5,0.5,0.5,0.5,0.5]
z = np.r_[0,0,0,1.,1.,1.]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEy).flatten()) == 0)
x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5]
y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75,0.25,0.25,0.25,0.5,0.75,0.75,0.75]
z = np.r_[0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEy).flatten()) == 0)
def test_gridEz(self):
x = np.r_[0,0.5,1,0,0.5,1]
y = np.r_[0,0,0,1.,1.,1.]
z = np.r_[0.5,0.5,0.5,0.5,0.5,0.5]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEz).flatten()) == 0)
x = np.r_[0,0.25,0.5,1,0 ,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0 ,0.25,0.5,0 ,0.25,0.5]
y = np.r_[0,0 ,0 ,0,0.5,0.5 ,0.5,1,1 ,1 ,1,0,0 ,0 ,0.5,0.5 ,0.5,1 ,1 ,1 ]
z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75]
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEz).flatten()) == 0)
class TestQuadTreeObjects(unittest.TestCase):
def setUp(self):
self.M = TreeMesh([2,1])
self.Mr = TreeMesh([2,1])
self.Mr.children[0,0].refine()
self.Mr.number()
# self.Mr.plotGrid(showIt=True)
def test_pointersM(self):
c0 = self.M.children[0,0]
c0fXm = c0.fXm
c0fXp = c0.fXp
c0fYm = c0.fYm
c0fYp = c0.fYp
c1 = self.M.children[1,0]
c1fXm = c1.fXm
c1fXp = c1.fXp
c1fYm = c1.fYm
c1fYp = c1.fYp
self.assertTrue(c0fXp is c1fXm)
self.assertTrue(c0fYp is not c1fYm)
self.assertTrue(c0fXm is not c1fXm)
self.assertTrue(c0fXm.area == 1)
self.assertTrue(c0fYm.area == 0.5)
self.assertTrue(c0.node1 is c1.node0)
self.assertTrue(c0.node3 is c1.node2)
self.assertTrue(self.M.nN == 6)
def test_pointersMr(self):
c0 = self.Mr.sortedCells[0]
c0fXm = c0.fXm
c0fXp = c0.fXp
c0fYm = c0.fYm
c0fYp = c0.fYp
c1 = self.Mr.sortedCells[1]
c1fXm = c1.fXm
c1fXp = c1.fXp
c1fYm = c1.fYm
c1fYp = c1.fYp
c2 = self.Mr.sortedCells[2]
c2fXm = c2.fXm
c2fXp = c2.fXp
c2fYm = c2.fYm
c2fYp = c2.fYp
c4 = self.Mr.sortedCells[4]
c4fXm = c4.fXm
c4fXp = c4.fXp
c4fYm = c4.fYm
c4fYp = c4.fYp
self.assertTrue(c0fXp is c1fXm)
self.assertTrue(c1fXp.node0 is c2fXm.node0)
self.assertTrue(c1fXp.node0 is c2fXm.node0)
self.assertTrue(c4fYm is c1fYp)
self.assertTrue(c4fXp.node1 is c2fXm.node1)
self.assertTrue(c4fXp.node0 is c1fYp.node1)
self.assertTrue(c0fXp.node1 is c4fYm.node0)
self.assertTrue(self.Mr.nN == 11)
self.assertTrue(np.all(c1fXp.node0.x0 == np.r_[0.5,0]))
self.assertTrue(np.all(c1fYp.node0.x0 == np.r_[0.25,0.5]))
class TestQuadTreeMesh(unittest.TestCase):
def setUp(self):
M = TreeMesh([np.ones(x) for x in [3,2]])
for ii in range(1):
M.children[ii,ii].refine()
self.M = M
M.number()
# M.plotGrid(showIt=True)
def test_MeshSizes(self):
self.assertTrue(self.M.nC==9)
self.assertTrue(self.M.nF==25)
self.assertTrue(self.M.nFx==12)
self.assertTrue(self.M.nFy==13)
self.assertTrue(self.M.nE==25)
self.assertTrue(self.M.nEx==13)
self.assertTrue(self.M.nEy==12)
def test_gridCC(self):
x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.5,1.5,2.5]
y = np.r_[0.25,0.25,0.5,0.5,0.75,0.75,1.5,1.5,1.5]
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridCC).flatten()) == 0)
def test_gridN(self):
x = np.r_[0,0.5,1,2,3,0,0.5,1,0,0.5,1,2,3,0,1,2,3]
y = np.r_[0,0,0,0,0,.5,.5,.5,1,1,1,1,1,2,2,2,2]
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridN).flatten()) == 0)
def test_gridFx(self):
x = np.r_[0.0,0.5,1.0,2.0,3.0,0.0,0.5,1.0,0.0,1.0,2.0,3.0]
y = np.r_[0.25,0.25,0.25,0.5,0.5,0.75,0.75,0.75,1.5,1.5,1.5,1.5]
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridFx).flatten()) == 0)
def test_gridFy(self):
x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.25,0.75,1.5,2.5,0.5,1.5,2.5]
y = np.r_[0,0,0,0,0.5,0.5,1,1,1,1,2,2,2]
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridFy).flatten()) == 0)
def test_gridEx(self):
x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.25,0.75,1.5,2.5,0.5,1.5,2.5]
y = np.r_[0,0,0,0,0.5,0.5,1,1,1,1,2,2,2]
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridEx).flatten()) == 0)
def test_gridEy(self):
x = np.r_[0.0,0.5,1.0,2.0,3.0,0.0,0.5,1.0,0.0,1.0,2.0,3.0]
y = np.r_[0.25,0.25,0.25,0.5,0.5,0.75,0.75,0.75,1.5,1.5,1.5,1.5]
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridEy).flatten()) == 0)
class SimpleOctreeOperatorTests(unittest.TestCase):
def setUp(self):
h1 = np.random.rand(5)
h2 = np.random.rand(7)
h3 = np.random.rand(3)
self.tM = TensorMesh([h1,h2,h3])
self.oM = TreeMesh([h1,h2,h3])
self.tM2 = TensorMesh([h1,h2])
self.oM2 = TreeMesh([h1,h2])
def test_faceDiv(self):
self.assertTrue((self.tM.faceDiv - self.oM.faceDiv).toarray().sum() == 0)
self.assertTrue((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum() == 0)
def test_nodalGrad(self):
self.assertTrue((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum() == 0)
self.assertTrue((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum() == 0)
def test_edgeCurl(self):
self.assertTrue((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum() == 0)
# self.assertTrue((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum() == 0)
def test_InnerProducts(self):
self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() == 0)
self.assertTrue((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum() == 0)
self.assertTrue((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum() == 0)
self.assertTrue((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum() == 0)
if __name__ == '__main__':
unittest.main()
SimPEG/Tests/test_forward_DCproblem.py
-85
View File
@@ -1,85 +0,0 @@
# import numpy as np
# import unittest
# from SimPEG.mesh import TensorMesh
# from SimPEG.Utils import ModelBuilder, sdiag
# from SimPEG.forward import Problem
# from SimPEG.examples.DC import *
# from TestUtils import checkDerivative
# from scipy.sparse.linalg import dsolve
# from SimPEG import inverse
# class DCProblemTests(unittest.TestCase):
# def setUp(self):
# # Create the mesh
# h1 = np.ones(20)
# h2 = np.ones(20)
# mesh = TensorMesh([h1,h2])
# # Create some parameters for the model
# sig1 = 1
# sig2 = 0.01
# # Create a synthetic model from a block in a half-space
# p0 = [2, 2]
# p1 = [5, 5]
# condVals = [sig1, sig2]
# mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
# # Set up the projection
# nelec = 10
# spacelec = 2
# surfloc = 0.5
# elecini = 0.5
# elecend = 0.5+spacelec*(nelec-1)
# elecLocR = np.linspace(elecini, elecend, nelec)
# rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
# q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
# P = Q.T
# # Create some data
# problem = DCProblem(mesh)
# problem.P = P
# problem.RHS = q
# data = problem.createSyntheticData(mSynth, std=0.05)
# # Now set up the problem to do some minimization
# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
# reg = inverse.Regularization(mesh)
# inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4)
# self.inv = inv
# self.reg = reg
# self.p = problem
# self.mesh = mesh
# self.m0 = mSynth
# self.data = data
# def test_misfit(self):
# derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)]
# passed = checkDerivative(derChk, self.m0, plotIt=False)
# self.assertTrue(passed)
# def test_adjoint(self):
# # Adjoint Test
# u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1])
# v = np.random.rand(self.mesh.nC)
# w = np.random.rand(self.data.dobs.shape[0])
# wtJv = w.dot(self.p.J(self.m0, v, u=u))
# vtJtw = v.dot(self.p.Jt(self.m0, w, u=u))
# passed = (wtJv - vtJtw) < 1e-10
# self.assertTrue(passed)
# def test_dataObj(self):
# derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)]
# checkDerivative(derChk, self.m0, plotIt=False)
# def test_modelObj(self):
# derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)]
# checkDerivative(derChk, self.m0, plotIt=False)
# if __name__ == '__main__':
# unittest.main()
SimPEG/Tests/test_interpolation.py
+2 -2
View File
@@ -4,7 +4,7 @@ from TestUtils import OrderTest
from SimPEG.Utils import mkvc
MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
TOLERANCES = [0.9, 0.55]
TOLERANCES = [0.9, 0.5]
call1 = lambda fun, xyz: fun(xyz)
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
@@ -23,7 +23,7 @@ class TestInterpolation1D(OrderTest):
meshTypes = MESHTYPES
tolerance = TOLERANCES
meshDimension = 1
meshSizes = [8, 16, 32]
meshSizes = [8, 16, 32, 64, 128]
def getError(self):
funX = lambda x: np.cos(2*np.pi*x)
SimPEG/Tests/test_model.py
+13 -7
View File
@@ -11,17 +11,23 @@ class ModelTests(unittest.TestCase):
a = np.array([1, 1, 1])
b = np.array([1, 2])
c = np.array([1, 4])
self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5]))
def test_modelTransforms(self):
print 'SimPEG.Model.BaseModel: Testing Model Transform'
for M in dir(Model):
if 'Model' not in M: continue
model = getattr(Model, M)(self.mesh2)
m = model.example()
passed = checkDerivative(lambda m : [model.transform(m), model.transformDeriv(m)], m, plotIt=False)
self.assertTrue(passed)
try:
model = getattr(Model, M)(self.mesh2)
assert isinstance(model, Model.BaseModel)
except Exception, e:
continue
self.assertTrue(model.test())
def test_comboModels(self):
combos = [(Model.LogModel, Model.Vertical1DModel)]
for combo in combos:
model = Model.ComboModel(self.mesh2, combo)
self.assertTrue(model.test())
if __name__ == '__main__':
unittest.main()
SimPEG/Tests/test_operators.py
+37 -1
View File
@@ -338,6 +338,15 @@ class TestAveraging2D(OrderTest):
self.getAve = lambda M: M.aveF2CC
self.orderTest()
def test_orderF2CCV(self):
self.name = "Averaging 2D: F2CCV"
funX = lambda x, y: (np.cos(x)+np.sin(y))
funY = lambda x, y: (np.cos(y)*np.sin(x))
self.getHere = lambda M: np.r_[call2(funX, M.gridFx), call2(funY, M.gridFy)]
self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)]
self.getAve = lambda M: M.aveF2CCV
self.orderTest()
def test_orderCC2F(self):
self.name = "Averaging 2D: CC2F"
fun = lambda x, y: (np.cos(x)+np.sin(y))
@@ -348,7 +357,6 @@ class TestAveraging2D(OrderTest):
self.orderTest()
self.expectedOrders = 2
def test_orderE2CC(self):
self.name = "Averaging 2D: E2CC"
fun = lambda x, y: (np.cos(x)+np.sin(y))
@@ -357,6 +365,15 @@ class TestAveraging2D(OrderTest):
self.getAve = lambda M: M.aveE2CC
self.orderTest()
def test_orderE2CCV(self):
self.name = "Averaging 2D: E2CCV"
funX = lambda x, y: (np.cos(x)+np.sin(y))
funY = lambda x, y: (np.cos(y)*np.sin(x))
self.getHere = lambda M: np.r_[call2(funX, M.gridEx), call2(funY, M.gridEy)]
self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)]
self.getAve = lambda M: M.aveE2CCV
self.orderTest()
class TestAveraging3D(OrderTest):
name = "Averaging 3D"
@@ -400,6 +417,15 @@ class TestAveraging3D(OrderTest):
self.getAve = lambda M: M.aveF2CC
self.orderTest()
def test_orderF2CCV(self):
self.name = "Averaging 3D: F2CCV"
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
funZ = lambda x, y, z: (np.cos(x)*np.sin(y)+np.exp(z))
self.getHere = lambda M: np.r_[call3(funX, M.gridFx), call3(funY, M.gridFy), call3(funZ, M.gridFz)]
self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)]
self.getAve = lambda M: M.aveF2CCV
self.orderTest()
def test_orderE2CC(self):
self.name = "Averaging 3D: E2CC"
@@ -409,6 +435,16 @@ class TestAveraging3D(OrderTest):
self.getAve = lambda M: M.aveE2CC
self.orderTest()
def test_orderE2CCV(self):
self.name = "Averaging 3D: E2CCV"
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
funZ = lambda x, y, z: (np.cos(x)*np.sin(y)+np.exp(z))
self.getHere = lambda M: np.r_[call3(funX, M.gridEx), call3(funY, M.gridEy), call3(funZ, M.gridEz)]
self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)]
self.getAve = lambda M: M.aveE2CCV
self.orderTest()
def test_orderCC2F(self):
self.name = "Averaging 3D: CC2F"
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
SimPEG/Tests/test_tensorMesh.py
+3 -3
View File
@@ -2,7 +2,7 @@ import numpy as np
import unittest
from SimPEG.Mesh import TensorMesh
from TestUtils import OrderTest
from scipy.sparse.linalg import dsolve
from SimPEG import Solver
class BasicTensorMeshTests(unittest.TestCase):
@@ -58,7 +58,7 @@ class BasicTensorMeshTests(unittest.TestCase):
class TestPoissonEqn(OrderTest):
name = "Poisson Equation"
meshSizes = [16, 20, 24]
meshSizes = [10, 16, 20]
def getError(self):
# Create some functions to integrate
@@ -75,7 +75,7 @@ class TestPoissonEqn(OrderTest):
err = np.linalg.norm((sA - sN), np.inf)
else:
fA = fun(self.M.gridCC)
fN = dsolve.spsolve(D*G, sol(self.M.gridCC))
fN = Solver(D*G).solve(sol(self.M.gridCC))
err = np.linalg.norm((fA - fN), np.inf)
return err
SimPEG/Utils/Geophysics/__init__.py
-1
View File
@@ -1 +0,0 @@
import emSources
SimPEG/Utils/Geophysics/emSources/__init__.py
-1
View File
@@ -1 +0,0 @@
from emSources import MagneticDipoleVectorPotential
SimPEG/Utils/Geophysics/emSources/emSources.py
-40
View File
@@ -1,40 +0,0 @@
import numpy as np
from scipy.constants import mu_0, pi
def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
"""
Calculate the vector potential of a set of magnetic dipoles
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
:param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z)
:param str component: The component to calculate - 'x', 'y', or 'z'
:param numpy.ndarray dipoleMoment: The vector dipole moment
:rtype: numpy.ndarray
:return: The vector potential each dipole at each observation location
"""
if component=='x':
dimInd = 0
elif component=='y':
dimInd = 1
elif component=='z':
dimInd = 2
else:
raise ValueError('Invalid component')
txLoc = np.atleast_2d(txLoc)
obsLoc = np.atleast_2d(obsLoc)
dipoleMoment = np.atleast_2d(dipoleMoment)
nEdges = obsLoc.shape[0]
nTx = txLoc.shape[0]
m = np.array(dipoleMoment).repeat(nEdges, axis=0)
A = np.empty((nEdges, nTx))
for i in range(nTx):
dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0)
mCr = np.cross(m, dR)
r = np.sqrt((dR**2).sum(axis=1))
A[:, i] = -(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
return A
SimPEG/Utils/Save.py
-352
View File
@@ -1,352 +0,0 @@
import numpy as np
import time
import re
try:
import h5py
except Exception, e:
print 'Warning: SimPEG.Utils.Save needs h5py to be installed.'
SAVEABLES = {}
def natural_keys(text):
'''
alist.sort(key=natural_keys) sorts in human order
http://nedbatchelder.com/blog/200712/human_sorting.html
(See Toothy's implementation in the comments)
'''
atoi = lambda text: int(text) if text.isdigit() else text
return [ atoi(c) for c in re.split('(\d+)', text) ]
def preIteration(group):
group.attrs['complete'] = False
group.attrs['time'] = time.time()
def postIteration(group):
group.attrs['time'] = time.time() - group.attrs['time']
group.attrs['date'] = time.ctime()
group.attrs['complete'] = True
class SimPEGTable:
"""
This is a wrapper class on the HDF5 file.
"""
def __init__(self, name, mode='a'):
if '.hdf5' not in name:
name += '.hdf5'
self.f = h5py.File(name, mode)
self.root = hdf5Group(self,self.f)
self.inversions = hdf5InversionGroup(self,self.root.addGroup('inversions',soft=True))
def show(self): self.root.show()
def saveInversion(self, invObj):
# Create a new inversion anytime this is run.
def _startup_hdf5_inv(invObj, m0):
node = self.inversions.addGroup('%d'%self.inversions.numChildren)
saveSavable(invObj,node.addGroup('rebuild'))
results = node.addGroup('results')
preIteration(results)
invObj._invNode = results
self.f.flush()
invObj.hook(_startup_hdf5_inv, overwrite=True)
# At the start of every iteration we will create a inversion iteration node.
def _doStartIteration_hdf5_inv(invObj):
invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj.iter+1))
preIteration(invObj._invNodeIt)
invObj.hook(_doStartIteration_hdf5_inv, overwrite=True)
def _doEndIteration_hdf5_inv(invObj):
invObj.save(invObj._invNodeIt)
postIteration(invObj._invNodeIt)
self.f.flush()
invObj.hook(_doEndIteration_hdf5_inv, overwrite=True)
# Delete all iterates that did not finish.
def _finish_hdf5_inv(invObj):
postIteration(invObj._invNode)
for it in invObj._invNode:
if not it.attrs['complete']:
del self.f[it.path]
del invObj._invNode
self.f.flush()
invObj.hook(_finish_hdf5_inv, overwrite=True)
def _doStartIteration_hdf5_opt(optObj):
optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent.iter, optObj.iter))
preIteration(optObj._optNodeIt)
invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True)
def _doEndIteration_hdf5_opt(optObj, xt):
optObj.save(optObj._optNodeIt)
postIteration(optObj._optNodeIt)
self.f.flush()
invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True)
class hdf5Group(object):
"""Has some low level support for wrapping the native HDF5-Group class"""
def __init__(self, T, groupNode):
self.T = T
# check if you are inputing a hdf5Group rather than a raw node, and act accordingly
if issubclass(groupNode.__class__, hdf5Group):
self.node = groupNode.node
else:
self.node = groupNode
self.childClass = hdf5Group
self.parentClass = hdf5Group
@property
def children(self):
"""Children names in a list
Use obj[name] to return the actual node.
"""
myChildren = [c for c in self.node]
myChildren.sort(key=natural_keys)
return myChildren
@property
def numChildren(self):
"""Returns the len(children)"""
return len(self.children)
@property
def parent(self):
"""Returns parent node"""
return self.parentClass(self.T, self.node.parent)
@property
def name(self):
return self.path.split('/')[-1]
@property
def path(self):
"""Returns the root path of the group"""
return self.node.name
@property
def attrs(self):
"""Returns a list of attributes in the group"""
return self.node.attrs
def addGroup(self, name, soft=False):
"""Adds a child group to the current node."""
if name in self.children and soft:
return self[name]
assert name not in self.children, 'Already a child called: '+self.path+'/'+name
return self.childClass(self.T, self.node.create_group(name))
def setArray(self, name, data):
a = self.node.create_dataset(name, data.shape)
a[...] = data
return a
def __getitem__(self, val):
if type(val) is int:
val = self.children[val]
child = self.node[val]
if type(child) is h5py.Group:
child = self.childClass(self.T, child)
return child
def __contains__(self, key):
return key in self.children
def show(self, pad='', maxDepth=1, depth=0):
"""
Recursively show the structure of the database.
For example::
<hdf5InversionGroup group "/inversions" (1 member)>
- <hdf5Inversion group "/inversions/0" (4 members)>
- <hdf5InversionIteration group "/inversions/0/0.0" (3 members)>
- <hdf5InversionIteration group "/inversions/0/0.1" (3 members)>
- <hdf5InversionIteration group "/inversions/0/0.2" (3 members)>
- <hdf5InversionIteration group "/inversions/0/0.3" (3 members)>
"""
s = self.__str__()
pad += ' '*4
if maxDepth <= 0: print s
if depth >= maxDepth: return s
for c in self.children:
if issubclass(self[c].__class__, hdf5Group):
s += '\n%s- %s' % (pad, self[c].show(pad=pad,depth=depth+1,maxDepth=maxDepth))
else:
s += '\n%s- %s' % (pad, self[c].__str__())
if depth is 0:
print s
else:
return s
def __str__(self):
return '<%s "%s" (%d member%s, attrs=[%s])>' % (self.__class__.__name__, self.path, self.numChildren, '' if self.numChildren == 1 else 's',', '.join([a for a in self.attrs]))
class hdf5InversionGroup(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.childClass = hdf5Inversion
class hdf5Inversion(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5InversionGroup
self.childClass = hdf5InversionResults
def rebuild(self):
return loadSavable(self['rebuild'])
@property
def results(self): return self['results']
class hdf5InversionResults(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5Inversion
self.childClass = hdf5InversionIteration
class hdf5InversionIteration(hdf5Group):
def __init__(self, T, groupNode):
hdf5Group.__init__(self, T, groupNode)
self.parentClass = hdf5InversionResults
class Savable(type):
def __new__(cls, name, bases, attrs):
__init__ = attrs['__init__']
def init_wrapper(self, *args, **kwargs):
self._args_init = args
self._kwargs_init = kwargs
return __init__(self, *args,**kwargs)
attrs['__init__'] = init_wrapper
newClass = super(Savable, cls).__new__(cls, name, bases, attrs)
SAVEABLES[name] = newClass
return newClass
def saveSavable(obj, group, debug=False):
"""
This creates softlinks if _savable exists in children object.
The first object is always created.
"""
assert type(obj.__class__) is Savable, 'Can only save objects that are Savable objects.'
def doSave(grp, name, val):
if debug: print name, val
if type(val.__class__) is Savable:
link = getattr(val,'_savable',None)
if link is not None:
group.node[name] = h5py.SoftLink(link.path)
if debug: 'Created a softlink path to %s' % link.path
else:
subgrp = grp.addGroup(name)
saveSavable(val, subgrp, debug=debug)
elif type(val) in [list, tuple]:
# Split up, and save each element
for i, v in enumerate(val):
doSave(grp, name + '[%d]'%i, v)
elif type(val) is np.ndarray:
grp.setArray(name, val)
elif val is None:
grp.attrs[name] = 'None'
else:
# just try saving it as an attr
try:
grp.attrs[name] = val
except Exception, e:
print 'Warning: Could not save %s, problems may arise in loading.' % name
group.attrs['__class__'] = obj.__class__.__name__
for arg in obj._kwargs_init:
doSave(group, '_kwarg_'+arg, obj._kwargs_init[arg])
for i, arg in enumerate(obj._args_init):
doSave(group, '_arg%d'%i, arg)
obj._savable = group
def loadSavable(node, pointers=None):
"""
pointers allow things that point to the same node in the h5py file to
be returned as the same object, if they have already been created.
"""
if pointers is None: pointers = []
for pointer in pointers:
if pointer._savable.node == node.node: return pointer
args = ([a for a in node.attrs if '_arg' in a] + [a for a in node.children if '_arg' in a])
kwargs = ([a for a in node.attrs if '_kwarg' in a] + [a for a in node.children if '_kwarg' in a])
args.sort(key=natural_keys)
kwargs.sort(key=natural_keys)
def get(node,key):
if key in node.children: return node[key]
elif key in node.attrs: return node.attrs[key]
ARGS = []
for name in args:
val = get(node, name)
if val.__class__ is h5py.Dataset: val = val[:]
if val is 'None': val = None
if '[' in name: # We are reloading a list
ind = int(name[4:name.index('[')])
if len(ARGS) is ind: # Create the list
ARGS.append([val])
else:
ARGS[ind].append(val)
elif issubclass(val.__class__,hdf5Group):
ARGS.append(loadSavable(val,pointers=pointers))
else:
ind = int(name[4:])
ARGS.append(val)
KWARGS = {}
for name in kwargs:
val = get(node, name)
if val.__class__ is h5py.Dataset: val = val[:]
if val is 'None': val = None
if '[' in name: # We are reloading a list
key = name[7:name.index('[')]
if key not in KWARGS: # Create the list
KWARGS[key] = [val]
else:
KWARGS[key].append(val)
elif issubclass(val.__class__,hdf5Group):
key = name[7:]
KWARGS[key] = loadSavable(val,pointers=pointers)
else:
key = name[7:]
KWARGS[key] = val
cls = get(node, '__class__')
if cls in SAVEABLES:
try:
out = SAVEABLES[cls](*ARGS, **KWARGS)
out._savable = node
pointers.append(out) # Because this is recursive.
return out
except Exception, e:
print 'Warning: %s Class could not be initiated.' % cls
print 'ARGS: ', ARGS
print 'KWARGS: ', KWARGS
return (cls, ARGS, KWARGS, node)
else:
print 'Warning: %s Class not found in SimPEG.Utils.Save.SAVABLES' % cls
return (cls, ARGS, KWARGS, node)
SimPEG/Utils/__init__.py
+8 -3
View File
@@ -1,11 +1,9 @@
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap
from sputils import spzeros, kron3, speye, sdiag, sdInv, ddx, av, avExtrap
from meshutils import exampleLomGird, meshTensors
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube
from interputils import interpmat
from ipythonutils import easyAnimate as animate
import Save
import Geophysics
import ModelBuilder
import types
@@ -13,6 +11,12 @@ import time
import numpy as np
from functools import wraps
class SimPEGMetaClass(type):
def __new__(cls, name, bases, attrs):
return super(SimPEGMetaClass, cls).__new__(cls, name, bases, attrs)
def hook(obj, method, name=None, overwrite=False, silent=False):
"""
This dynamically binds a method to the instance of the class.
@@ -132,6 +136,7 @@ def callHooks(match, mainFirst=False):
def dependentProperty(name, value, children, doc):
def fget(self): return getattr(self,name,value)
def fset(self, val):
if getattr(self,name,value) == val: return # it is the same!
for child in children:
if hasattr(self, child):
delattr(self, child)
SimPEG/Utils/sputils.py
+3
View File
@@ -7,6 +7,9 @@ def sdiag(h):
"""Sparse diagonal matrix"""
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
def sdInv(M):
"Inverse of a sparse diagonal matrix"
return sdiag(1/M.diagonal())
def speye(n):
"""Sparse identity"""
SimPEG/__init__.py
-1
View File
@@ -11,7 +11,6 @@ import ObjFunction
import Optimization
import Inversion
import Parameters
import Examples
import Tests
SimPEG/visualize/__init__.py
-2
View File
@@ -1,2 +0,0 @@
import vtk
#import mpl
SimPEG/visualize/vtk/__init__.py
-2
View File
@@ -1,2 +0,0 @@
from vtkTools import vtkTools
from vtkView import vtkView
SimPEG/visualize/vtk/vtkTools.py
-385
View File
@@ -1,385 +0,0 @@
import numpy as np
try:
import vtk, vtk.util.numpy_support as npsup, pdb
except Exception, e:
print 'VTK import error. Please ensure you have VTK installed to use this visualization package.'
from SimPEG.Utils import mkvc
class vtkTools(object):
"""
Class that interacts with VTK visulization toolkit.
"""
def __init__(self):
""" Initializes the VTK vtkTools.
"""
pass
@staticmethod
def makeCellVTKObject(mesh,model):
"""
Make and return a cell based VTK object for a simpeg mesh and model.
Input:
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells
Output:
:rtype: vtkRecilinearGrid object
:return: vtkObj
"""
# Deal with dimensionalities
if mesh.dim >= 1:
vX = mesh.vectorNx
xD = mesh.nNx
yD,zD = 1,1
vY, vZ = np.array([0,0])
if mesh.dim >= 2:
vY = mesh.vectorNy
yD = mesh.nNy
if mesh.dim == 3:
vZ = mesh.vectorNz
zD = mesh.nNz
# Use rectilinear VTK grid.
# Assign the spatial information.
vtkObj = vtk.vtkRectilinearGrid()
vtkObj.SetDimensions(xD,yD,zD)
vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1))
vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1))
vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1))
# Assign the model('s) to the object
for item in model.iteritems():
# Convert numpy array
vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
vtkObj.Update()
return vtkObj
@staticmethod
def makeFaceVTKObject(mesh,model):
"""
Make and return a face based VTK object for a simpeg mesh and model.
Input:
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param model, dictionary of numpy.array - Name('s) and array('s).
Property array must be order hstack(Fx,Fy,Fz)
Output:
:rtype: vtkUnstructuredGrid object
:return: vtkObj
"""
## Convert simpeg mesh to VTK properties
# Convert mesh nodes to vtkPoints
vtkPts = vtk.vtkPoints()
vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1))
# Define the face "cells"
# Using VTK_QUAD cell for faces (see VTK file format)
nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M')
def faceR(mat,length):
return mat.T.reshape((length,1))
# First direction
nTFx = np.prod(mesh.nFx)
FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] )
FyCellBlock = np.array([],dtype='int64')
FzCellBlock = np.array([],dtype='int64')
# Second direction
if mesh.dim >= 2:
nTFy = np.prod(mesh.nFy)
FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] )
# Third direction
if mesh.dim == 3:
nTFz = np.prod(mesh.nFz)
FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] )
# Cells -cell array
FCellArr = vtk.vtkCellArray()
FCellArr.SetNumberOfCells(mesh.nF)
FCellArr.SetCells(mesh.nF,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1))
# Cell type
FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1)
# Cell location
FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nF*5,5,dtype='int64'),deep=1)
## Make the object
vtkObj = vtk.vtkUnstructuredGrid()
# Set the objects properties
vtkObj.SetPoints(vtkPts)
vtkObj.SetCells(FCellType,FCellLoc,FCellArr)
# Assign the model('s) to the object
for item in model.iteritems():
# Convert numpy array
vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
vtkObj.Update()
return vtkObj
@staticmethod
def makeEdgeVTKObject(mesh,model):
"""
Make and return a edge based VTK object for a simpeg mesh and model.
Input:
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param model, dictionary of numpy.array - Name('s) and array('s).
Property array must be order hstack(Ex,Ey,Ez)
Output:
:rtype: vtkUnstructuredGrid object
:return: vtkObj
"""
## Convert simpeg mesh to VTK properties
# Convert mesh nodes to vtkPoints
vtkPts = vtk.vtkPoints()
vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1))
# Define the face "cells"
# Using VTK_QUAD cell for faces (see VTK file format)
nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M')
def edgeR(mat,length):
return mat.T.reshape((length,1))
# First direction
nTEx = np.prod(mesh.nEx)
ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)])
# Second direction
if mesh.dim >= 2:
nTEy = np.prod(mesh.nEy)
EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)])
# Third direction
if mesh.dim == 3:
nTEz = np.prod(mesh.nEz)
EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)])
# Cells -cell array
ECellArr = vtk.vtkCellArray()
ECellArr.SetNumberOfCells(mesh.nE)
ECellArr.SetCells(mesh.nE,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1))
# Cell type
ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1)
# Cell location
ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nE*3,3,dtype='int64'),deep=1)
## Make the object
vtkObj = vtk.vtkUnstructuredGrid()
# Set the objects properties
vtkObj.SetPoints(vtkPts)
vtkObj.SetCells(ECellType,ECellLoc,ECellArr)
# Assign the model('s) to the object
for item in model.iteritems():
# Convert numpy array
vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
vtkObj.Update()
return vtkObj
@staticmethod
def makeRenderWindow(ren):
renwin = vtk.vtkRenderWindow()
renwin.AddRenderer(ren)
iren = vtk.vtkRenderWindowInteractor()
iren.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
iren.SetRenderWindow(renwin)
return iren, renwin
@staticmethod
def closeRenderWindow(iren):
renwin = iren.GetRenderWindow()
renwin.Finalize()
iren.TerminateApp()
del iren, renwin
@staticmethod
def makeVTKActor(vtkObj):
""" Makes a vtk mapper and Actor"""
mapper = vtk.vtkDataSetMapper()
mapper.SetInput(vtkObj)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(0,0,0)
actor.GetProperty().SetRepresentationToWireframe()
return actor
@staticmethod
def makeVTKLODActor(vtkObj,clipper):
"""Make LOD vtk Actor"""
selectMapper = vtk.vtkDataSetMapper()
selectMapper.SetInputConnection(clipper.GetOutputPort())
selectMapper.SetScalarVisibility(1)
selectMapper.SetColorModeToMapScalars()
selectMapper.SetScalarModeToUseCellData()
selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange())
selectActor = vtk.vtkLODActor()
selectActor.SetMapper(selectMapper)
selectActor.GetProperty().SetEdgeColor(1,0.5,0)
selectActor.GetProperty().SetEdgeVisibility(0)
selectActor.VisibilityOn()
selectActor.SetScale(1.01, 1.01, 1.01)
return selectActor
@staticmethod
def setScalar2View(vtkObj,scalarName):
""" Sets the sclar to view """
useArr = vtkObj.GetCellData().GetArray(scalarName)
if useArr == None:
raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName))
vtkObj.GetCellData().SetActiveScalars(scalarName)
@staticmethod
def makeRectiVTKVOIThres(vtkObj,VOI,limits):
"""Make volume of interest and threshold for rectilinear grid."""
# Check for the input
cellCore = vtk.vtkExtractRectilinearGrid()
cellCore.SetVOI(VOI)
cellCore.SetInput(vtkObj)
cellThres = vtk.vtkThreshold()
cellThres.AllScalarsOn()
cellThres.SetInputConnection(cellCore.GetOutputPort())
cellThres.ThresholdBetween(limits[0],limits[1])
cellThres.Update()
return cellThres.GetOutput(), cellCore.GetOutput()
@staticmethod
def makeUnstructVTKVOIThres(vtkObj,extent,limits):
"""Make volume of interest and threshold for rectilinear grid."""
# Check for the input
cellCore = vtk.vtkExtractUnstructuredGrid()
cellCore.SetExtent(extent)
cellCore.SetInput(vtkObj)
cellThres = vtk.vtkThreshold()
cellThres.AllScalarsOn()
cellThres.SetInputConnection(cellCore.GetOutputPort())
cellThres.ThresholdBetween(limits[0],limits[1])
cellThres.Update()
return cellThres.GetOutput(), cellCore.GetOutput()
@staticmethod
def makePlaneClipper(vtkObj):
"""Makes a plane and clipper """
plane = vtk.vtkPlane()
clipper = vtk.vtkClipDataSet()
clipper.SetInputConnection(vtkObj.GetProducerPort())
clipper.SetClipFunction(plane)
clipper.InsideOutOff()
return clipper, plane
@staticmethod
def makePlaneWidget(vtkObj,iren,plane,actor):
"""Make an interactive planeWidget"""
# Callback function
def movePlane(obj, events):
obj.GetPlane(intPlane)
intActor.VisibilityOn()
# Associate the line widget with the interactor
planeWidget = vtk.vtkImplicitPlaneWidget()
planeWidget.SetInteractor(iren)
planeWidget.SetPlaceFactor(1.25)
planeWidget.SetInput(vtkObj)
planeWidget.PlaceWidget()
#planeWidget.AddObserver("InteractionEvent", movePlane)
planeWidget.SetScaleEnabled(0)
planeWidget.SetEnabled(1)
planeWidget.SetOutlineTranslation(0)
planeWidget.GetPlaneProperty().SetOpacity(0.1)
return planeWidget
@staticmethod
def startRenderWindow(iren):
""" Start a vtk rendering window"""
iren.Initialize()
renwin = iren.GetRenderWindow()
renwin.Render()
iren.Start()
# Simple write/read VTK xml model functions.
@staticmethod
def writeVTPFile(fileName,vtkPolyObject):
'''Function to write vtk polydata file (vtp).'''
polyWriter = vtk.vtkXMLPolyDataWriter()
polyWriter.SetInput(vtkPolyObject)
polyWriter.SetFileName(fileName)
polyWriter.Update()
@staticmethod
def writeVTUFile(fileName,vtkUnstructuredGrid):
'''Function to write vtk unstructured grid (vtu).'''
Writer = vtk.vtkXMLUnstructuredGridWriter()
Writer.SetInput(vtkUnstructuredGrid)
Writer.SetFileName(fileName)
Writer.Update()
@staticmethod
def writeVTRFile(fileName,vtkRectilinearGrid):
'''Function to write vtk rectilinear grid (vtr).'''
Writer = vtk.vtkXMLRectilinearGridWriter()
Writer.SetInput(vtkRectilinearGrid)
Writer.SetFileName(fileName)
Writer.Update()
@staticmethod
def writeVTSFile(fileName,vtkStructuredGrid):
'''Function to write vtk structured grid (vts).'''
Writer = vtk.vtkXMLStructuredGridWriter()
Writer.SetInput(vtkStructuredGrid)
Writer.SetFileName(fileName)
Writer.Update()
@staticmethod
def readVTSFile(fileName):
'''Function to read vtk structured grid (vts) and return a grid object.'''
Reader = vtk.vtkXMLStructuredGridReader()
Reader.SetFileName(fileName)
Reader.Update()
return Reader.GetOutput()
@staticmethod
def readVTUFile(fileName):
'''Function to read vtk structured grid (vtu) and return a grid object.'''
Reader = vtk.vtkXMLUnstructuredGridReader()
Reader.SetFileName(fileName)
Reader.Update()
return Reader.GetOutput()
@staticmethod
def readVTRFile(fileName):
'''Function to read vtk structured grid (vtr) and return a grid object.'''
Reader = vtk.vtkXMLRectilinearGridReader()
Reader.SetFileName(fileName)
Reader.Update()
return Reader.GetOutput()
@staticmethod
def readVTPFile(fileName):
'''Function to read vtk structured grid (vtp) and return a grid object.'''
Reader = vtk.vtkXMLPolyDataReader()
Reader.SetFileName(fileName)
Reader.Update()
return Reader.GetOutput()
SimPEG/visualize/vtk/vtkView.py
-350
View File
@@ -1,350 +0,0 @@
import numpy as np, matplotlib as mpl
try:
import vtk, vtk.util.numpy_support as npsup
#import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import
except Exception, e:
print 'VTK import error. Please ensure you have VTK installed to use this visualization package.'
import SimPEG as simpeg
class vtkView(object):
"""
Class for storing and view of SimPEG models in VTK (visualization toolkit).
Inputs:
:param mesh, SimPEG mesh.
:param propdict, dictionary of property models.
Can have these dictionary names:
'C' - cell model; 'F' - face model; 'E' - edge model; ('V' - vector field : NOT SUPPORTED)
The dictionary values are given as dictionaries with:
{'NameOfThePropertyModel': np.array of the properties}.
The property np.array has to be ordered in compliance with SimPEG standards.
::
Example of usages.
ToDo
"""
def __init__(self,mesh,propdict):
"""
"""
# Setup hidden properties, used for the visualization
self._ren = None
self._iren = None
self._renwin = None
self._core = None
self._viewobj = None
self._plane = None
self._clipper = None
self._widget = None
self._actor = None
self._lut = None
# Set vtk object containers
self._cells = None
self._faces = None
self._edges = None
self._vectors = None # Not implemented
# Set default values
self.name = 'VTK figure of SimPEG model'
# Error check the input mesh
if type(mesh).__name__ != 'TensorMesh':
raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh))
# Set the mesh
self._mesh = mesh
# Read the property dictionary
self._readPropertyDictionary(propdict)
# Set/Get properties
@property
def cmap(self):
''' Colormap to use in vtkView. Colormap is a matplotlib cmap(cm) array, has to be uint8(use flag bytes=True during cmap generation).'''
if getattr(self,'_cmap',None) is None:
# Set default
self._cmap = mpl.cm.hsv(np.arange(0.,1.,0.05),bytes=True)
return self._cmap
@cmap.setter
def cmap(self,value):
if value.min() > 0 or value.max() < 255 or value.shape[1] != 4 or value.dtype != np.uint8:
raise Exception('Input not an allowed array.\n Use matplotlib.cm to generate an array of size [nrColors,4] and dtype = uint8(flag bytes=True).')
self._cmap = value
@property
def range(self):
''' Range of the colors in vtkView.'''
if getattr(self,'_range',None) is None:
self._range = np.array(self._getActiveVTKobj().GetArray(self.viewprop.values()[0]).GetRange())
return self._range
@range.setter
def range(self,value):
if type(value) not in [tuple, list, np.ndarray] or len(value) != 2 or np.array(value).dtype is not np.dtype('float'):
raise Exception('Input not in correct format. \n Has to be a list, tuple or np.arry of 2 floats.')
self._range = np.array(value)
@property
def extent(self):
''' Extent of the sub-domain of the model to view'''
if getattr(self,'_extent',None) is None:
self._extent = [0,self._mesh.nCx-1,0,self._mesh.nCy-1,0,self._mesh.nCz-1]
return self._extent
@extent.setter
def extent(self,value):
import warnings
# Error check
valnp = np.array(value,dtype=int)
if valnp.dtype != int or len(valnp) != 6:
raise Exception('.extent has to be list or nparray of 6 integers.')
# Test the range of the values
loB = np.zeros(3,dtype=int)
upB = np.array(self._mesh.nCv - np.ones(3),dtype=int)
# Test the bounds
change = 0
# Test for lower bounds, can't be smaller the 0
tlb = valnp[::2] < loB
if tlb.any():
valnp[::2][tlb] = loB[tlb]
change = 1
warnings.warn('Lower bounds smaller then 0')
# Test for lower bounds, can't be larger then upB
tlub = valnp[::2] > upB
if tlub.any():
valnp[::2][tlub] = upB[tlub] - 1
change = 1
warnings.warn('Lower bounds larger then uppermost bounds')
# Test for upper bounds, can't be larger the extent of the mesh
tub = valnp[1::2] > upB
if tub.any():
valnp[1::2][tub] = upB[tub]
change = 1
warnings.warn('Upper bounds greater then number of cells')
# Test if lower is smaller the upper
tgt = valnp[::2] > valnp[1::2]
if tgt.any():
valnp[1::2][tgt] = valnp[::2][tgt] + 1
change = 1
warnings.warn('Lower bounds greater the Upper bounds')
# Print a warning
if change:
warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))
# Set extent
self._extent = valnp
@property
def limits(self):
''' Lower and upper limits (cutoffs) of the values to view. '''
return getattr(self,'_limits',None)
@limits.setter
def limits(self,value):
if value is None:
self._limits = None
else:
valnp = np.array(value)
if valnp.dtype != float or len(valnp) != 2:
raise Exception('.limits has to be list or numpy array of 2 floats.')
self._limits = valnp
@property
def viewprop(self):
''' Controls the property that will be viewed.'''
if getattr(self,'_viewprop',None) is None:
self._viewprop = {'C':0} # Name of the type and Int order of the array or name of the vector.
return self._viewprop
@viewprop.setter
def viewprop(self,value):
if type(value) != dict:
raise Exception('{:s} has to be a python dictionary containing property type and name index. ')
if len(value) > 1:
raise Exception('Too many input items in the viewprop dictionary')
if value.keys()[0] not in ['C','F','E']:
raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0]))
if not(type(self.viewprop.values()[0]) is int or type(self.viewprop.values()[0]) is str):
raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be integer or a string with the index.')
self._viewprop = value
def _getActiveVTKobj(self):
"""
Finds the active VTK object.
"""
if self.viewprop.keys()[0] is 'C':
vtkCellData = self._cells.GetCellData()
elif self.viewprop.keys()[0] is 'F':
vtkCellData = self._faces.GetCellData()
elif self.viewprop.keys()[0] is 'E':
vtkCellData = self._edges.GetCellData()
return vtkCellData
def _getActiveArrayName(self):
"""
Finds the name of the active array.
"""
actArr = self.viewprop.values()[0]
if type(actArr) is str:
activeName = actArr
elif type(actArr) is int:
activeName = self._getActiveVTKobj().GetArrayName(actArr)
return activeName
def _readPropertyDictionary(self,propdict):
"""
Reads the property and assigns to the object
"""
import SimPEG.visualize.vtk.vtkTools as vtkSP
# Test the property dictionary
if type(propdict) != dict:
raise Exception('{:s} has to be a python dictionary containing property models. ')
if len(propdict) > 4:
raise Exception('Too many input items in the property dictionary')
for propitem in propdict.iteritems():
if propitem[0] in ['C','F','E']:
if propitem[0] == 'C':
self._cells = vtkSP.makeCellVTKObject(self._mesh,propitem[1])
if propitem[0] == 'F':
self._faces = vtkSP.makeFaceVTKObject(self._mesh,propitem[1])
if propitem[0] == 'E':
self._edges = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1])
else:
raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0]))
def Show(self):
"""
Open the VTK figure window and show the mesh.
"""
#vtkSP = simpeg.visualize.vtk.vtkTools
import SimPEG.visualize.vtk.vtkTools as vtkSP
# Make a renderer
self._ren = vtk.vtkRenderer()
# Make renderwindow. Returns the interactor.
self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren)
# Set the active scalar.
if type(self.viewprop.values()[0]) == int:
actScalar = self._getActiveVTKobj().GetArrayName(self.viewprop.values()[0])
elif type(self.viewprop.values()[0]) == str:
actScalar = self.viewprop.values()[0]
else :
raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be interger or a string.')
self._getActiveVTKobj().SetActiveScalars(actScalar)
# Sort out the actor
imageType = self.viewprop.keys()[0]
if imageType == 'C':
if self.limits is None:
self.limits = self._cells.GetCellData().GetArray(self.viewprop.values()[0]).GetRange()
self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cells,self.extent,self.limits)
elif imageType == 'F':
if self.limits is None:
self.limits = self._faces.GetCellData().GetArray(self.viewprop.values()[0]).GetRange()
extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ]
self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._faces,extent,self.limits)
elif imageType == 'E':
if self.limits is None:
self.limits = self._edges.GetCellData().GetArray(self.viewprop.values()[0]).GetRange()
extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ]
self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edges,extent,self.limits)
else:
raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType))
#self._vtkobj.GetCellData().SetActiveScalars(actScalar)
# Set up the plane, clipper and the user interaction.
global intPlane, intActor
self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj)
intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper)
self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor)
# Callback function
self._plane = intPlane
self._actor = intActor
def movePlane(obj, events):
global intPlane, intActor
obj.GetPlane(intPlane)
intActor.VisibilityOn()
self._widget.AddObserver("InteractionEvent",movePlane)
lut = vtk.vtkLookupTable()
lut.SetNumberOfColors(len(self.cmap))
lut.SetTable(npsup.numpy_to_vtk(self.cmap))
lut.Build()
self._lut = lut
scalarBar = vtk.vtkScalarBarActor()
scalarBar.SetLookupTable(lut)
scalarBar.SetTitle(self._getActiveArrayName())
scalarBar.GetPositionCoordinate().SetCoordinateSystemToNormalizedViewport()
scalarBar.GetPositionCoordinate().SetValue(0.1,0.01)
scalarBar.SetOrientationToHorizontal()
scalarBar.SetWidth(0.8)
scalarBar.SetHeight(0.17)
self._actor.GetMapper().SetScalarRange(self.range)
self._actor.GetMapper().SetLookupTable(lut)
# Set renderer options
self._ren.SetBackground(.5,.5,.5)
self._ren.AddActor(self._actor)
self._ren.AddActor2D(scalarBar)
self._renwin.SetSize(450,450)
# Start the render Window
vtkSP.startRenderWindow(self._iren)
# Close the window when exited
vtkSP.closeRenderWindow(self._iren)
del self._iren, self._renwin
if __name__ == '__main__':
#Make a mesh and model
x0 = np.zeros(3)
h1 = np.ones(60)*50
h2 = np.ones(60)*100
h3 = np.ones(50)*200
mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)
# Make a models that correspond to the cells, faces and edges.
t = np.ones(mesh.nC)
t[10000:50000] = 100
t[100000:120000] = 100
t[100000:120000] = 50
models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t, 'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,mesh.nF),'AllOnce':np.ones(mesh.nF)},'E':{'Test':np.arange(0,mesh.nE),'AllOnce':np.ones(mesh.nE)}}
# Make the vtk viewer object.
vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)
# Set the .viewprop for which model to view
vtkViewer.viewprop = {'F':'Test'}
# Show the image
vtkViewer.Show()
# Set subset of the mesh to view (remove padding)
vtkViewer.extent = [4,14,0,7,0,3]
vtkViewer.Show()
# Change viewing property
vtkViewer.viewprop = {'C':'Model'}
# Set the color range
# Reset extent.
vtkViewer.extent = [-1,1000,-1,1000,-1,1000]
vtkViewer.range = [0.,100.]
vtkViewer.Show()
# Change color scale, has to be set to bytes=True.
vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)
vtkViewer.Show()
# Set limits of values to view
vtkViewer.limits = [5.0,100.0]
vtkViewer.Show()
SimPEG/Examples/Linear.py → Tutorials/Linear.py
+1 -1
View File
@@ -9,7 +9,7 @@ class LinearProblem(Problem.BaseProblem):
Problem.BaseProblem.__init__(self, mesh, model, **kwargs)
self.G = G
def field(self, m, u=None):
def fields(self, m, u=None):
return self.G.dot(m)
def J(self, m, v, u=None):