mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-04 02:28:46 +08:00
mass matrices for anisotropic sigma
This commit is contained in:
ehaber99
@@ -0,0 +1,181 @@
|
||||
from scipy.sparse import linalg
|
||||
from scipy import sparse
|
||||
from sputils import *
|
||||
from utils import *
|
||||
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)
|
||||
# / /
|
||||
# / / |
|
||||
# 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
|
||||
|
||||
|
||||
def subarray(T,i1,i2,i3):
|
||||
return take(take(take(T,i1,0),i2,1),i3,2)
|
||||
|
||||
|
||||
def getEdgeInnerProduct(mesh,sigma):
|
||||
|
||||
h = mesh.h
|
||||
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 = 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])
|
||||
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])
|
||||
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])
|
||||
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])
|
||||
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])
|
||||
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])
|
||||
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])
|
||||
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])
|
||||
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])))
|
||||
Sigma = sp.vstack((row1, row2, row3))
|
||||
|
||||
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 = 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])]
|
||||
mesh = TensorMesh(h)
|
||||
sigma = ones((mesh.nC,6))
|
||||
A = getEdgeInnerProduct(mesh,sigma)
|
||||
@@ -0,0 +1,94 @@
|
||||
import numpy as np
|
||||
from scipy import sparse as sp
|
||||
from sputils import sdiag, speye, kron3, spzeros
|
||||
from utils import mkvc
|
||||
|
||||
|
||||
|
||||
def getEdgeMassMatrix(sigma,mesh):
|
||||
"""Get anisotropic mass matrix"""
|
||||
|
||||
n = array([size(mesh.h[0]),size(mesh.h[1]),size(mesh.h[2])])
|
||||
nx = prod(n + [1, 0, 0])
|
||||
ex = reshape(arange(0,nx),[n[0]+1,n[1],n[2]])
|
||||
ny = prod(n + [0, 1, 0])
|
||||
ey = reshape(arange(0,ny),[n[0],n[1]+1,n[2]])
|
||||
nz = prod(n + [0, 0, 1]);
|
||||
ez = reshape(arange(0,nz),[n[0],n[1],n[2]+1])
|
||||
|
||||
|
||||
i = arange(0,n[0]-1); j = arange(0,n[1]-1); k = arange(0,n[2]-1)
|
||||
|
||||
# corner i,j,k
|
||||
Px1 = take(ex,[i,j,k]); Py1 = take(ey,[i,j,k]); Pz1 = take(ez,[i,j,k])
|
||||
# corner i+1,j,k
|
||||
Px2 = take(ex,[i,j,k]); Py2 = take(ey,[i+1,j,k]); Pz2 = take(ez,[i+1,j,k])
|
||||
# corner i,j+1,k
|
||||
Px3 = take(ex,[i,j+1,k]); Py3 = take(ey,[i,j,k]); Pz3 = take(ez,[i,j+1,k])
|
||||
# corner i+1,j+1,k
|
||||
Px4 = take(ex,[i,j+1,k]); Py4 = take(ey,[i+1,j,k]); Pz4 = take(ez,[i+1,j+1,k]);
|
||||
|
||||
# corner i,j,k+1
|
||||
Px5 = take(ex,[i,j,k+1]); Py5 = take(ey,[i,j,k+1]); Pz5 = take(ez,[i,j,k])
|
||||
# corner i+1,j,k+1
|
||||
Px6 = take(ex,[i,j,k+1]); Py6 = take(ey,[i+1,j,k+1]); Pz6 = take(ez,[i+1,j,k])
|
||||
# corner i,j+1,k+1
|
||||
Px7 = take(ex,[i,j+1,k+1]); Py7 = take(ey,[i,j,k+1]); Pz7 = take(ez,[i,j+1,k])
|
||||
# corner i+1,j+1,k+1
|
||||
Px8 = take(ex,[i,j+1,k+1]); Py8 = take(ey,[i+1,j,k+1]); Pz8 = take(ez,[i+1,j+1,k])
|
||||
|
||||
|
||||
nx1 = size(Px1); ny1 = size(Py1); nz1 = size(Pz1)
|
||||
#sparse.coo_matrix((V,(I,J)),shape=(4,4))
|
||||
P1 = block_diag(( sparse.coo_matrix(arange(0,nx1),Px1(:), e(nx1), nx1,nx),
|
||||
sparse.coo_matrix(arange(0,ny1),Py1(:),e(ny1), ny1,ny),
|
||||
sparse.coo_matrix(arange(0,nz1),Pz1(:),e(nz1), nz1,nz)))
|
||||
|
||||
nx2 = numel(Px2); ny2 = numel(Py2); nz2 = numel(Pz2);
|
||||
P2 = blkdiag( sparse(1:nx2,Px2(:), e(nx2), nx2,nx) , ...
|
||||
sparse(1:ny2,Py2(:),e(ny2), ny2,ny), ...
|
||||
sparse(1:nz2,Pz2(:),e(nz2), nz2,nz));
|
||||
|
||||
nx3 = numel(Px3); ny3 = numel(Py3); nz3 = numel(Pz3);
|
||||
P3 = blkdiag( sparse(1:nx3,Px3(:), e(nx3), nx3,nx) , ...
|
||||
sparse(1:ny3,Py3(:),e(ny3), ny3,ny), ...
|
||||
sparse(1:nz3,Pz3(:),e(nz3), nz3,nz));
|
||||
|
||||
nx4 = numel(Px4); ny4 = numel(Py4); nz4 = numel(Pz4);
|
||||
P4 = blkdiag( sparse(1:nx4,Px4(:), e(nx4), nx4,nx) , ...
|
||||
sparse(1:ny4,Py4(:), e(ny4), ny4,ny), ...
|
||||
sparse(1:nz4,Pz4(:), e(nz4), nz4,nz));
|
||||
|
||||
nx5 = numel(Px5); ny5 = numel(Py5); nz5 = numel(Pz5);
|
||||
P5 = blkdiag( sparse(1:nx5,Px5(:), e(nx5), nx5,nx) , ...
|
||||
sparse(1:ny5,Py5(:), e(ny5), ny5,ny), ...
|
||||
sparse(1:nz5,Pz5(:), e(nz5), nz5,nz));
|
||||
|
||||
nx6 = numel(Px6); ny6 = numel(Py6); nz6 = numel(Pz6);
|
||||
P6 = blkdiag( sparse(1:nx6,Px6(:), e(nx6), nx6,nx) , ...
|
||||
sparse(1:ny6,Py6(:), e(ny6), ny6,ny), ...
|
||||
sparse(1:nz6,Pz6(:), e(nz6), nz6,nz));
|
||||
|
||||
nx7 = numel(Px7); ny7 = numel(Py7); nz7 = numel(Pz7);
|
||||
P7 = blkdiag( sparse(1:nx7,Px7(:), e(nx7), nx7,nx) , ...
|
||||
sparse(1:ny7,Py7(:), e(ny7), ny7,ny), ...
|
||||
sparse(1:nz7,Pz7(:), e(nz7), nz7,nz));
|
||||
|
||||
nx8 = numel(Px8); ny8 = numel(Py8); nz8 = numel(Pz8);
|
||||
P8 = blkdiag( sparse(1:nx8,Px8(:), e(nx8), nx8,nx) , ...
|
||||
sparse(1:ny8,Py8(:), e(ny8), ny8,ny), ...
|
||||
sparse(1:nz8,Pz8(:), e(nz8), nz8,nz));
|
||||
|
||||
V = sdiag(sqrt([v(:); v(:); v(:)]));
|
||||
|
||||
# generate the conductivity
|
||||
S = [sdiag(sig(:,1)) , sdiag(sig(:,4)) , sdiag(sig(:,5)); ...
|
||||
sdiag(sig(:,4)) , sdiag(sig(:,2)) , sdiag(sig(:,6)); ...
|
||||
sdiag(sig(:,5)) , sdiag(sig(:,6)) , sdiag(sig(:,3))];
|
||||
|
||||
# scale by the volume
|
||||
S = V*S*V;
|
||||
|
||||
M = 1/8*(P1'*S*P1 + P2'*S*P2 + P3'*S*P3 + P4'*S*P4 + ...
|
||||
P5'*S*P5 + P6'*S*P6 + P7'*S*P7 + P8'*S*P8);
|
||||
|
||||
@@ -0,0 +1,8 @@
|
||||
import numpy as np
|
||||
|
||||
def getSubArray(A,ind):
|
||||
"""subArray"""
|
||||
i = ind[0]; j = ind[1]; k = ind[2]
|
||||
|
||||
return A[i,:,:][:,j,:][:,:,k]
|
||||
|
||||
@@ -0,0 +1,48 @@
|
||||
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
|
||||
|
||||
|
||||
class TestNodalGrad(OrderTest):
|
||||
name = "Nodal Gradient"
|
||||
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
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 = 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]
|
||||
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)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
+30
-1
@@ -1,5 +1,5 @@
|
||||
import numpy as np
|
||||
|
||||
from numpy import *
|
||||
|
||||
def reshapeF(x, size):
|
||||
return np.reshape(x, size, order='F')
|
||||
@@ -97,3 +97,32 @@ def ndgrid(*args, **kwargs):
|
||||
return np.c_[X1, X2, X3]
|
||||
else:
|
||||
return XYZ[2], XYZ[1], XYZ[0]
|
||||
|
||||
|
||||
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