mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-06 05:16:51 +08:00
Merge branch 'eldadswork' of https://bitbucket.org/rcockett/simpeg into LOM
Conflicts: SimPEG/utils.py
This commit is contained in:
Rowan Cockett
@@ -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)
|
||||
@@ -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()
|
||||
+35
-2
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user