```from __future__ import print_function
from SimPEG import Mesh, Utils, Tests
import numpy as np
import sympy
from sympy.abc import r, t, z
import unittest
import scipy.sparse as sp

TOL = 1e-1
TOLD = 0.7  # tolerance on deriv checks

class FaceInnerProductFctsIsotropic(object):
""" Some made up face functions to test the face inner product """
def fcts(self):
j = sympy.Matrix([
r**2 * z,
r * z**2
])

# Create an isotropic sigma vector
Sig = sympy.Matrix([
[420/(sympy.pi)*(r*z)**2, 0],
[0, 420/(sympy.pi)*(r*z)**2],
])

return j, Sig

def sol(self):
# Do the inner product! - we are in cyl coordinates!
j, Sig = self.fcts()
jTSj = j.T*Sig*j
# we are integrating in cyl coordinates
ans  = sympy.integrate(sympy.integrate(sympy.integrate(r * jTSj,
(r, 0, 1)),
(t, 0, 2*sympy.pi)),
(z, 0, 1))[0] # The `[0]` is to make it an int.

return ans

def vectors(self, mesh):
""" Get Vectors sig, sr. jx from sympy"""
j, Sig = self.fcts()

f_jr = sympy.lambdify((r, z), j[0], 'numpy')
f_jz = sympy.lambdify((r, z), j[1], 'numpy')
f_sigr = sympy.lambdify((r, z), Sig[0], 'numpy')
# f_sigz = sympy.lambdify((r,z), Sig[1], 'numpy')

jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 2])
jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 2])
sigr = f_sigr(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

return sigr, np.r_[jr, jz]

class FaceInnerProductFunctionsDiagAnisotropic(FaceInnerProductFctsIsotropic):
"""
Some made up face functions to test the diagonally anisotropic face
inner product
"""

def fcts(self):
j = sympy.Matrix([
r**2 * z,
r * z**2
])

# Create an isotropic sigma vector
Sig = sympy.Matrix([
[120/(sympy.pi)*(r*z)**2, 0],
[0, 420/(sympy.pi)*(r*z)**2],
])

return j, Sig

def vectors(self, mesh):
""" Get Vectors sig, sr. jx from sympy"""
j, Sig = self.fcts()

f_jr = sympy.lambdify((r, z), j[0], 'numpy')
f_jz = sympy.lambdify((r, z), j[1], 'numpy')
f_sigr = sympy.lambdify((r, z), Sig[0], 'numpy')
f_sigz = sympy.lambdify((r, z), Sig[3], 'numpy')

jr = f_jr(mesh.gridFx[:, 0], mesh.gridFx[:, 2])
jz = f_jz(mesh.gridFz[:, 0], mesh.gridFz[:, 2])
sigr = f_sigr(mesh.gridCC[:, 0], mesh.gridCC[:, 2])
sigz = f_sigz(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

return np.c_[sigr, sigr, sigz], np.r_[jr, jz]

class EdgeInnerProductFctsIsotropic(object):
""" Some made up edge functions to test the edge inner product """

def fcts(self):
h = sympy.Matrix([r**2 * z])

# Create an isotropic sigma vector
Sig = sympy.Matrix([200/(sympy.pi)*(r*z)**2])

return h, Sig

def sol(self):
h, Sig = self.fcts()
# Do the inner product! - we are in cyl coordinates!
hTSh = h.T*Sig*h
ans  = sympy.integrate(sympy.integrate(sympy.integrate(r * hTSh,
(r, 0, 1)),
(t, 0, 2*sympy.pi)),
(z, 0, 1))[0] # The `[0]` is to make it an int.
return ans

def vectors(self, mesh):
""" Get Vectors sig, sr. jx from sympy"""
h, Sig = self.fcts()

f_h = sympy.lambdify((r, z), h[0], 'numpy')
f_sig = sympy.lambdify((r, z), Sig[0], 'numpy')

ht = f_h(mesh.gridEy[:, 0], mesh.gridEy[:, 2])
sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

return sig, np.r_[ht]

class EdgeInnerProductFunctionsDiagAnisotropic(EdgeInnerProductFctsIsotropic):
"""
Some made up edge functions to test the diagonally anisotropic edge
inner product
"""

def vectors(self, mesh):
h, Sig = self.fcts()

f_h = sympy.lambdify((r, z), h[0], 'numpy')
f_sig = sympy.lambdify((r, z), Sig[0], 'numpy')

ht = f_h(mesh.gridEy[:, 0], mesh.gridEy[:, 2])
sig = f_sig(mesh.gridCC[:, 0], mesh.gridCC[:, 2])

return np.c_[sig, sig, sig], np.r_[ht]

class TestCylInnerProducts_simple(unittest.TestCase):

def setUp(self):
n = 100.
self.mesh = Mesh.CylMesh([n, 1, n])

def test_FaceInnerProductIsotropic(self):
# Here we will make up some j vectors that vary in space
# j = [j_r, j_z] - to test face inner products

fcts = FaceInnerProductFctsIsotropic()
sig, jv = fcts.vectors(self.mesh)
MfSig = self.mesh.getFaceInnerProduct(sig)
numeric_ans = jv.T.dot(MfSig.dot(jv))

ans = fcts.sol()

print('------ Testing Face Inner Product-----------')
print(' Analytic: {analytic}, Numeric: {numeric}, '
'ratio (num/ana): {ratio}'.format(
analytic=ans, numeric=numeric_ans,
ratio=float(numeric_ans)/ans))
assert(np.abs(ans-numeric_ans) < TOL)

def test_FaceInnerProductDiagAnisotropic(self):
# Here we will make up some j vectors that vary in space
# j = [j_r, j_z] - to test face inner products

fcts = FaceInnerProductFunctionsDiagAnisotropic()
sig, jv = fcts.vectors(self.mesh)
MfSig = self.mesh.getFaceInnerProduct(sig)
numeric_ans = jv.T.dot(MfSig.dot(jv))

ans = fcts.sol()

print('------ Testing Face Inner Product Anisotropic -----------')
print(' Analytic: {analytic}, Numeric: {numeric}, '
'ratio (num/ana): {ratio}'.format(
analytic=ans, numeric=numeric_ans,
ratio=float(numeric_ans)/ans))
assert(np.abs(ans-numeric_ans) < TOL)

def test_EdgeInnerProduct(self):
# Here we will make up some j vectors that vary in space
# h = [h_t] - to test edge inner products

fcts = EdgeInnerProductFctsIsotropic()
sig, hv = fcts.vectors(self.mesh)
MeSig = self.mesh.getEdgeInnerProduct(sig)
numeric_ans = hv.T.dot(MeSig.dot(hv))

ans = fcts.sol()

print('------ Testing Edge Inner Product-----------')
print(' Analytic: {analytic}, Numeric: {numeric}, '
'ratio (num/ana): {ratio}'.format(
analytic=ans, numeric=numeric_ans,
ratio=float(numeric_ans)/ans))
assert(np.abs(ans-numeric_ans) < TOL)

def test_EdgeInnerProductDiagAnisotropic(self):
# Here we will make up some j vectors that vary in space
# h = [h_t] - to test edge inner products

fcts = EdgeInnerProductFunctionsDiagAnisotropic()

sig, hv = fcts.vectors(self.mesh)
MeSig = self.mesh.getEdgeInnerProduct(sig)
numeric_ans = hv.T.dot(MeSig.dot(hv))

ans = fcts.sol()

print('------ Testing Edge Inner Product Anisotropic -----------')
print(' Analytic: {analytic}, Numeric: {numeric}, '
'ratio (num/ana): {ratio}'.format(
analytic=ans, numeric=numeric_ans,
ratio=float(numeric_ans)/ans))
assert(np.abs(ans-numeric_ans) < TOL)

class TestCylFaceInnerProducts_Order(Tests.OrderTest):

meshTypes = ['uniformCylMesh']
meshDimension = 2

def getError(self):
fct = FaceInnerProductFctsIsotropic()
sig, jv = fct.vectors(self.M)
Msig = self.M.getFaceInnerProduct(sig)
return float(fct.sol()) - jv.T.dot(Msig.dot(jv))

def test_order(self):
self.orderTest()

class TestCylEdgeInnerProducts_Order(Tests.OrderTest):

meshTypes = ['uniformCylMesh']
meshDimension = 2

def getError(self):
fct = EdgeInnerProductFctsIsotropic()
sig, ht = fct.vectors(self.M)
Msig = self.M.getEdgeInnerProduct(sig)
return float(fct.sol()) - ht.T.dot(Msig.dot(ht))

def test_order(self):
self.orderTest()

class TestCylFaceInnerProductsDiagAnisotropic_Order(Tests.OrderTest):

meshTypes = ['uniformCylMesh']
meshDimension = 2

def getError(self):
fct = FaceInnerProductFunctionsDiagAnisotropic()
sig, jv = fct.vectors(self.M)
Msig = self.M.getFaceInnerProduct(sig)
return float(fct.sol()) - jv.T.dot(Msig.dot(jv))

def test_order(self):
self.orderTest()

class TestCylEdgeInnerProducts_Order(Tests.OrderTest):

meshTypes = ['uniformCylMesh']
meshDimension = 2

def getError(self):
fct = EdgeInnerProductFunctionsDiagAnisotropic()
sig, ht = fct.vectors(self.M)
Msig = self.M.getEdgeInnerProduct(sig)
return float(fct.sol()) - ht.T.dot(Msig.dot(ht))

def test_order(self):
self.orderTest()

class TestCylInnerProducts_Deriv(unittest.TestCase):

def setUp(self):
n = 2
self.mesh = Mesh.CylMesh([n, 1, n])
self.face_vec = np.random.rand(self.mesh.nF)
self.edge_vec = np.random.rand(self.mesh.nE)
# make up a smooth function
self.x0 = 2*self.mesh.gridCC[:, 0]**2 + self.mesh.gridCC[:, 2]**4

def test_FaceInnerProductIsotropicDeriv(self):

def fun(x):
MfSig = self.mesh.getFaceInnerProduct(x)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(self.x0)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec)

print('Testing FaceInnerProduct Isotropic')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD, plotIt=False))

def test_FaceInnerProductIsotropicDerivInvProp(self):

def fun(x):
MfSig = self.mesh.getFaceInnerProduct(x, invProp=True)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(self.x0,
invProp=True)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec)

print('Testing FaceInnerProduct Isotropic InvProp')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_FaceInnerProductIsotropicDerivInvMat(self):

def fun(x):
MfSig = self.mesh.getFaceInnerProduct(x, invMat=True)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(self.x0,
invMat=True)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec)

print('Testing FaceInnerProduct Isotropic InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_FaceInnerProductIsotropicDerivInvPropInvMat(self):

def fun(x):
MfSig = self.mesh.getFaceInnerProduct(x, invProp=True, invMat=True)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(self.x0,
invProp=True,
invMat=True)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec)

print('Testing FaceInnerProduct Isotropic InvProp InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductIsotropicDeriv(self):

def fun(x):
MeSig = self.mesh.getEdgeInnerProduct(x)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(self.x0)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec)

print('Testing EdgeInnerProduct Isotropic')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductIsotropicDerivInvProp(self):

def fun(x):
MeSig = self.mesh.getEdgeInnerProduct(x, invProp=True)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(self.x0,
invProp=True)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec)

print('Testing EdgeInnerProduct Isotropic InvProp')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductIsotropicDerivInvMat(self):

def fun(x):
MeSig = self.mesh.getEdgeInnerProduct(x, invMat=True)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(self.x0,
invMat=True)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec)

print('Testing EdgeInnerProduct Isotropic InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductIsotropicDerivInvPropInvMat(self):

def fun(x):
MeSig = self.mesh.getEdgeInnerProduct(x, invProp=True, invMat=True)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(self.x0,
invProp=True,
invMat=True)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec)

print('Testing EdgeInnerProduct Isotropic InvProp InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

class TestCylInnerProductsAnisotropic_Deriv(unittest.TestCase):

def setUp(self):
n = 60
self.mesh = Mesh.CylMesh([n, 1, n])
self.face_vec = np.random.rand(self.mesh.nF)
self.edge_vec = np.random.rand(self.mesh.nE)
# make up a smooth function
self.x0 = np.array([2*self.mesh.gridCC[:, 0]**2 + self.mesh.gridCC[:, 2]**4
])

def test_FaceInnerProductAnisotropicDeriv(self):

def fun(x):
# fake anisotropy (testing anistropic implementation with isotropic
# vector). First order behavior expected for fully anisotropic
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Eye, Zero, Eye])])

MfSig = self.mesh.getFaceInnerProduct(x)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(x0)
return MfSig*self.face_vec ,  MfSigDeriv(self.face_vec) * P.T

print('Testing FaceInnerProduct Anisotropic')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD, plotIt=False))

def test_FaceInnerProductAnisotropicDerivInvProp(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Eye, Zero, Eye])])

MfSig = self.mesh.getFaceInnerProduct(x, invProp=True)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(x0,
invProp=True)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec) * P.T

print('Testing FaceInnerProduct Anisotropic InvProp')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_FaceInnerProductAnisotropicDerivInvMat(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Eye, Zero, Eye])])

MfSig = self.mesh.getFaceInnerProduct(x, invMat=True)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(x0, invMat=True)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec) * P.T

print('Testing FaceInnerProduct Anisotropic InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_FaceInnerProductAnisotropicDerivInvPropInvMat(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Eye, Zero, Eye])])

MfSig = self.mesh.getFaceInnerProduct(x, invProp=True, invMat=True)
MfSigDeriv = self.mesh.getFaceInnerProductDeriv(x0,
invProp=True,
invMat=True)
return MfSig*self.face_vec, MfSigDeriv(self.face_vec) * P.T

print('Testing FaceInnerProduct Anisotropic InvProp InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductAnisotropicDeriv(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Zero, Eye, Zero])])

MeSig = self.mesh.getEdgeInnerProduct(x.reshape(self.mesh.nC, 3))
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(x0)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec) * P.T

print('Testing EdgeInnerProduct Anisotropic')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductAnisotropicDerivInvProp(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Zero, Eye, Zero])])

MeSig = self.mesh.getEdgeInnerProduct(x, invProp=True)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(x0, invProp=True)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec) * P.T

print('Testing EdgeInnerProduct Anisotropic InvProp')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductAnisotropicDerivInvMat(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Zero, Eye, Zero])])

MeSig = self.mesh.getEdgeInnerProduct(x, invMat=True)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(x0, invMat=True)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec) * P.T

print('Testing EdgeInnerProduct Anisotropic InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

def test_EdgeInnerProductAnisotropicDerivInvPropInvMat(self):

def fun(x):
x = np.repeat(np.atleast_2d(x), 3, axis=0).T
x0 = np.repeat(self.x0, 3, axis=0).T

Zero = sp.csr_matrix((self.mesh.nC, self.mesh.nC))
Eye = sp.eye(self.mesh.nC)
P = sp.vstack([sp.hstack([Zero, Eye, Zero])])

MeSig = self.mesh.getEdgeInnerProduct(x, invProp=True, invMat=True)
MeSigDeriv = self.mesh.getEdgeInnerProductDeriv(x0,
invProp=True,
invMat=True)
return MeSig*self.edge_vec, MeSigDeriv(self.edge_vec) * P.T

print('Testing EdgeInnerProduct Anisotropic InvProp InvMat')
return self.assertTrue(Tests.checkDerivative(fun, self.x0, num=7,
tolerance=TOLD,
plotIt=False))

if __name__ == '__main__':
unittest.main()
```