From 9a53545d0a5e20efc104e82ce51c642dce2149ba Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 19:10:58 -0800 Subject: [PATCH] testFaceInnerProduct1D --- SimPEG/Tests/test_massMatrices.py | 130 ++++++++++++++++++++++-------- 1 file changed, 98 insertions(+), 32 deletions(-) diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 79626ca1..0db8c007 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -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))