From d6213db486fb7ed81d1d3c6f6d849260d1f6b91e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 7 Feb 2014 14:03:37 -0800 Subject: [PATCH 1/6] Initial commit for slicing in matplotlib --- SimPEG/Mesh/BaseMesh.py | 2 +- SimPEG/Mesh/TensorView.py | 86 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+), 1 deletion(-) diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index 2dd82808..57d6f482 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -44,7 +44,7 @@ class BaseMesh(object): For example, you have a face variable, and you want the x component of it reshaped to a 3D matrix. - Mesh.r can fulfil your dreams:: + Mesh.r can fulfill your dreams:: mesh.r(V, 'F', 'Fx', 'M') | | | { How: 'M' or ['V'] for a matrix (ndgrid style) or a vector (n x dim) } diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index 7446aecd..8e44ee10 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -238,6 +238,92 @@ class TensorView(object): if showIt: plt.show() return ph + def plotSlice(self, v, vType='CC', + normal='Z', ind=None, grid=False, view='real', + ax=None, clim=None, showIt=False, + gridOpts={'color':'k'} + ): + viewOpts = ['real','imag','abs','vec'] + normalOpts = ['X', 'Y', 'Z'] + vTypeOpts = ['CC','F','E'] + + # Some user error checking + assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) + assert self.dim == 3, 'Must be a 3D mesh.' + assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts) + assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts) + assert type(grid) is bool, 'grid must be a boolean' + + szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension + if ind is None: ind = int(szSliceDim/2) + assert type(ind) in [int, long], 'ind must be an integer' + + if ax is None: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + else: + assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes" + fig = ax.figure + + # The slicing and plotting code!! + + def getIndSlice(v): + if normal == 'X': v = v[ind,:,:] + elif normal == 'Y': v = v[:,ind,:] + elif normal == 'Z': v = v[:,:,ind] + return v + + def doSlice(v): + if vType == 'CC': + return getIndSlice(self.r(v,'CC','CC','M')) + # Now just deal with 'F' and 'E' + aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC') + v = getattr(self,aveOp)*v # average to cell centers (might be a vector) + if view == 'vec': + v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M') + outSlice = [] + if 'X' not in normal: outSlice.append(getIndSlice(v[0])) + if 'Y' not in normal: outSlice.append(getIndSlice(v[1])) + if 'Z' not in normal: outSlice.append(getIndSlice(v[2])) + return outSlice + else: + return getIndSlice(self.r(v,'CC','CC','M')) + + h2d = [] + if 'X' not in normal: h2d.append(self.hx) + if 'Y' not in normal: h2d.append(self.hy) + if 'Z' not in normal: h2d.append(self.hz) + tM = Mesh.TensorMesh(h2d) #: Temp Mesh + + out = () + if view in ['real','imag','abs']: + v = getattr(np,view)(v) # e.g. np.real(v) + v = doSlice(v) + if clim is None: + clim = [v.min(),v.max()] + out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1]),) + elif view in ['vec']: + U, V = doSlice(v) + if clim is None: + clim = [v.min(),v.max()] + out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),) + out += (quiver( tM.gridCC[:,0], tM.gridCC[:,1], U, V),) + + if grid: + xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten() + xYGrid = np.c_[tM.vectorNy[0]*np.ones(tM.nNx),tM.vectorNy[-1]*np.ones(tM.nNx),np.nan*np.ones(tM.nNx)].flatten() + yXGrid = np.c_[tM.vectorNx[0]*np.ones(tM.nNy),tM.vectorNx[-1]*np.ones(tM.nNy),np.nan*np.ones(tM.nNy)].flatten() + yYGrid = np.c_[tM.vectorNy,tM.vectorNy,np.nan*np.ones(tM.nNy)].flatten() + out += (ax.plot(np.r_[xXGrid,yXGrid],np.r_[xYGrid,yYGrid],**gridOpts)[0],) + + ax.set_xlabel('y' if normal == 'X' else 'x') + ax.set_ylabel('y' if normal == 'Z' else 'z') + ax.set_title('Slice %d' % ind) + + if showIt: plt.show() + return out + def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. From 1c7c79344ba811f51118fdc8e13942060ebff237 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 7 Feb 2014 14:54:15 -0800 Subject: [PATCH 2/6] Option for normalized vector plot --- SimPEG/Mesh/TensorView.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index 8e44ee10..b0ae3915 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -243,7 +243,7 @@ class TensorView(object): ax=None, clim=None, showIt=False, gridOpts={'color':'k'} ): - viewOpts = ['real','imag','abs','vec'] + viewOpts = ['real','imag','abs','vec','|vec|'] normalOpts = ['X', 'Y', 'Z'] vTypeOpts = ['CC','F','E'] @@ -258,6 +258,11 @@ class TensorView(object): if ind is None: ind = int(szSliceDim/2) assert type(ind) in [int, long], 'ind must be an integer' + normalizeVector = False + if view == '|vec|': + view = 'vec' + normalizeVector = True + if ax is None: fig = plt.figure(1) fig.clf() @@ -308,6 +313,9 @@ class TensorView(object): if clim is None: clim = [v.min(),v.max()] out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),) + if normalizeVector: + length = np.sqrt(U**2+V**2) + U, V = U/length, V/length out += (quiver( tM.gridCC[:,0], tM.gridCC[:,1], U, V),) if grid: From 530de6cb5c7bb6c61672dbe66efbe0fafcb4e219 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 18:31:30 -0800 Subject: [PATCH 3/6] FaceInnerProduct for 1D --- SimPEG/Mesh/BaseMesh.py | 4 +- SimPEG/Mesh/InnerProducts.py | 65 ++++++++++++++++++++++---- SimPEG/Mesh/LogicallyOrthogonalMesh.py | 1 + SimPEG/Mesh/TreeMesh.py | 2 + 4 files changed, 60 insertions(+), 12 deletions(-) diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index aa487415..2d70093c 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -243,7 +243,7 @@ class BaseMesh(object): :return: projected face vector """ assert type(fV) == np.ndarray, 'fV must be an ndarray' - assert len(fV.shape) == 2 and fV.shape[0] == np.sum(self.nF) and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)' + assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)' return np.sum(fV*self.normals, 1) def projectEdgeVector(self, eV): @@ -255,7 +255,7 @@ class BaseMesh(object): :return: projected edge vector """ assert type(eV) == np.ndarray, 'eV must be an ndarray' - assert len(eV.shape) == 2 and eV.shape[0] == np.sum(self.nE) and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)' + assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)' return np.sum(eV*self.tangents, 1) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 427e5ee2..86afecaf 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -142,7 +142,14 @@ class InnerProducts(object): Note that this is completed for each cell in the mesh at the same time. """ - if M.dim == 2: + if M.dim == 1: + v = np.sqrt(0.5*M.vol) + V1 = sdiag(v) # We will multiply on each side to keep symmetry + + Px = _getFacePx(M) + P000 = V2*Px('fXm') + P100 = V2*Px('fXp') + elif M.dim == 2: # Square root of cell volume multiplied by 1/4 v = np.sqrt(0.25*M.vol) V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry @@ -168,9 +175,13 @@ class InnerProducts(object): P111 = V3*Pxxx('fXp', 'fYp', 'fZp') Mu = _makeTensor(M, mu) - A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110 - P = [P000, P100, P010, P110] - if M.dim == 3: + A = P000.T*Mu*P000 + P100.T*Mu*P100 + P = [P000, P100] + + if M.dim > 1: + A = A + P010.T*Mu*P010 + P110.T*Mu*P110 + P += [P010, P110] + if M.dim > 2: A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111 P += [P001, P101, P011, P111] if returnP: @@ -246,8 +257,10 @@ class InnerProducts(object): Note that this is completed for each cell in the mesh at the same time. """ + if M.dim == 1: + raise NotImplementedError('getEdgeInnerProduct not implemented for 1D') # We will multiply by V on each side to keep symmetry - if M.dim == 2: + elif M.dim == 2: # Square root of cell volume multiplied by 1/4 v = np.sqrt(0.25*M.vol) V = sdiag(np.r_[v, v]) @@ -304,7 +317,13 @@ def _makeTensor(M, sigma): if sigma is None: # default is ones sigma = np.ones((M.nC, 1)) - if M.dim == 2: + if M.dim == 1: + if sigma.size == M.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(sigma) + else: + raise Exception('Unexpected shape of sigma') + elif M.dim == 2: if sigma.size == M.nC: # Isotropic! sigma = mkvc(sigma) # ensure it is a vector. Sigma = sdiag(np.r_[sigma, sigma]) @@ -314,6 +333,8 @@ def _makeTensor(M, sigma): row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) Sigma = sp.vstack((row1, row2)) + else: + raise Exception('Unexpected shape of sigma') elif M.dim == 3: if sigma.size == M.nC: # Isotropic! sigma = mkvc(sigma) # ensure it is a vector. @@ -325,8 +346,15 @@ def _makeTensor(M, sigma): row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) Sigma = sp.vstack((row1, row2, row3)) + else: + raise Exception('Unexpected shape of sigma') return Sigma + +def _getFacePx(M): + assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh' + return _getFacePx_Rectangular(M) + def _getFacePxx(M): if M._meshType == 'TREE': return M._getFacePxx @@ -351,6 +379,23 @@ def _getEdgePxxx(M): return _getEdgePxxx_Rectangular(M) +def _getFacePx_Rectangular(M): + """Returns a function for creating projection matrices + + """ + i = np.int64(range(M.nCx)) + + def Px(xFace): + """ + xFace is 'fXp' or 'fXm' + """ + posFx = 0 if xFace == 'fXm' else 1 + IND = ii + posFx + PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF)) + return PX + + return Px + def _getFacePxx_Rectangular(M): """returns a function for creating projection matrices @@ -408,7 +453,7 @@ def _getFacePxx_Rectangular(M): IND = np.r_[ind1, ind2].flatten() - PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nF))) + PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF)) if M._meshType == 'LOM': I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]), @@ -465,7 +510,7 @@ def _getFacePxxx_Rectangular(M): IND = np.r_[ind1, ind2, ind3].flatten() - PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nF))).tocsr() + PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr() if M._meshType == 'LOM': I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]), @@ -500,7 +545,7 @@ def _getEdgePxx_Rectangular(M): IND = np.r_[ind1, ind2].flatten() - PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nE))).tocsr() + PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr() if M._meshType == 'LOM': I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]), @@ -543,7 +588,7 @@ def _getEdgePxxx_Rectangular(M): IND = np.r_[ind1, ind2, ind3].flatten() - PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nE))).tocsr() + PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr() if M._meshType == 'LOM': I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]), diff --git a/SimPEG/Mesh/LogicallyOrthogonalMesh.py b/SimPEG/Mesh/LogicallyOrthogonalMesh.py index 480c4f06..7b4ccd73 100644 --- a/SimPEG/Mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/Mesh/LogicallyOrthogonalMesh.py @@ -32,6 +32,7 @@ class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts, def __init__(self, nodes): assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray" + assert len(nodes) > 1, "len(node) must be greater than 1" for i, nodes_i in enumerate(nodes): assert type(nodes_i) == np.ndarray, ("nodes[%i] is not a numpy array." % i) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 3d2c2794..ae83eae8 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -679,6 +679,8 @@ class TreeMesh(InnerProducts, BaseMesh): def __init__(self, h_in, x0=None): assert type(h_in) is list, 'h_in must be a list' + assert len(h_in) > 1, "len(h_in) must be greater than 1" + h = range(len(h_in)) for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: From abd4f47dfcb863579e1b9034bfdddbb3e1e3e4e1 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 18:39:14 -0800 Subject: [PATCH 4/6] bug fix. --- SimPEG/Mesh/InnerProducts.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 86afecaf..17747f10 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -147,8 +147,8 @@ class InnerProducts(object): V1 = sdiag(v) # We will multiply on each side to keep symmetry Px = _getFacePx(M) - P000 = V2*Px('fXm') - P100 = V2*Px('fXp') + P000 = V1*Px('fXm') + P100 = V1*Px('fXp') elif M.dim == 2: # Square root of cell volume multiplied by 1/4 v = np.sqrt(0.25*M.vol) From 5b09891601ad295ed045f22066290d240d4d53eb Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 18:40:32 -0800 Subject: [PATCH 5/6] bug fix. --- SimPEG/Mesh/InnerProducts.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 17747f10..7a1a7869 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -383,7 +383,7 @@ def _getFacePx_Rectangular(M): """Returns a function for creating projection matrices """ - i = np.int64(range(M.nCx)) + ii = np.int64(range(M.nCx)) def Px(xFace): """ From 9a53545d0a5e20efc104e82ce51c642dce2149ba Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 19:10:58 -0800 Subject: [PATCH 6/6] 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))