mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 02:47:11 +08:00
Rowan Cockett
403 lines
13 KiB
Python
403 lines
13 KiB
Python
import numpy as np
|
|
import unittest
|
|
from SimPEG import Utils, Tests
|
|
|
|
|
|
class TestInnerProducts(Tests.OrderTest):
|
|
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
|
|
|
meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
|
|
meshDimension = 3
|
|
meshSizes = [16, 32]
|
|
|
|
def getError(self):
|
|
|
|
call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
|
|
|
ex = lambda x, y, z: x**2+y*z
|
|
ey = lambda x, y, z: (z**2)*x+y*z
|
|
ez = lambda x, y, z: y**2+x*z
|
|
|
|
sigma1 = lambda x, y, z: x*y+1
|
|
sigma2 = lambda x, y, z: x*z+2
|
|
sigma3 = lambda x, y, z: 3+z*y
|
|
sigma4 = lambda x, y, z: 0.1*x*y*z
|
|
sigma5 = lambda x, y, z: 0.2*x*y
|
|
sigma6 = lambda x, y, z: 0.1*z
|
|
|
|
Gc = self.M.gridCC
|
|
if self.sigmaTest == 1:
|
|
sigma = np.c_[call(sigma1, Gc)]
|
|
analytic = 647./360 # Found using sympy.
|
|
elif self.sigmaTest == 3:
|
|
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),
|
|
call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)]
|
|
analytic = 69881./21600 # Found using sympy.
|
|
|
|
if self.location == 'edges':
|
|
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
|
Ec = np.vstack((cart(self.M.gridEx),
|
|
cart(self.M.gridEy),
|
|
cart(self.M.gridEz)))
|
|
E = self.M.projectEdgeVector(Ec)
|
|
|
|
if self.invProp:
|
|
A = self.M.getEdgeInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
|
else:
|
|
A = self.M.getEdgeInnerProduct(sigma)
|
|
numeric = E.T.dot(A.dot(E))
|
|
elif self.location == 'faces':
|
|
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
|
Fc = np.vstack((cart(self.M.gridFx),
|
|
cart(self.M.gridFy),
|
|
cart(self.M.gridFz)))
|
|
F = self.M.projectFaceVector(Fc)
|
|
|
|
if self.invProp:
|
|
A = self.M.getFaceInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
|
else:
|
|
A = self.M.getFaceInnerProduct(sigma)
|
|
numeric = F.T.dot(A.dot(F))
|
|
|
|
err = np.abs(numeric - analytic)
|
|
return err
|
|
|
|
def test_order1_edges(self):
|
|
self.name = "Edge Inner Product - Isotropic"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 1
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order1_edges_invProp(self):
|
|
self.name = "Edge Inner Product - Isotropic - invProp"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 1
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order3_edges(self):
|
|
self.name = "Edge Inner Product - Anisotropic"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 3
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order3_edges_invProp(self):
|
|
self.name = "Edge Inner Product - Anisotropic - invProp"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 3
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order6_edges(self):
|
|
self.name = "Edge Inner Product - Full Tensor"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 6
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order6_edges_invProp(self):
|
|
self.name = "Edge Inner Product - Full Tensor - invProp"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 6
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order1_faces(self):
|
|
self.name = "Face Inner Product - Isotropic"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 1
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order1_faces_invProp(self):
|
|
self.name = "Face Inner Product - Isotropic - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 1
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order3_faces(self):
|
|
self.name = "Face Inner Product - Anisotropic"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 3
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order3_faces_invProp(self):
|
|
self.name = "Face Inner Product - Anisotropic - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 3
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order6_faces(self):
|
|
self.name = "Face Inner Product - Full Tensor"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 6
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order6_faces_invProp(self):
|
|
self.name = "Face Inner Product - Full Tensor - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 6
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
|
|
class TestInnerProducts2D(Tests.OrderTest):
|
|
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
|
|
|
meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
|
|
meshDimension = 2
|
|
meshSizes = [4, 8, 16, 32, 64, 128]
|
|
|
|
def getError(self):
|
|
|
|
z = 5 # Because 5 is just such a great number.
|
|
|
|
call = lambda fun, xy: fun(xy[:, 0], xy[:, 1])
|
|
|
|
ex = lambda x, y: x**2+y*z
|
|
ey = lambda x, y: (z**2)*x+y*z
|
|
|
|
sigma1 = lambda x, y: x*y+1
|
|
sigma2 = lambda x, y: x*z+2
|
|
sigma3 = lambda x, y: 3+z*y
|
|
|
|
Gc = self.M.gridCC
|
|
if self.sigmaTest == 1:
|
|
sigma = np.c_[call(sigma1, Gc)]
|
|
analytic = 144877./360 # Found using sympy. z=5
|
|
elif self.sigmaTest == 2:
|
|
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
|
analytic = 189959./120 # Found using sympy. z=5
|
|
elif self.sigmaTest == 3:
|
|
sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
|
analytic = 781427./360 # Found using sympy. z=5
|
|
|
|
if self.location == 'edges':
|
|
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
|
Ec = np.vstack((cart(self.M.gridEx),
|
|
cart(self.M.gridEy)))
|
|
E = self.M.projectEdgeVector(Ec)
|
|
if self.invProp:
|
|
A = self.M.getEdgeInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
|
else:
|
|
A = self.M.getEdgeInnerProduct(sigma)
|
|
numeric = E.T.dot(A.dot(E))
|
|
elif self.location == 'faces':
|
|
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
|
Fc = np.vstack((cart(self.M.gridFx),
|
|
cart(self.M.gridFy)))
|
|
F = self.M.projectFaceVector(Fc)
|
|
|
|
if self.invProp:
|
|
A = self.M.getFaceInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
|
else:
|
|
A = self.M.getFaceInnerProduct(sigma)
|
|
numeric = F.T.dot(A.dot(F))
|
|
|
|
err = np.abs(numeric - analytic)
|
|
return err
|
|
|
|
def test_order1_edges(self):
|
|
self.name = "2D Edge Inner Product - Isotropic"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 1
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order1_edges_invProp(self):
|
|
self.name = "2D Edge Inner Product - Isotropic - invProp"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 1
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order3_edges(self):
|
|
self.name = "2D Edge Inner Product - Anisotropic"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 2
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order3_edges_invProp(self):
|
|
self.name = "2D Edge Inner Product - Anisotropic - invProp"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 2
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order6_edges(self):
|
|
self.name = "2D Edge Inner Product - Full Tensor"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 3
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order6_edges_invProp(self):
|
|
self.name = "2D Edge Inner Product - Full Tensor - invProp"
|
|
self.location = 'edges'
|
|
self.sigmaTest = 3
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order1_faces(self):
|
|
self.name = "2D Face Inner Product - Isotropic"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 1
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order1_faces_invProp(self):
|
|
self.name = "2D Face Inner Product - Isotropic - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 1
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order2_faces(self):
|
|
self.name = "2D Face Inner Product - Anisotropic"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 2
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order2_faces_invProp(self):
|
|
self.name = "2D Face Inner Product - Anisotropic - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 2
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
def test_order3_faces(self):
|
|
self.name = "2D Face Inner Product - Full Tensor"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 3
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order3_faces_invProp(self):
|
|
self.name = "2D Face Inner Product - Full Tensor - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 3
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
|
|
|
|
class TestInnerProducts1D(Tests.OrderTest):
|
|
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
|
|
|
meshTypes = ['uniformTensorMesh']
|
|
meshDimension = 1
|
|
meshSizes = [4, 8, 16, 32, 64, 128]
|
|
|
|
def getError(self):
|
|
|
|
y = 12 # Because 12 is just such a great number.
|
|
z = 5 # Because 5 is just such a great number as well!
|
|
|
|
call = lambda fun, x: fun(x)
|
|
|
|
ex = lambda x: x**2+y*z
|
|
|
|
sigma1 = lambda x: x*y+1
|
|
|
|
Gc = self.M.gridCC
|
|
sigma = call(sigma1, Gc)
|
|
analytic = 128011./5 # Found using sympy. y=12, z=5
|
|
|
|
if self.location == 'faces':
|
|
F = call(ex, self.M.gridFx)
|
|
if self.invProp:
|
|
A = self.M.getFaceInnerProduct(1/sigma, invProp=True)
|
|
else:
|
|
A = self.M.getFaceInnerProduct(sigma)
|
|
numeric = F.T.dot(A.dot(F))
|
|
|
|
err = np.abs(numeric - analytic)
|
|
return err
|
|
|
|
def test_order1_faces(self):
|
|
self.name = "1D Face Inner Product"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 1
|
|
self.invProp = False
|
|
self.orderTest()
|
|
|
|
def test_order1_faces_invProp(self):
|
|
self.name = "1D Face Inner Product - invProp"
|
|
self.location = 'faces'
|
|
self.sigmaTest = 1
|
|
self.invProp = True
|
|
self.orderTest()
|
|
|
|
|
|
if __name__ == '__main__':
|
|
unittest.main()
|
|
|
|
###################################################
|
|
#### Uncomment to Reevaluate the InnerProducts ####
|
|
###################################################
|
|
|
|
# if __name__ == '__main__':
|
|
# import sympy
|
|
|
|
# x,y,z = sympy.symbols(['x','y','z'])
|
|
# ex = x**2+y*z
|
|
# ey = (z**2)*x+y*z
|
|
# ez = y**2+x*z
|
|
# e = sympy.Matrix([ex,ey,ez])
|
|
|
|
# sigma1 = x*y+1
|
|
# sigma2 = x*z+2
|
|
# sigma3 = 3+z*y
|
|
# sigma4 = 0.1*x*y*z
|
|
# sigma5 = 0.2*x*y
|
|
# sigma6 = 0.1*z
|
|
|
|
# S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
|
|
# S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
|
|
# S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
|
|
|
|
# print '3D'
|
|
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
|
|
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
|
|
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
|
|
|
|
|
|
# z = 5
|
|
# ex = x**2+y*z
|
|
# ey = (z**2)*x+y*z
|
|
# e = sympy.Matrix([ex,ey])
|
|
|
|
# sigma1 = x*y+1
|
|
# sigma2 = x*z+2
|
|
# sigma3 = 3+z*y
|
|
|
|
# S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
|
|
# S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
|
|
# S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
|
|
|
|
# print '2D'
|
|
# print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
|
|
# print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
|
|
# print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
|
|
|
|
# y = 12
|
|
# z = 5
|
|
# ex = x**2+y*z
|
|
# e = ex
|
|
|
|
# sigma1 = x*y+1
|
|
|
|
# print '1D'
|
|
# print sympy.integrate(e*sigma1*e, (x,0,1))
|