diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py new file mode 100644 index 00000000..68423a3f --- /dev/null +++ b/SimPEG/getEdgeInnerProducts.py @@ -0,0 +1,87 @@ +from scipy import sparse as sp +from sputils import sdiag +from utils import sub2ind, ndgrid, mkvc +import numpy as np + + +def getEdgeInnerProduct(mesh, sigma): + + m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) + nc = mesh.nC + + i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) + + iijjkk = ndgrid(i, j, k) + ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] + + def Pxxx(pos): + ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) + ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0] + ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1] + + IND = np.r_[ind1, ind2, ind3].flatten() + + return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr() + + # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) + # / / + # / / | + # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) + # / / | + # / / | + # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) + # | | | + # | | node(i+1,j+1,k+1) + # | | / + # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) + # | | / + # | | / + # | |/ + # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) + + # no | node | e1 | e2 | e3 + # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k + # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k + # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k + # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k + # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k + # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k + # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k + # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k + P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) + P100 = Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]]) + P010 = Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]]) + P110 = Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]]) + P001 = Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]]) + P101 = Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]]) + P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]]) + P111 = Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) + + if sigma.size == mesh.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(np.r_[sigma, sigma, sigma]) + elif sigma.shape[1] == 3: # Diagonal tensor + Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) + elif sigma.shape[1] == 6: # Fully anisotropic + row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) + 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)) + + # Cell volume + v = np.sqrt(mesh.vol) + v3 = np.r_[v, v, v] + V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry + + A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111 + + A = 0.125*A + + return A + +if __name__ == '__main__': + from TensorMesh import TensorMesh + h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])] + mesh = TensorMesh(h) + sigma = np.ones((mesh.nC, 6)) + A = getEdgeInnerProduct(mesh, sigma) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py new file mode 100644 index 00000000..0aae0615 --- /dev/null +++ b/SimPEG/tests/test_massMatrices.py @@ -0,0 +1,49 @@ +import numpy as np +import unittest +from OrderTest import OrderTest +import sys +sys.path.append('../') +from getEdgeInnerProducts import * + + +class TestEdgeInnerProduct(OrderTest): + """Integrate an edge function over a unit cube domain using edgeInnerProducts.""" + + name = "Edge Inner Product" + + 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 + + Ex = call(ex, self.M.gridEx) + Ey = call(ey, self.M.gridEy) + Ez = call(ez, self.M.gridEz) + + E = np.matrix(np.r_[Ex, Ey, Ez]).T + Gc = self.M.gridCC + sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), + call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)] + + A = getEdgeInnerProduct(self.M, sigma) + numeric = E.T*A*E + analytic = 69881./21600 # Found using matlab symbolic toolbox. + err = np.abs(numeric - analytic) + return err + + def test_order(self): + self.orderTest() + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/utils.py b/SimPEG/utils.py index 84bf11bf..55acf9b0 100644 --- a/SimPEG/utils.py +++ b/SimPEG/utils.py @@ -1,4 +1,5 @@ import numpy as np +from numpy import * def mkvc(x, numDims=1): @@ -119,7 +120,39 @@ def volTetra(xyz, A, B, C, D): BD = xyz[B, :] - xyz[D, :] CD = xyz[C, :] - xyz[D, :] - - V = (BD[:, 0]*CD[:, 1] - BD[:, 1]*CD[:, 0])*AD[:, 2] - (BD[:, 0]*CD[:, 2] - BD[:, 2]*CD[:, 0])*AD[:, 1] + (BD[:, 1]*CD[:, 2] - BD[:, 2]*CD[:, 1])*AD[:, 0] return V/6 + + +def getSubArray(A, ind): + """subArray""" + return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]] + + +def ind2sub(shape, ind): + """From the given shape, returns the subscrips of the given index""" + revshp = [] + revshp.extend(shape) + mult = [1] + for i in range(0, len(revshp)-1): + mult.extend([mult[i]*revshp[i]]) + mult = array(mult).reshape(len(mult)) + + sub = [] + + for i in range(0, len(shape)): + sub.extend([math.floor(ind / mult[i])]) + ind = ind - (math.floor(ind/mult[i]) * mult[i]) + return sub + + +def sub2ind(shape, subs): + """From the given shape, returns the index of the given subscript""" + revshp = list(shape) + mult = [1] + for i in range(0, len(revshp)-1): + mult.extend([mult[i]*revshp[i]]) + mult = array(mult).reshape(len(mult), 1) + + idx = dot((subs), (mult)) + return idx