diff --git a/.travis.yml b/.travis.yml index 1918d175..1094912b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -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 diff --git a/SimPEG/Data.py b/SimPEG/Data.py index fa37b09a..0290be36 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Data.py @@ -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 diff --git a/SimPEG/Examples/DC.py b/SimPEG/Examples/DC.py deleted file mode 100644 index 074f5e3d..00000000 --- a/SimPEG/Examples/DC.py +++ /dev/null @@ -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() - - - - - - diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py deleted file mode 100644 index a5c37345..00000000 --- a/SimPEG/Examples/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -import DC -import Linear diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 56a00f8f..64445ab2 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -7,7 +7,7 @@ class BaseInversion(object): """BaseInversion(objFunc, opt, **kwargs) """ - __metaclass__ = Utils.Save.Savable + __metaclass__ = Utils.SimPEGMetaClass name = 'BaseInversion' diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index 672b5c02..2dd82808 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -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): """ diff --git a/SimPEG/Mesh/Cyl1DMesh.py b/SimPEG/Mesh/Cyl1DMesh.py index 2b291c54..418a866b 100644 --- a/SimPEG/Mesh/Cyl1DMesh.py +++ b/SimPEG/Mesh/Cyl1DMesh.py @@ -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 #################################################### diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 492d628e..87b6aac3 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -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) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 9e7c1d93..43d6c31f 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -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) diff --git a/SimPEG/Mesh/LogicallyOrthogonalMesh.py b/SimPEG/Mesh/LogicallyOrthogonalMesh.py index e39c0f54..a8b4207f 100644 --- a/SimPEG/Mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/Mesh/LogicallyOrthogonalMesh.py @@ -26,7 +26,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): M.plotGrid(showIt=True) """ - __metaclass__ = Utils.Save.Savable + __metaclass__ = Utils.SimPEGMetaClass _meshType = 'LOM' diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 341cfba7..b2981bdd 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -33,7 +33,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): """ - __metaclass__ = Utils.Save.Savable + __metaclass__ = Utils.SimPEGMetaClass _meshType = 'TENSOR' diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index 7446aecd..a664501a 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -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) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py new file mode 100644 index 00000000..3d2c2794 --- /dev/null +++ b/SimPEG/Mesh/TreeMesh.py @@ -0,0 +1,1116 @@ +from SimPEG import np, sp, Utils, Solver +from BaseMesh import BaseMesh +from InnerProducts import InnerProducts +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.colors as colors +import matplotlib.cm as cmx + + + +def SortByX0(): + eps = 1e-7 + def mycmp(c1,c2): + if c1.x0.size == 2: + if np.abs(c1.x0[1] - c2.x0[1]) < eps: + return c1.x0[0] - c2.x0[0] + return c1.x0[1] - c2.x0[1] + elif c1.x0.size == 3: + if np.abs(c1.x0[2] - c2.x0[2]) < eps: + if np.abs(c1.x0[1] - c2.x0[1]) < eps: + return c1.x0[0] - c2.x0[0] + return c1.x0[1] - c2.x0[1] + return c1.x0[2] - c2.x0[2] + + class K(object): + def __init__(self, obj, *args): + self.obj = obj + def __lt__(self, other): + return mycmp(self.obj, other.obj) < 0 + def __gt__(self, other): + return mycmp(self.obj, other.obj) > 0 + def __eq__(self, other): + return mycmp(self.obj, other.obj) == 0 + def __le__(self, other): + return mycmp(self.obj, other.obj) <= 0 + def __ge__(self, other): + return mycmp(self.obj, other.obj) >= 0 + def __ne__(self, other): + return mycmp(self.obj, other.obj) != 0 + return K + + +class TreeNode(object): + """docstring for TreeNode""" + + __slots__ = ['x0', 'num'] + + def __init__(self, mesh, x0=[0,0]): + self.x0 = np.array(x0, dtype=float) + mesh.nodes.add(self) + + @property + def center(self): return self.x0 + +class TreeEdge(object): + """docstring for TreeEdge""" + + __slots__ = ['mesh', 'children', 'depth', 'x0', 'num', 'edgeType', 'sz', 'node0', 'node1'] + + def __init__(self, mesh, x0=[0,0], edgeType=None, sz=[1,], depth=0, + node0=None, node1=None): + self.mesh = mesh + self.depth = depth + + self.x0 = x0 + self.sz = sz + self.edgeType = edgeType + + mesh.edges.add(self) + if edgeType is 'x': mesh.edgesX.add(self) + elif edgeType is 'y': mesh.edgesY.add(self) + elif edgeType is 'z': mesh.edgesZ.add(self) + + self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=self.x0) + self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent*self.sz[0]) + + @property + def isleaf(self): return getattr(self, 'children', None) is None + + def refine(self): + if not self.isleaf: return + self.mesh.isNumbered = False + + self.children = np.empty(2,dtype=TreeFace) + # Create refined x0's + x0r_0 = self.x0 + x0r_1 = self.x0+0.5*self.tangent*self.sz + self.children[0] = TreeEdge(self.mesh, x0=x0r_0, edgeType=self.edgeType, sz=0.5*self.sz, depth=self.depth+1, node0=self.node0) + self.children[1] = TreeEdge(self.mesh, x0=x0r_1, edgeType=self.edgeType, sz=0.5*self.sz, depth=self.depth+1, node0=self.children[0].node1, node1=self.node1) + self.mesh.edges.remove(self) + if self.edgeType is 'x': + self.mesh.edgesX.remove(self) + elif self.edgeType is 'y': + self.mesh.edgesY.remove(self) + elif self.edgeType is 'z': + self.mesh.edgesZ.remove(self) + + @property + def tangent(self): + if self.edgeType is 'x': return np.r_[1.,0,0] + elif self.edgeType is 'y': return np.r_[0,1.,0] + elif self.edgeType is 'z': return np.r_[0,0,1.] + + def plotGrid(self, ax, text=False, lineOpts={'color':'r', 'ls': '-'}): + line = np.c_[self.node0.x0, self.node1.x0].T + ax.plot(line[:,0], line[:,1], zs=line[:,2], **lineOpts) + + @property + def center(self): + return 0.5*(self.node0.x0 + self.node1.x0) + + @property + def length(self): + return np.sqrt(((self.node1.x0 - self.node0.x0)**2).sum()) + + @property + def index(self): + if self.isleaf: return [self.num] + l = [edge.index for edge in self.children.flatten(order='F')] + # Flatten the list + # e.g. + # [[1,3],[4]] --> [1, 3, 4] + return [item for sublist in l for item in sublist] + +class TreeFace(object): + """docstring for TreeFace""" + + __slots__ = ['mesh', 'children', 'depth', 'num', 'faceType', 'sz', 'node0', 'node1', 'node2', 'node3', 'edge0', 'edge1', 'edge2', 'edge3', '_tangent0', '_tangent1'] + + def __init__(self, mesh, x0=[0,0], faceType=None, sz=[1,], depth=0, + node0=None, node1=None, + edge0=None, edge1=None, edge2=None, edge3=None): + + self.mesh = mesh + self.depth = depth + + self.faceType = faceType + self.sz = sz + + mesh.faces.add(self) + if faceType is 'x': mesh.facesX.add(self) + elif faceType is 'y': mesh.facesY.add(self) + elif faceType is 'z': mesh.facesZ.add(self) + if self.dim == 2: + # Add the nodes: + self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=x0) + self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=x0 + self.tangent0*self.sz[0]) + if self.dim == 3: + #TODO: Change this to edges + + # + # 2___________3 + # | e1 | + # | | + # e2| x |e3 t1 + # | | ^ + # |___________| |___> t0 + # 0 e0 1 + # + + N = {} + n0 = getattr(edge0, 'node0', None) or getattr(edge2, 'node0', None) + n1 = getattr(edge0, 'node1', None) or getattr(edge3, 'node0', None) + n2 = getattr(edge1, 'node0', None) or getattr(edge2, 'node1', None) + n3 = getattr(edge1, 'node1', None) or getattr(edge3, 'node1', None) + + eType = ['x', 'y'] if self.faceType == 'z' else ['x', 'z'] if self.faceType == 'y' else ['y', 'z'] + + e0 = edge0 if isinstance(edge0,TreeEdge) else TreeEdge(mesh, x0=x0, edgeType=eType[0], sz=np.r_[sz[0]], depth=depth, node0=n0, node1=n1) + n0, n1 = e0.node0, e0.node1 + + e1 = edge1 if isinstance(edge1,TreeEdge) else TreeEdge(mesh, x0=x0 + self.tangent1*self.sz[1], edgeType=eType[0], sz=np.r_[sz[0]], depth=depth, node0=n2, node1=n3) + n2, n3 = e1.node0, e1.node1 + + e2 = edge2 if isinstance(edge2,TreeEdge) else TreeEdge(mesh, x0=x0, edgeType=eType[1], sz=np.r_[sz[1]], depth=depth, node0=n0, node1=n2) + n0, n2 = e2.node0, e2.node1 + + e3 = edge3 if isinstance(edge3,TreeEdge) else TreeEdge(mesh, x0=x0 + self.tangent0*self.sz[0], edgeType=eType[1], sz=np.r_[sz[1]], depth=depth, node0=n1, node1=n3) + n1, n3 = e3.node0, e3.node1 + + # self.nodes = N + self.node0, self.node1, self.node2, self.node3 = n0, n1, n2, n3 + self.edge0, self.edge1, self.edge2, self.edge3 = e0, e1, e2, e3 + # self.edges = {'e0':e0, 'e1':e1, 'e2':e2, 'e3':e3} + + @property + def dim(self): return self.mesh.dim + + @property + def x0(self): return self.node0.x0 + + @property + def isleaf(self): return getattr(self, 'children', None) is None + + @property + def branchdepth(self): + if self.isleaf: + return self.depth + else: + return np.max([node.branchdepth for node in self.children.flatten('F')]) + + @property + def tangent0(self): + if getattr(self,'_tangent0',None) is None: + if self.faceType is 'x': t = np.r_[0,1.,0] + elif self.faceType is 'y': t = np.r_[1.,0,0] + elif self.faceType is 'z': t = np.r_[1.,0,0] + self._tangent0 = t[:self.dim] + return self._tangent0 + + @property + def tangent1(self): + if self.dim == 2: return + if getattr(self,'_tangent1',None) is None: + if self.faceType is 'x': t = np.r_[0,0,1.] + elif self.faceType is 'y': t = np.r_[0,0,1.] + elif self.faceType is 'z': t = np.r_[0,1.,0] + self._tangent1 = t + return self._tangent1 + + @property + def normal(self): + if self.faceType is 'x': n = np.r_[1.,0,0] + elif self.faceType is 'y': n = np.r_[0,1.,0] + elif self.faceType is 'z': n = np.r_[0,0,1.] + return n[:self.dim] + + @property + def index(self): + if self.isleaf: return [self.num] + l = [face.index for face in self.children.flatten(order='F')] + # Flatten the list + # e.g. + # [[1,3],[4]] --> [1, 3, 4] + return [item for sublist in l for item in sublist] + + @property + def area(self): + """area of the face""" + return self.sz.prod() + + @property + def length(self): + if self.dim == 3: raise Exception('face.length is not defined for 2D face') + return np.sqrt(((self.node1.x0 - self.node0.x0)**2).sum()) + + def refine(self): + if not self.isleaf: return + self.mesh.isNumbered = False + if self.dim == 2: + self._refine2D() + elif self.dim == 3: + self._refine3D() + + def _refine2D(self): + self.children = np.empty(2,dtype=TreeFace) + # Create refined x0's + x0r_0 = self.x0 + x0r_1 = self.x0+0.5*self.tangent0*self.sz + self.children[0] = TreeFace(self.mesh, x0=x0r_0, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, node0=self.node0) + self.children[1] = TreeFace(self.mesh, x0=x0r_1, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, node0=self.children[0].node1, node1=self.node1) + self.mesh.faces.remove(self) + if self.faceType is 'x': + self.mesh.facesX.remove(self) + elif self.faceType is 'y': + self.mesh.facesY.remove(self) + + def _refine3D(self): + # + # 2_______________3 _______________ + # | e1--> | | | | + # ^ | | ^ | (0,1) | (1,1) | + # | | | | | | | + # | | x | | ---> |-------+-------| + # e2 | | e3 | | | + # | | | (0,0) | (1,0) | + # |_______________| |_______|_______| + # 0 e0--> 1 + + + order = [{'c':[0,0], + 'e0': ('p', 'e0', [0]), 'e1': 'new' , + 'e2': ('p', 'e2', [0]), 'e3': 'new' }, + {'c':[1,0], + 'e0': ('p', 'e0', [1]), 'e1': 'new' , + 'e2': ('c', 'e3', [0,0]), 'e3': ('p', 'e3', [0])}, + {'c':[0,1], + 'e0': ('c', 'e1', [0,0]), 'e1': ('p', 'e1', [0]), + 'e2': ('p', 'e2', [1]), 'e3': 'new' }, + {'c':[1,1], + 'e0': ('c', 'e1', [1,0]), 'e1': ('p', 'e1', [1]), + 'e2': ('c', 'e3', [0,1]), 'e3': ('p', 'e3', [1])}] + + def getEdge(pointer): + if pointer is 'new': return + if pointer[0] == 'p': + return getattr(self, 'edg' + pointer[1]).children[pointer[2][0]] + if pointer[0] == 'c': + f = self.children[pointer[2][0],pointer[2][1]] + return getattr(f, 'edg' + pointer[1]) + + self.children = np.empty((2,2), dtype=TreeFace) + + for edge in [self.edge0, self.edge1, self.edge2, self.edge3]: + edge.refine() + + for O in order: + i, j = O['c'] + x0r = self.x0 + 0.5*i*self.tangent0*self.sz[0] + 0.5*j*self.tangent1*self.sz[1] + e0, e1, e2, e3 = getEdge(O['e0']), getEdge(O['e1']), getEdge(O['e2']), getEdge(O['e3']) + self.children[i,j] = TreeFace(self.mesh, x0=x0r, faceType=self.faceType, depth=self.depth+1, sz=0.5*self.sz, edge0=e0, edge1=e1, edge2=e2, edge3=e3) + + self.mesh.faces.remove(self) + if self.faceType is 'x': + self.mesh.facesX.remove(self) + elif self.faceType is 'y': + self.mesh.facesY.remove(self) + elif self.faceType is 'z': + self.mesh.facesZ.remove(self) + + def plotGrid(self, ax, text=True): + if not self.isleaf: return + if self.dim == 2: + line = np.c_[self.node0.x0, self.node1.x0].T + ax.plot(line[:,0], line[:,1],'r-') + if text: ax.text(self.center[0], self.center[1],self.num) + elif self.dim == 3: + if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) + + @property + def center(self): + if self.dim == 2: + return self.x0 + 0.5*self.tangent0*self.sz[0] + elif self.dim == 3: + return self.x0 + 0.5*self.tangent0*self.sz[0] + 0.5*self.tangent1*self.sz[1] + + +class TreeCell(object): + + __slots__ = ['mesh', 'children', 'depth', 'num', 'sz', + 'node0', 'node1', 'node2', 'node3', + 'node4', 'node5', 'node6', 'node7', + 'fXm', 'fXp', 'fYm', 'fYp', 'fZm', 'fZp', + 'eX0','eX1','eX2','eX3', + 'eY0','eY1','eY2','eY3', + 'eZ0','eZ1','eZ2','eZ3'] + + def __init__(self, mesh, x0=[0,0], depth=0, sz=[1,1], + fXm=None, fXp=None, + fYm=None, fYp=None, + fZm=None, fZp=None): + + self.mesh = mesh + self.depth = depth + + self.sz = np.array(sz, dtype=float) + if self.dim == 2: + # + # 2___________3 + # | fYp | + # | | + # fXm| x |fXp y + # | | ^ + # |___________| |___> x + # 0 fYm 1 + # + n0 = getattr(fXm, 'node0', None) or getattr(fYm, 'node0', None) + n1 = getattr(fXp, 'node0', None) or getattr(fYm, 'node1', None) + n2 = getattr(fXm, 'node1', None) or getattr(fYp, 'node0', None) + n3 = getattr(fXp, 'node1', None) or getattr(fYp, 'node1', None) + + self.fXm = fXm if isinstance(fXm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='x', sz=np.r_[sz[1]], depth=depth, node0=n0, node1=n2) + n0, n2 = self.fXm.node0, self.fXm.node1 + + self.fXp = fXp if isinstance(fXp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] ], faceType='x', sz=np.r_[sz[1]], depth=depth, node0=n1, node1=n3) + n1, n3 = self.fXp.node0, self.fXp.node1 + + self.fYm = fYm if isinstance(fYm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='y', sz=np.r_[sz[0]], depth=depth, node0=n0, node1=n1) + n0, n1 = self.fYm.node0, self.fYm.node1 + + self.fYp = fYp if isinstance(fYp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1]], faceType='y', sz=np.r_[sz[0]], depth=depth, node0=n2, node1=n3) + n2, n3 = self.fYp.node0, self.fYp.node1 + + self.node0, self.node1, self.node2, self.node3 = n0, n1, n2, n3 + + elif self.dim == 3: + # fZp + # | + # 6 ------eX3------ 7 + # /| | / | + # /eZ2 . / eZ3 + # eY2 | fYp eY3 | + # / | / fXp| + # 4 ------eX2----- 5 | + # |fXm 2 -----eX1--|---- 3 z + # eZ0 / | eY1 ^ y + # | eY0 . fYm eZ1 / | / + # | / | | / | / + # 0 ------eX0------1 o----> x + # | + # fZm + # + # + # fX fY fZ + # 2___________3 2___________3 2___________3 + # | e1 | | e1 | | e1 | + # | | | | | | + # e2 | x | e3 z e2 | x | e3 z e2 | x | e3 y + # | | ^ | | ^ | | ^ + # |___________| |___> y |___________| |___> x |___________| |___> x + # 0 e0 1 0 e0 1 0 e0 1 + # + # Mapping Nodes: numOnFace > numOnCell + # + # fXm 0>0, 1>2, 2>4, 3>6 fYm 0>0, 1>1, 2>4, 3>5 fZm 0>0, 1>1, 2>2, 3>3 + # fXp 0>1, 1>3, 2>5, 3>7 fYp 0>2, 1>3, 2>6, 3>7 fZp 0>4, 1>5, 2>6, 3>7 + + def getEdge(face, key): + if face is None: return + return getattr(face, key) + + E = {} + eX0 = getEdge(fYm, 'edge0') or getEdge(fZm, 'edge0') + eX1 = getEdge(fYp, 'edge0') or getEdge(fZm, 'edge1') + eX2 = getEdge(fYm, 'edge1') or getEdge(fZp, 'edge0') + eX3 = getEdge(fYp, 'edge1') or getEdge(fZp, 'edge1') + + eY0 = getEdge(fXm, 'edge0') or getEdge(fZm, 'edge2') + eY1 = getEdge(fXp, 'edge0') or getEdge(fZm, 'edge3') + eY2 = getEdge(fXm, 'edge1') or getEdge(fZp, 'edge2') + eY3 = getEdge(fXp, 'edge1') or getEdge(fZp, 'edge3') + + eZ0 = getEdge(fXm, 'edge2') or getEdge(fYm, 'edge2') + eZ1 = getEdge(fXp, 'edge2') or getEdge(fYm, 'edge3') + eZ2 = getEdge(fXm, 'edge3') or getEdge(fYp, 'edge2') + eZ3 = getEdge(fXp, 'edge3') or getEdge(fYp, 'edge3') + + + self.fXm = fXm if isinstance(fXm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, edge0=eY0, edge1=eY2, edge2=eZ0, edge3=eZ2) + eY0, eY2, eZ0, eZ2 = self.fXm.edge0, self.fXm.edge1, self.fXm.edge2, self.fXm.edge3 + + self.fXp = fXp if isinstance(fXp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, edge0=eY1, edge1=eY3, edge2=eZ1, edge3=eZ3) + eY1, eY3, eZ1, eZ3 = self.fXp.edge0, self.fXp.edge1, self.fXp.edge2, self.fXp.edge3 + + self.fYm = fYm if isinstance(fYm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, edge0=eX0, edge1=eX2, edge2=eZ0, edge3=eZ1) + eX0, eX2, eZ0, eZ1 = self.fYm.edge0, self.fYm.edge1, self.fYm.edge2, self.fYm.edge3 + + self.fYp = fYp if isinstance(fYp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1], x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, edge0=eX1, edge1=eX3, edge2=eZ2, edge3=eZ3) + eX1, eX3, eZ2, eZ3 = self.fYp.edge0, self.fYp.edge1, self.fYp.edge2, self.fYp.edge3 + + self.fZm = fZm if isinstance(fZm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, edge0=eX0, edge1=eX1, edge2=eY0, edge3=eY1) + eX0, eX1, eY0, eY1 = self.fZm.edge0, self.fZm.edge1, self.fZm.edge2, self.fZm.edge3 + + self.fZp = fZp if isinstance(fZp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2]+sz[2]], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, edge0=eX2, edge1=eX3, edge2=eY2, edge3=eY3) + eX2, eX3, eY2, eY3 = self.fZp.edge0, self.fZp.edge1, self.fZp.edge2, self.fZp.edge3 + + self.eX0, self.eX1, self.eX2, self.eX3, self.eY0, self.eY1, self.eY2, self.eY3, self.eZ0, self.eZ1, self.eZ2, self.eZ3 = eX0, eX1, eX2, eX3, eY0, eY1, eY2, eY3, eZ0, eZ1, eZ2, eZ3 + self.node0, self.node1, self.node2, self.node3, self.node4, self.node5, self.node6, self.node7 = self.fZm.node0, self.fZm.node1, self.fZm.node2, self.fZm.node3, self.fZp.node0, self.fZp.node1, self.fZp.node2, self.fZp.node3 + + mesh.cells.add(self) + + @property + def x0(self): return self.node0.x0 + + @property + def center(self): return self.x0 + 0.5*self.sz + + @property + def dim(self): return self.mesh.dim + + @property + def faceDict(self): + d = {"fXm":self.fXm, "fXp":self.fXp, "fYm":self.fYm, "fYp":self.fYp} + if self.dim == 3: + d["fZm"] = self.fZm + d["fZp"] = self.fZp + return d + + @property + def edgeDict(self): + if self.dim == 2: return None + return {'eX0': self.eX0, 'eX1': self.eX1, 'eX2': self.eX2, 'eX3': self.eX3, 'eY0': self.eY0, 'eY1': self.eY1, 'eY2': self.eY2, 'eY3': self.eY3, 'eZ0': self.eZ0, 'eZ1': self.eZ1, 'eZ2': self.eZ2, 'eZ3': self.eZ3} + + @property + def faceList(self): + l = [self.fXm, self.fXp, self.fYm, self.fYp] + if self.dim == 3: + l += [self.fZm, self.fZp] + return l + + @property + def edgeList(self): + if self.dim == 2: return None + return [self.eX0, self.eX1, self.eX2, self.eX3, self.eY0, self.eY1, self.eY2, self.eY3, self.eZ0, self.eZ1, self.eZ2, self.eZ3] + + @property + def isleaf(self): return getattr(self, 'children', None) is None + + def refine(self, function=None): + if not self.isleaf and function is None: return + + if function is not None: + do = function(self.center) > self.depth + if not do: return + + if self.dim == 2: + self._refine2D() + elif self.dim == 3: + self._refine3D() + + # pass the refine function to the children + if function is not None: + for child in self.children.flatten(): + child.refine(function) + + def _refine2D(self): + + self.mesh.isNumbered = False + + self.children = np.empty((2,2), dtype=TreeCell) + x0, sz = self.x0, self.sz + + for face in self.faceList: + face.refine() + + order = [{'c':[0,0], + 'fXm': ('p', 'fXm', [0]), 'fXp': 'new' , + 'fYm': ('p', 'fYm', [0]), 'fYp': 'new' }, + {'c':[1,0], + 'fXm': ('c', 'fXp', [0,0]), 'fXp': ('p', 'fXp', [0]), + 'fYm': ('p', 'fYm', [1]), 'fYp': 'new' }, + {'c':[0,1], + 'fXm': ('p', 'fXm', [1]), 'fXp': 'new' , + 'fYm': ('c', 'fYp', [0,0]), 'fYp': ('p', 'fYp', [0])}, + {'c':[1,1], + 'fXm': ('c', 'fXp', [0,1]), 'fXp': ('p', 'fXp', [1]), + 'fYm': ('c', 'fYp', [1,0]), 'fYp': ('p', 'fYp', [1])}] + + def getFace(pointer): + if pointer is 'new': return None + if pointer[0] == 'p': + return self.faceDict[pointer[1]].children[pointer[2][0],] + if pointer[0] == 'c': + return self.children[pointer[2][0],pointer[2][1]].faceDict[pointer[1]] + + for O in order: + i, j = O['c'] + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1]] + fXm, fXp, fYm, fYp = getFace(O['fXm']), getFace(O['fXp']), getFace(O['fYm']), getFace(O['fYp']) + self.children[i,j] = TreeCell(self.mesh, x0=x0r, depth=self.depth+1, sz=0.5*sz, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp) + + self.mesh.cells.remove(self) + + + def _refine3D(self): + # .----------------.----------------. + # /| /| /| + # / | / | / | + # / | 011 / | 111 / | + # / | / | / | + # .----------------.----+-----------. | + # /| . ---------/|----.----------/|----. + # / | /| / | /| / | /| + # / | / | 001 / | / | 101 / | / | + # / | / | / | / | / | / | + # . -------------- .----------------. |/ | + # | . ---+------|----.----+------|----. | + # | /| .______|___/|____.______|___/|____. + # | / | / 010 | / | / 110| / | / + # | / | / | / | / | / | / + # . ---+---------- . ---+---------- . | / + # | |/ | |/ | |/ z + # | . ----------|----.-----------|----. ^ y + # | / 000 | / 100 | / | / + # | / | / | / | / + # | / | / | / o----> x + # . -------------- . -------------- . + # + # + # Face Refinement: + # + # 2_______________3 _______________ + # | | | | | + # ^ | | | (0,1) | (1,1) | + # | | | | | | + # | | x | ---> |-------+-------| + # t1 | | | | | + # | | | (0,0) | (1,0) | + # |_______________| |_______|_______| + # 0 t0--> 1 + + + order = [{'c':[0,0,0], + 'fXm': ('p', 'fXm', [0,0]), 'fXp': 'new' , + 'fYm': ('p', 'fYm', [0,0]), 'fYp': 'new' , + 'fZm': ('p', 'fZm', [0,0]), 'fZp': 'new' ,}, + {'c':[1,0,0], + 'fXm': ('c', 'fXp', [0,0,0]), 'fXp': ('p', 'fXp', [0,0]), + 'fYm': ('p', 'fYm', [1,0]), 'fYp': 'new' , + 'fZm': ('p', 'fZm', [1,0]), 'fZp': 'new' }, + {'c':[0,1,0], + 'fXm': ('p', 'fXm', [1,0]), 'fXp': 'new' , + 'fYm': ('c', 'fYp', [0,0,0]), 'fYp': ('p', 'fYp', [0,0]), + 'fZm': ('p', 'fZm', [0,1]), 'fZp': 'new' }, + {'c':[1,1,0], + 'fXm': ('c', 'fXp', [0,1,0]), 'fXp': ('p', 'fXp', [1,0]), + 'fYm': ('c', 'fYp', [1,0,0]), 'fYp': ('p', 'fYp', [1,0]), + 'fZm': ('p', 'fZm', [1,1]), 'fZp': 'new' }, + {'c':[0,0,1], + 'fXm': ('p', 'fXm', [0,1]), 'fXp': 'new' , + 'fYm': ('p', 'fYm', [0,1]), 'fYp': 'new' , + 'fZm': ('c', 'fZp', [0,0,0]), 'fZp': ('p', 'fZp', [0,0])}, + {'c':[1,0,1], + 'fXm': ('c', 'fXp', [0,0,1]), 'fXp': ('p', 'fXp', [0,1]), + 'fYm': ('p', 'fYm', [1,1]), 'fYp': 'new' , + 'fZm': ('c', 'fZp', [1,0,0]), 'fZp': ('p', 'fZp', [1,0])}, + {'c':[0,1,1], + 'fXm': ('p', 'fXm', [1,1]), 'fXp': 'new' , + 'fYm': ('c', 'fYp', [0,0,1]), 'fYp': ('p', 'fYp', [0,1]), + 'fZm': ('c', 'fZp', [0,1,0]), 'fZp': ('p', 'fZp', [0,1])}, + {'c':[1,1,1], + 'fXm': ('c', 'fXp', [0,1,1]), 'fXp': ('p', 'fXp', [1,1]), + 'fYm': ('c', 'fYp', [1,0,1]), 'fYp': ('p', 'fYp', [1,1]), + 'fZm': ('c', 'fZp', [1,1,0]), 'fZp': ('p', 'fZp', [1,1])}] + + self.mesh.isNumbered = False + + self.children = np.empty((2,2,2), dtype=TreeCell) + x0, sz = self.x0, self.sz + + for face in self.faceList: + face.refine() + + def getFace(pointer): + if pointer is 'new': return None + if pointer[0] == 'p': + return self.faceDict[pointer[1]].children[pointer[2][0],pointer[2][1]] + if pointer[0] == 'c': + return self.children[pointer[2][0],pointer[2][1],pointer[2][2]].faceDict[pointer[1]] + + for O in order: + i, j, k = O['c'] + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1], x0[2] + 0.5*k*sz[2]] + fXm, fXp, fYm, fYp, fZm, fZp = getFace(O['fXm']), getFace(O['fXp']), getFace(O['fYm']), getFace(O['fYp']), getFace(O['fZm']), getFace(O['fZp']) + self.children[i,j,k] = TreeCell(self.mesh, x0=x0r, depth=self.depth+1, sz=0.5*sz, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp, fZm=fZm, fZp=fZp) + + self.mesh.cells.remove(self) + + @property + def faceIndex(self): + F = {} + for face in self.faces: + F[face] = self.faces[face].index + return F + + @property + def vol(self): return self.sz.prod() + + def viz(self, ax, color='none', text=False): + if not self.isleaf: return + x0, sz = self.x0, self.sz + ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=color, edgecolor='k')) + if text: ax.text(self.center[0],self.center[1],self.num) + + def plotGrid(self, ax, text=False): + if not self.isleaf: return + if self.dim == 2: + ax.plot(self.center[0],self.center[1],'b.') + if text: ax.text(self.center[0],self.center[1],self.num) + elif self.dim == 3: + ax.plot([self.center[0]],[self.center[1]],'b.', zs=[self.center[2]]) + if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) + + +class TreeMesh(InnerProducts, BaseMesh): + """TreeMesh""" + + _meshType = 'TREE' + + def __init__(self, h_in, x0=None): + assert type(h_in) is list, 'h_in must be a list' + h = range(len(h_in)) + for i, h_i in enumerate(h_in): + if type(h_i) in [int, long, float]: + # This gives you something over the unit cube. + h_i = np.ones(int(h_i))/int(h_i) + assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) + assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) + h[i] = h_i[:] # make a copy. + self.h = h + + if x0 is None: + x0 = np.zeros(self.dim) + else: + assert type(x0) in [list, np.ndarray], 'x0 must be a numpy array or a list' + x0 = np.array(x0, dtype=float) + assert len(x0) == self.dim, 'x0 must have the same dimensions as the mesh' + + # TODO: this has a lot of stuff which doesn't work for this style of mesh... + BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + + # set the sets for holding the cells, nodes, faces, and edges + self.cells = set() + self.nodes = set() + self.faces = set() + self.facesX = set() + self.facesY = set() + if self.dim == 3: + self.facesZ = set() + self.edges = set() + self.edgesX = set() + self.edgesY = set() + self.edgesZ = set() + + self.children = np.empty([hi.size for hi in h],dtype=TreeCell) + + if self.dim == 2: + for i in range(h[0].size): + for j in range(h[1].size): + fXm = None if i is 0 else self.children[i-1][j].fXp + fYm = None if j is 0 else self.children[i][j-1].fYp + x0i = (np.r_[x0[0], h[0][:i]]).sum() + x0j = (np.r_[x0[1], h[1][:j]]).sum() + self.children[i][j] = TreeCell(self, x0=[x0i, x0j], depth=0, sz=[h[0][i], h[1][j]], fXm=fXm, fYm=fYm) + + elif self.dim == 3: + for i in range(h[0].size): + for j in range(h[1].size): + for k in range(h[2].size): + fXm = None if i is 0 else self.children[i-1][j][k].fXp + fYm = None if j is 0 else self.children[i][j-1][k].fYp + fZm = None if k is 0 else self.children[i][j][k-1].fZp + x0i = (np.r_[x0[0], h[0][:i]]).sum() + x0j = (np.r_[x0[1], h[1][:j]]).sum() + x0k = (np.r_[x0[2], h[2][:k]]).sum() + self.children[i][j][k] = TreeCell(self, x0=[x0i, x0j, x0k], depth=0, sz=[h[0][i], h[1][j], h[2][k]], fXm=fXm, fYm=fYm, fZm=fZm) + + isNumbered = Utils.dependentProperty('_isNumbered', False, ['_faceDiv'], 'Setting this to False will delete all operators.') + + @property + def branchdepth(self): + return np.max([node.branchdepth for node in self.children.flatten('F')]) + + def refine(self, function): + for cell in self.children.flatten(): + cell.refine(function) + + def number(self): + if self.isNumbered: return + + self.sortedCells = sorted(self.cells,key=SortByX0()) + for i, sC in enumerate(self.sortedCells): sC.num = i + + self.sortedNodes = sorted(self.nodes,key=SortByX0()) + for i, sN in enumerate(self.sortedNodes): sN.num = i + + self.sortedFaceX = sorted(self.facesX,key=SortByX0()) + for i, sFx in enumerate(self.sortedFaceX): sFx.num = i + + self.sortedFaceY = sorted(self.facesY,key=SortByX0()) + for i, sFy in enumerate(self.sortedFaceY): sFy.num = i + self.nFx + + if self.dim == 3: + self.sortedFaceZ = sorted(self.facesZ,key=SortByX0()) + for i, sFz in enumerate(self.sortedFaceZ): sFz.num = i + self.nFx + self.nFy + + self.sortedEdgeX = sorted(self.edgesX,key=SortByX0()) + for i, sEx in enumerate(self.sortedEdgeX): sEx.num = i + + self.sortedEdgeY = sorted(self.edgesY,key=SortByX0()) + for i, sEy in enumerate(self.sortedEdgeY): sEy.num = i + self.nEx + + self.sortedEdgeZ = sorted(self.edgesZ,key=SortByX0()) + for i, sEz in enumerate(self.sortedEdgeZ): sEz.num = i + self.nEx + self.nEy + + self.isNumbered = True + + @property + def dim(self): return len(self.h) + + @property + def nC(self): return len(self.cells) + + @property + def nN(self): return len(self.nodes) + + @property + def nF(self): return len(self.faces) + + @property + def nFx(self): return len(self.facesX) + + @property + def nFy(self): return len(self.facesY) + + @property + def nFz(self): return None if self.dim < 3 else len(self.facesZ) + + @property + def nE(self): + if self.dim == 2: + return len(self.faces) + elif self.dim == 3: + return len(self.edges) + + @property + def nEx(self): + if self.dim == 2: + return len(self.facesY) + elif self.dim == 3: + return len(self.edgesX) + + @property + def nEy(self): + if self.dim == 2: + return len(self.facesX) + elif self.dim == 3: + return len(self.edgesY) + + @property + def nEz(self): return None if self.dim < 3 else len(self.edgesZ) + + def _grid(self, key): + self.number() + sObjs = {'CC':self.sortedCells, + 'N':self.sortedNodes, + 'Fx': self.sortedFaceX, + 'Fy': self.sortedFaceY, + 'Fz': getattr(self,'sortedFaceZ', None), + 'Ex': getattr(self,'sortedEdgeX', self.sortedFaceY), + 'Ey': getattr(self,'sortedEdgeY', self.sortedFaceX), + 'Ez': getattr(self,'sortedEdgeZ', None)}[key] + G = np.empty((len(sObjs),self.dim)) + for ii, obj in enumerate(sObjs): + G[ii,:] = obj.center + return G + + @property + def gridCC(self): + if getattr(self, '_gridCC', None) is None: + self._gridCC = self._grid('CC') + return self._gridCC + + @property + def gridN(self): + if getattr(self, '_gridN', None) is None: + self._gridN = self._grid('N') + return self._gridN + + @property + def gridFx(self): + if getattr(self, '_gridFx', None) is None: + self._gridFx = self._grid('Fx') + return self._gridFx + + @property + def gridFy(self): + if getattr(self, '_gridFy', None) is None: + self._gridFy = self._grid('Fy') + return self._gridFy + + @property + def gridFz(self): + if self.dim == 2: return None + if getattr(self, '_gridFz', None) is None: + self._gridFz = self._grid('Fz') + return self._gridFz + + @property + def gridEx(self): + if self.dim == 2: return self.gridFy + if getattr(self, '_gridEx', None) is None: + self._gridEx = self._grid('Ex') + return self._gridEx + + @property + def gridEy(self): + if self.dim == 2: return self.gridFx + if getattr(self, '_gridEy', None) is None: + self._gridEy = self._grid('Ey') + return self._gridEy + + @property + def gridEz(self): + if self.dim == 2: return None + if getattr(self, '_gridEz', None) is None: + self._gridEz = self._grid('Ez') + return self._gridEz + + @property + def vol(self): + self.number() + return np.array([cell.vol for cell in self.sortedCells]) + + @property + def area(self): + self.number() + faces = self.sortedFaceX + self.sortedFaceY + if self.dim == 3: + faces += self.sortedFaceZ + return np.array([face.area for face in faces], dtype=float) + + @property + def edge(self): + self.number() + if self.dim == 2: + edges = self.sortedFaceY + self.sortedFaceX + elif self.dim == 3: + edges = self.sortedEdgeX + self.sortedEdgeY + self.sortedEdgeZ + return np.array([e.length for e in edges], dtype=float) + + @property + def faceDiv(self): + if getattr(self, '_faceDiv', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + for cell in self.sortedCells: + faces = cell.faceDict + for face in faces: + j = faces[face].index + I += [cell.num]*len(j) + J += j + V += [-1 if 'm' in face else 1]*len(j) + VOL = self.vol + D = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + S = self.area + self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S) + return self._faceDiv + + @property + def edgeCurl(self): + """Construct the 3D curl operator.""" + assert self.dim > 2, "Edge Curl only programed for 3D." + + if getattr(self, '_edgeCurl', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + for face in self.faces: + for ii, edge in enumerate([face.edge0, face.edge1, face.edge2, face.edge3]): + j = edge.index + I += [face.num]*len(j) + J += j + isNeg = [True, False, True, False] + V += [-1 if isNeg[ii] else 1]*len(j) + C = sp.csr_matrix((V,(I,J)), shape=(self.nF, self.nE)) + S = self.area + L = self.edge + self._edgeCurl = Utils.sdiag(1/S)*C*Utils.sdiag(L) + return self._edgeCurl + + @property + def nodalGrad(self): + if getattr(self, '_nodalGrad', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + # kinda a hack for the 2D gradient + # because edges are not stored + edges = self.faces if self.dim == 2 else self.edges + for edge in edges: + if self.dim == 3: + I += [edge.num, edge.num] + elif self.dim == 2 and edge.faceType == 'x': + I += [edge.num + self.nFy, edge.num + self.nFy] + elif self.dim == 2 and edge.faceType == 'y': + I += [edge.num - self.nFx, edge.num - self.nFx] + J += [edge.node0.num, edge.node1.num] + V += [-1, 1] + G = sp.csr_matrix((V,(I,J)), shape=(self.nE, self.nN)) + L = self.edge + self._nodalGrad = Utils.sdiag(1/L)*G + return self._nodalGrad + + def _getFaceP(self, face0, face1, face2): + I, J, V = [], [], [] + for cell in self.sortedCells: + face = cell.faceDict[face0] + if face.isleaf: + j = face.index + elif self.dim == 2: + j = face.children[0 if 'm' in face1 else 1].index + elif self.dim == 3: + j = face.children[0 if 'm' in face1 else 1, + 0 if 'm' in face2 else 1].index + lenj = len(j) + I += [cell.num]*lenj + J += j + V += [1./lenj]*lenj + return sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + + def _getEdgeP(self, edge0, edge1, edge2): + I, J, V = [], [], [] + for cell in self.sortedCells: + if self.dim == 2: + e2f = lambda e: ('f' + {'X':'Y','Y':'X'}[e[1]] + + {'0':'m','1':'p'}[e[2]]) + face = cell.faceDict[e2f(edge0)] + if face.isleaf: + j = face.index + else: + j = face.children[0 if 'm' in e2f(edge1) else 1].index + # Need to flip the numbering for edges + if 'X' in edge0: + j = [jj - self.nFx for jj in j] + elif 'Y' in edge0: + j = [jj + self.nFy for jj in j] + elif self.dim == 3: + edge = cell.edgeDict[edge0] + if edge.isleaf: + j = edge.index + else: + mSide = lambda e: {'0':True,'1':True,'2':False,'3':False}[e[2]] + j = edge.children[0 if mSide(edge1) else 1, + 0 if mSide(edge2) else 1].index + lenj = len(j) + I += [cell.num]*lenj + J += j + V += [1./lenj]*lenj + return sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nE)) + + def _getFacePxx(self, xFace, yFace): + self.number() + xP = self._getFaceP(xFace, yFace, None) + yP = self._getFaceP(yFace, xFace, None) + return sp.vstack((xP, yP)) + + def _getEdgePxx(self, xEdge, yEdge): + self.number() + xP = self._getEdgeP(xEdge, yEdge, None) + yP = self._getEdgeP(yEdge, xEdge, None) + return sp.vstack((xP, yP)) + + def _getFacePxxx(self, xFace, yFace, zFace): + self.number() + xP = self._getFaceP(xFace, yFace, zFace) + yP = self._getFaceP(yFace, xFace, zFace) + zP = self._getFaceP(zFace, xFace, yFace) + return sp.vstack((xP, yP, zP)) + + def _getEdgePxxx(self, xEdge, yEdge, zEdge): + self.number() + xP = self._getEdgeP(xEdge, yEdge, zEdge) + yP = self._getEdgeP(yEdge, xEdge, zEdge) + zP = self._getEdgeP(zEdge, xEdge, yEdge) + return sp.vstack((xP, yP, zP)) + + def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, plotE=False, plotEx=False, plotEy=False, plotEz=False, showIt=False): + axOpts = {'projection':'3d'} if self.dim == 3 else {} + if ax is None: ax = plt.subplot(111, **axOpts) + + if plotC: [c.plotGrid(ax, text=text) for c in self.cells] + if plotF: [f.plotGrid(ax, text=text) for f in self.faces] + if plotE and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edges] + if plotEx and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesX] + if plotEy and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesY] + if plotEz and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesZ] + + ax.set_xlim((self.x0[0], self.h[0].sum())) + ax.set_ylim((self.x0[1], self.h[1].sum())) + if self.dim == 3: + ax.set_zlim((self.x0[2], self.h[2].sum())) + if showIt: plt.show() + + def plotImage(self, I, ax=None, showIt=True): + if self.dim == 2: + self._plotImage2D(I, ax=ax, showIt=showIt) + elif self.dim == 3: + raise NotImplementedError('3D visualization is not yet implemented.') + + def _plotImage2D(self, I, ax=None, showIt=True): + if ax is None: ax = plt.subplot(111) + jet = cm = plt.get_cmap('jet') + cNorm = colors.Normalize(vmin=I.min(), vmax=I.max()) + scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet) + ax.set_xlim((self.x0[0], self.h[0].sum())) + ax.set_ylim((self.x0[1], self.h[1].sum())) + for ii, node in enumerate(self.sortedCells): + node.viz(ax=ax, color=scalarMap.to_rgba(I[ii])) + scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots + plt.colorbar(scalarMap) + if showIt: plt.show() + + + +if __name__ == '__main__': + M = TreeMesh([np.ones(x) for x in [4,10]]) + + def function(xc): + r = xc - np.r_[2.,6.] + dist = np.sqrt(r.dot(r)) + if dist < 1.0: + return 3 + if dist < 1.5: + return 2 + else: + return 1 + + M.refine(function) + + DIV = M.faceDiv + Mf = M.getFaceInnerProduct() + # plt.subplot(211) + # plt.spy(DIV) + M.plotGrid(ax=plt.subplot(111),text=True,showIt=True) + + q = np.zeros(M.nC) + q[208] = -1.0 + q[291] = 1.0 + b = Solver(-DIV*Mf*DIV.T).solve(q) + plt.figure() + M.plotImage(b) + # plt.gca().invert_yaxis() + print M.vol + plt.show() diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 3b8e1eef..3da22a01 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -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 diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 33cfeef8..6941607d 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -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)) diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index ea2fe831..7d08b9fc 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -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) diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 4e5963f6..8e08250b 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -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 diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py index 2a07aa69..9b12d8d1 100644 --- a/SimPEG/Parameters.py +++ b/SimPEG/Parameters.py @@ -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)) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 386fac9d..314af81f 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -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. diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index aa62449a..b360d726 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -10,7 +10,7 @@ class BaseRegularization(object): """ - __metaclass__ = Utils.Save.Savable + __metaclass__ = Utils.SimPEGMetaClass modelPair = Model.BaseModel #: Some regularizations only work on specific models diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 3e9c2d02..60ca5558 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -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) diff --git a/SimPEG/Tests/test_TreeMesh.py b/SimPEG/Tests/test_TreeMesh.py new file mode 100644 index 00000000..f8120859 --- /dev/null +++ b/SimPEG/Tests/test_TreeMesh.py @@ -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() diff --git a/SimPEG/Tests/test_forward_DCproblem.py b/SimPEG/Tests/test_forward_DCproblem.py deleted file mode 100644 index 1a44b350..00000000 --- a/SimPEG/Tests/test_forward_DCproblem.py +++ /dev/null @@ -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() diff --git a/SimPEG/Tests/test_interpolation.py b/SimPEG/Tests/test_interpolation.py index d01dcbc3..018f9091 100644 --- a/SimPEG/Tests/test_interpolation.py +++ b/SimPEG/Tests/test_interpolation.py @@ -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) diff --git a/SimPEG/Tests/test_model.py b/SimPEG/Tests/test_model.py index 0c3b2c54..e6a8844b 100644 --- a/SimPEG/Tests/test_model.py +++ b/SimPEG/Tests/test_model.py @@ -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() diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index d1c38443..67741aa7 100644 --- a/SimPEG/Tests/test_operators.py +++ b/SimPEG/Tests/test_operators.py @@ -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)) diff --git a/SimPEG/Tests/test_tensorMesh.py b/SimPEG/Tests/test_tensorMesh.py index 3e01181b..45c03399 100644 --- a/SimPEG/Tests/test_tensorMesh.py +++ b/SimPEG/Tests/test_tensorMesh.py @@ -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 diff --git a/SimPEG/Utils/Geophysics/__init__.py b/SimPEG/Utils/Geophysics/__init__.py deleted file mode 100644 index bcf01943..00000000 --- a/SimPEG/Utils/Geophysics/__init__.py +++ /dev/null @@ -1 +0,0 @@ -import emSources \ No newline at end of file diff --git a/SimPEG/Utils/Geophysics/emSources/__init__.py b/SimPEG/Utils/Geophysics/emSources/__init__.py deleted file mode 100644 index c05fda69..00000000 --- a/SimPEG/Utils/Geophysics/emSources/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from emSources import MagneticDipoleVectorPotential \ No newline at end of file diff --git a/SimPEG/Utils/Geophysics/emSources/emSources.py b/SimPEG/Utils/Geophysics/emSources/emSources.py deleted file mode 100644 index db492a94..00000000 --- a/SimPEG/Utils/Geophysics/emSources/emSources.py +++ /dev/null @@ -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. ' - - :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 \ No newline at end of file diff --git a/SimPEG/Utils/Save.py b/SimPEG/Utils/Save.py deleted file mode 100644 index a9c77191..00000000 --- a/SimPEG/Utils/Save.py +++ /dev/null @@ -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:: - - - - - - - - - - - - - """ - 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) - diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 1daa2e29..b2d0d56a 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -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) diff --git a/SimPEG/Utils/sputils.py b/SimPEG/Utils/sputils.py index 06eef80d..f671be43 100644 --- a/SimPEG/Utils/sputils.py +++ b/SimPEG/Utils/sputils.py @@ -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""" diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 4967608e..17f64d54 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -11,7 +11,6 @@ import ObjFunction import Optimization import Inversion import Parameters -import Examples import Tests diff --git a/SimPEG/visualize/__init__.py b/SimPEG/visualize/__init__.py deleted file mode 100644 index 485f9782..00000000 --- a/SimPEG/visualize/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -import vtk -#import mpl diff --git a/SimPEG/visualize/vtk/__init__.py b/SimPEG/visualize/vtk/__init__.py deleted file mode 100644 index 6d60ed5c..00000000 --- a/SimPEG/visualize/vtk/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -from vtkTools import vtkTools -from vtkView import vtkView \ No newline at end of file diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py deleted file mode 100644 index 68326662..00000000 --- a/SimPEG/visualize/vtk/vtkTools.py +++ /dev/null @@ -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() - diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py deleted file mode 100644 index 3b52b992..00000000 --- a/SimPEG/visualize/vtk/vtkView.py +++ /dev/null @@ -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() \ No newline at end of file diff --git a/SimPEG/Examples/Linear.py b/Tutorials/Linear.py similarity index 97% rename from SimPEG/Examples/Linear.py rename to Tutorials/Linear.py index a83b26f1..e18750c1 100644 --- a/SimPEG/Examples/Linear.py +++ b/Tutorials/Linear.py @@ -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):