wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-03 19:14:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

testFaceInnerProduct1D

Browse Source
This commit is contained in:
rowanc1
2014-02-16 19:10:58 -08:00
parent 5b09891601
commit 9a53545d0a
1 changed files with 98 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Tests/test_massMatrices.py
+98 -32
View File
@@ -3,32 +3,6 @@ import unittest
from TestUtils import OrderTest
# MATLAB code:
# syms x y z
# ex = x.^2+y.*z;
# ey = (z.^2).*x+y.*z;
# ez = y.^2+x.*z;
# e = [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 = [sigma1,0,0;0,sigma1,0;0,0,sigma1];
# S2 = [sigma1,0,0;0,sigma2,0;0,0,sigma3];
# S3 = [sigma1,sigma4,sigma5;sigma4,sigma2,sigma6;sigma5,sigma6,sigma3];
# i1 = int(int(int(e.'*S1*e,x,0,1),y,0,1),z,0,1);
# i2 = int(int(int(e.'*S2*e,x,0,1),y,0,1),z,0,1);
# i3 = int(int(int(e.'*S3*e,x,0,1),y,0,1),z,0,1);
class TestInnerProducts(OrderTest):
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
@@ -54,14 +28,14 @@ class TestInnerProducts(OrderTest):
Gc = self.M.gridCC
if self.sigmaTest == 1:
sigma = np.c_[call(sigma1, Gc)]
analytic = 647./360 # Found using matlab symbolic toolbox.
analytic = 647./360 # Found using sympy.
elif self.sigmaTest == 3:
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
analytic = 37./12 # Found using matlab symbolic toolbox.
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 matlab symbolic toolbox.
analytic = 69881./21600 # Found using sympy.
if self.location == 'edges':
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
@@ -143,13 +117,13 @@ class TestInnerProducts2D(OrderTest):
Gc = self.M.gridCC
if self.sigmaTest == 1:
sigma = np.c_[call(sigma1, Gc)]
analytic = 144877./360 # Found using matlab symbolic toolbox. z=5
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 matlab symbolic toolbox. z=5
analytic = 189959./120 # Found using sympy. z=5
elif self.sigmaTest == 3:
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
analytic = 781427./360 # Found using matlab symbolic toolbox. z=5
analytic = 781427./360 # Found using sympy. z=5
if self.location == 'edges':
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
@@ -206,5 +180,97 @@ class TestInnerProducts2D(OrderTest):
self.orderTest()
class TestInnerProducts1D(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)
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.orderTest()
if __name__ == '__main__':
unittest.main()
if __name__ == '__main__' and False:
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))