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

Test for edge inner products working.

Browse Source
This commit is contained in:
Rowan Cockett
2013-07-26 12:11:45 -07:00
parent 87331c4c92
commit 150cbc7df3
2 changed files with 86 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/getEdgeInnerProducts.py
+61 -61
View File
@@ -6,7 +6,7 @@ from numpy import *
from TensorMesh import *
# [A] = getEdgeInnerProduct(X,Y,Z,sigma)
#
#
# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
# / /
@@ -35,147 +35,147 @@ from TensorMesh import *
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
def subarray(T,i1,i2,i3):
return take(take(take(T,i1,0),i2,1),i3,2)
def subarray(T, i1, i2, i3):
return take(take(take(T, i1, 0), i2, 1), i3, 2)
def getEdgeInnerProduct(mesh,sigma):
def getEdgeInnerProduct(mesh, sigma):
h = mesh.h
m = array([size(h[0]),size(h[1]),size(h[2])])
m = array([size(h[0]), size(h[1]), size(h[2])])
nc = prod(m)
me1 = m + array([0, 1, 1]); ne1 = prod(me1)
me2 = m + array([1, 0, 1]); ne2 = prod(me2)
me3 = m + array([1, 1, 0]); ne3 = prod(me3)
i = int64(linspace(0,m[0]-1,m[0]))
j = int64(linspace(0,m[1]-1,m[1]))
k = int64(linspace(0,m[2]-1,m[2]))
ii,jj,kk = ndgrid(i,j,k,vector=False)
ii,jj,kk = ndgrid(i,j,k,vector=False)
ii = mkvc(ii); jj = mkvc(jj); kk = mkvc(kk)
## --------
# no | node | e1 | e2 | e3
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
ind1 = sub2ind(me1,c_[ii,jj,kk])
ind1 = sub2ind(me1,c_[ii,jj,kk])
ind2 = sub2ind(me2,c_[ii,jj,kk]) + ne1
ind3 = sub2ind(me3,c_[ii,jj,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P000 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
## --------
# no | node | e1 | e2 | e3
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
ind1 = sub2ind(me1,c_[ii,jj,kk])
ind1 = sub2ind(me1,c_[ii,jj,kk])
ind2 = sub2ind(me2,c_[ii+1,jj,kk]) + ne1
ind3 = sub2ind(me3,c_[ii+1,jj,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P100 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
## --------
# no | node | e1 | e2 | e3
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
ind1 = sub2ind(me1,c_[ii,jj+1,kk])
ind1 = sub2ind(me1,c_[ii,jj+1,kk])
ind2 = sub2ind(me2,c_[ii,jj,kk]) + ne1
ind3 = sub2ind(me3,c_[ii,jj+1,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P010 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
## --------
# no | node | e1 | e2 | e3
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
ind1 = sub2ind(me1,c_[ii,jj+1,kk])
ind1 = sub2ind(me1,c_[ii,jj+1,kk])
ind2 = sub2ind(me2,c_[ii+1,jj,kk]) + ne1
ind3 = sub2ind(me3,c_[ii+1,jj+1,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P110 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
######
## --------
# no | node | e1 | e2 | e3
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
ind1 = sub2ind(me1,c_[ii,jj,kk+1])
ind1 = sub2ind(me1,c_[ii,jj,kk+1])
ind2 = sub2ind(me2,c_[ii,jj,kk+1]) + ne1
ind3 = sub2ind(me3,c_[ii,jj,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P001 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
## --------
# no | node | e1 | e2 | e3
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k+1
ind1 = sub2ind(me1,c_[ii,jj,kk+1])
ind1 = sub2ind(me1,c_[ii,jj,kk+1])
ind2 = sub2ind(me2,c_[ii+1,jj,kk+1]) + ne1
ind3 = sub2ind(me3,c_[ii+1,jj,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P101 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
## --------
# no | node | e1 | e2 | e3
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k+1
ind1 = sub2ind(me1,c_[ii,jj+1,kk+1])
ind1 = sub2ind(me1,c_[ii,jj+1,kk+1])
ind2 = sub2ind(me2,c_[ii,jj,kk+1]) + ne1
ind3 = sub2ind(me3,c_[ii,jj+1,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P011 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
## --------
# no | node | e1 | e2 | e3
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k+1
ind1 = sub2ind(me1,c_[ii,jj+1,kk+1])
ind1 = sub2ind(me1,c_[ii,jj+1,kk+1])
ind2 = sub2ind(me2,c_[ii+1,jj,kk+1]) + ne1
ind3 = sub2ind(me3,c_[ii+1,jj+1,kk]) + ne1 + ne2
IND = vstack((vstack((ind1,ind2)),ind3))
IND = array(IND).flatten()
P111 = sparse.coo_matrix((ones(3*nc),(linspace(0,3*nc-1,3*nc),IND)),shape=(3*nc,ne1+ne2+ne3)).tocsr()
# Cell volume
v = sqrt(mesh.vol)
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])))
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))
v3 = r_[v,v,v]
v = sqrt(mesh.vol)
v3 = r_[v, v, v]
V = sdiag(v3)*Sigma*sdiag(v3)
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 = 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__':
h = [array([1,2,3,4]),array([1,2,1,4,2]),array([1,1,4,1])]
h = [array([1, 2, 3, 4]), array([1, 2, 1, 4, 2]), array([1, 1, 4, 1])]
mesh = TensorMesh(h)
sigma = ones((mesh.nC,6))
A = getEdgeInnerProduct(mesh,sigma)
sigma = ones((mesh.nC, 6))
A = getEdgeInnerProduct(mesh, sigma)
SimPEG/tests/test_massMatrices.py
+25 -22
View File
@@ -2,18 +2,21 @@ import numpy as np
import unittest
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
from OrderTest import OrderTest
from scipy.sparse.linalg import dsolve
from getEdgeInnerProducts import getEdgeInnerProducts
from getEdgeInnerProducts import *
class TestNodalGrad(OrderTest):
name = "Nodal Gradient"
class TestEdgeInnerProduct(OrderTest):
"""Integrate a function over a unit cube domain."""
name = "Edge Inner Product"
meshSizes = [4, 8, 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
@@ -24,25 +27,25 @@ class TestNodalGrad(OrderTest):
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 = ex(self.M.gridEx[:, 0],self.M.gridEx[:, 1],self.M.gridEx[:, 2])
Ey = ey(self.M.gridEy[:, 0],self.M.gridEy[:, 1],self.M.gridEy[:, 2])
Ez = ez(self.M.gridEz[:, 0],self.M.gridEz[:, 1],self.M.gridEz[:, 2])
E = np.r_[Ex,Ey,Ez]
Ex = call(ex, self.M.gridEx)
Ey = call(ey, self.M.gridEy)
Ez = call(ez, self.M.gridEz)
E = np.matrix(mkvc(np.r_[Ex, Ey, Ez], 2))
Gc = self.M.gridCC
sigma = np.c_[sigma1(Gc[:,0],Gc[:,1],Gc[:,2]),
sigma2(Gc[:,0],Gc[:,1],Gc[:,2]),
sigma3(Gc[:,0],Gc[:,1],Gc[:,2]),
sigma4(Gc[:,0],Gc[:,1],Gc[:,2]),
sigma5(Gc[:,0],Gc[:,1],Gc[:,2]),
sigma6(Gc[:,0],Gc[:,1],Gc[:,2])]
A = getEdgeInnerProducts(self.M, sigma)
err = np.abs(E.T*A*E - 69881./21600)
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()