diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 87ec7094..034cdd01 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -10,7 +10,8 @@ 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, materialProperty=None, returnP=False, invertProperty=False): + def getFaceInnerProduct(self, materialProperty=None, returnP=False, + invertProperty=False): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param bool returnP: returns the projection matrices @@ -18,6 +19,14 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ + fast = None + + if returnP is False and hasattr(self, '_fastFaceInnerProduct'): + fast = self._fastFaceInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) + + if fast is not None: + return fast + if invertProperty: materialProperty = invPropertyTensor(self, materialProperty) @@ -62,6 +71,34 @@ class InnerProducts(object): else: return A + def getFaceInnerProductDeriv(self, materialProperty=None, P=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + + fast = None + + if hasattr(self, '_fastFaceInnerProductDeriv'): + fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty) + + if fast is not None: + return fast + + raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') + + if P is None: + M, P = getFaceInnerProduct(self, materialProperty=materialProperty, returnP=True) + + d = self.dim + + if d == 1: + P[0].T * sp.identity(n) * P[0] + + + + def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 8a80ffea..f55c6df6 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -39,6 +39,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): def __init__(self, h_in, x0=None): assert type(h_in) is list, 'h_in must be a list' + assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3' h = range(len(h_in)) for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: @@ -491,6 +492,49 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) return indxd, indxu, indyd, indyu, indzd, indzu + def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None: + materialProperty = np.ones(self.nC) + + if materialProperty.size == self.nC: + if invertProperty: + materialProperty = 1./materialProperty + Av = self.aveF2CC + V = Utils.sdiag(self.vol) + return self.dim * Utils.sdiag(Av.T * V * materialProperty) + if materialProperty.size == self.nC*self.dim: + if invertProperty: + materialProperty = 1./materialProperty + Av = self.aveF2CCV + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) + + + + def _fastFaceInnerProductDeriv(self, materialProperty=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None or materialProperty.size == self.nC: + Av = self.aveF2CC + return self.dim * Av.T * Utils.sdiag(self.vol) + if materialProperty.size == self.nC*self.dim: # anisotropic + Av = self.aveF2CCV + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + return Av.T * V + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 0db8c007..200d2884 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -30,7 +30,7 @@ class TestInnerProducts(OrderTest): sigma = np.c_[call(sigma1, Gc)] analytic = 647./360 # Found using sympy. elif self.sigmaTest == 3: - sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] + sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] analytic = 37./12 # Found using sympy. elif self.sigmaTest == 6: sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py new file mode 100644 index 00000000..f5484568 --- /dev/null +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -0,0 +1,40 @@ +import numpy as np +import unittest +from SimPEG import * +from TestUtils import checkDerivative + + +class TestInnerProductsDerivs(unittest.TestCase): + def setUp(self): + pass + + def test_FaceIP_derivs_isotropic(self): + for d in range(3): + mesh = Mesh.TensorMesh([10,5,4][d:]) + M,Ps = mesh.getFaceInnerProduct(returnP=True) + v = np.random.rand(mesh.nF) + def fun(sig): + M = mesh.getFaceInnerProduct(sig) + Md = mesh.getFaceInnerProductDeriv(sig) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC) + passed = checkDerivative(fun, sig, plotIt=False) + self.assertTrue(passed) + + def test_FaceIP_derivs_anisotropic(self): + for d in range(3): + mesh = Mesh.TensorMesh([10,5,4][d:]) + M,Ps = mesh.getFaceInnerProduct(returnP=True) + v = np.random.rand(mesh.nF) + def fun(sig): + M = mesh.getFaceInnerProduct(sig) + Md = mesh.getFaceInnerProductDeriv(sig) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC*mesh.dim) + passed = checkDerivative(fun, sig, plotIt=False) + self.assertTrue(passed) + + + +if __name__ == '__main__': + unittest.main() diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index dec9c4c8..06c0b2ed 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -131,6 +131,7 @@ If ``returnP=True`` is requested in any of these methods the projection matrices The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same. + Defining Tensor Properties -------------------------- @@ -158,6 +159,42 @@ Depending on the number of columns (either 1, 3, or 6) of mu, the material prope \vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{3} \\ \mu_{3} & \mu_{2} \end{matrix}\right] +Structure of Matrices +--------------------- + +Both the isotropic, and anisotropic material properties result in a diagonal mass matrix. +Which is nice and easy to invert if necessary, however, in the fully anisotropic case which is not aligned with the grid, the matrix is not diagonal. This can be seen for a 3D mesh in the figure below. + +.. plot:: + + from SimPEG import * + mesh = Mesh.TensorMesh([10,50,3]) + m1 = np.random.rand(mesh.nC) + m2 = np.random.rand(mesh.nC,3) + m3 = np.random.rand(mesh.nC,6) + M = range(3) + M[0] = mesh.getFaceInnerProduct(m1) + M[1] = mesh.getFaceInnerProduct(m2) + M[2] = mesh.getFaceInnerProduct(m3) + plt.figure(figsize=(13,5)) + for i, lab in enumerate(['Isotropic','Anisotropic','Tensor']): + plt.subplot(131 + i) + plt.spy(M[i],ms=0.5,color='k') + plt.tick_params(axis='both',which='both',labeltop='off',labelleft='off') + plt.title(lab + ' Material Property') + plt.show() + + +Taking Derivatives +------------------ + +TODO: Take the derivatives of the tensors. + +.. math:: + + \left[\begin{smallmatrix}0.5 \sigma_{1} & 0 & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0 & 0.5 \sigma_{1} & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0.5 \sigma_{2} & 0\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0 & 0.5 \sigma_{2}\end{smallmatrix}\right] + + The API -------