mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 22:23:39 +08:00
boundaryCondition initial work.
This commit is contained in:
rowanc1
@@ -460,6 +460,104 @@ class DiffOperators(object):
|
||||
_edgeCurl = None
|
||||
edgeCurl = property(**edgeCurl())
|
||||
|
||||
def getBCProjWF(self, BC, discretization='CC'):
|
||||
"""
|
||||
|
||||
The weak form boundary condition projection matrices.
|
||||
|
||||
Examples::
|
||||
|
||||
BC = 'neumann' # Neumann in all directions
|
||||
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
|
||||
# Dirichlet else
|
||||
|
||||
"""
|
||||
if(type(BC) is str):
|
||||
BC = [BC for _ in self.vnC] # Repeat the str self.dim times
|
||||
elif(type(BC) is list):
|
||||
assert len(BC) == self.dim, 'BC list must be the size of your mesh'
|
||||
else:
|
||||
raise Exception("BC must be a str or a list.")
|
||||
|
||||
for i, bc_i in enumerate(BC):
|
||||
BC[i] = checkBC(bc_i)
|
||||
|
||||
|
||||
def projDirichlet(n, bc):
|
||||
bc = checkBC(bc)
|
||||
ij = ([0,n], [0,1])
|
||||
vals = [0,0]
|
||||
if(bc[0] == 'dirichlet'):
|
||||
vals[0] = -1
|
||||
if(bc[1] == 'dirichlet'):
|
||||
vals[1] = 1
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1,2))
|
||||
|
||||
def projNeumannIn(n, bc):
|
||||
bc = checkBC(bc)
|
||||
P = sp.identity(n+1).tocsr()
|
||||
if(bc[0] == 'neumann'):
|
||||
P = P[1:,:]
|
||||
if(bc[0] == 'neumann'):
|
||||
P = P[:-1,:]
|
||||
return P
|
||||
|
||||
def projNeumannOut(n, bc):
|
||||
bc = checkBC(bc)
|
||||
ij = ([0, 1],[0, n])
|
||||
vals = [0,0]
|
||||
if(bc[0] == 'neumann'):
|
||||
vals[0] = 1
|
||||
if(bc[1] == 'neumann'):
|
||||
vals[1] = 1
|
||||
return sp.csr_matrix((vals, ij), shape=(2,n+1))
|
||||
|
||||
n = self.vnC
|
||||
indF = self.faceBoundaryInd
|
||||
if(self.dim == 1):
|
||||
Pbc = projDirichlet(n[0], BC[0])
|
||||
indF = indF[0] | indF[1]
|
||||
Pbc = Pbc*sdiag(self.area[indF])
|
||||
|
||||
Pin = projNeumannIn(n[0], BC[0])
|
||||
|
||||
Pout = projNeumannOut(n[0], BC[0])
|
||||
elif(self.dim == 2):
|
||||
Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
|
||||
Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
|
||||
Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
|
||||
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
|
||||
Pbc = Pbc*sdiag(self.area[indF])
|
||||
|
||||
P1 = sp.kron(speye(n[1]), projNeumannIn(n[0], BC[0]))
|
||||
P2 = sp.kron(projNeumannIn(n[1], BC[1]), speye(n[0]))
|
||||
Pin = sp.block_diag((P1, P2), format="csr")
|
||||
|
||||
P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
|
||||
P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
|
||||
Pout = sp.block_diag((P1, P2), format="csr")
|
||||
elif(self.dim == 3):
|
||||
Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
|
||||
Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
|
||||
Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
|
||||
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
|
||||
Pbc = Pbc*sdiag(self.area[indF])
|
||||
|
||||
P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
|
||||
P2 = kron3(speye(n[2]), projNeumannIn(n[1], BC[1]), speye(n[0]))
|
||||
P3 = kron3(projNeumannIn(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
Pin = sp.block_diag((P1, P2, P3), format="csr")
|
||||
|
||||
P1 = kron3(speye(n[2]), speye(n[1]), projNeumannOut(n[0], BC[0]))
|
||||
P2 = kron3(speye(n[2]), projNeumannOut(n[1], BC[1]), speye(n[0]))
|
||||
P3 = kron3(projNeumannOut(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
Pout = sp.block_diag((P1, P2, P3), format="csr")
|
||||
|
||||
return Pbc, Pin, Pout
|
||||
|
||||
|
||||
# --------------- Averaging ---------------------
|
||||
|
||||
@property
|
||||
|
||||
@@ -419,6 +419,56 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
|
||||
return Q
|
||||
|
||||
|
||||
@property
|
||||
def faceBoundaryInd(self):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
"""
|
||||
if self.dim==1:
|
||||
indxd = (self.gridFx==min(self.gridFx))
|
||||
indxu = (self.gridFx==max(self.gridFx))
|
||||
return indxd, indxu
|
||||
elif self.dim==2:
|
||||
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
|
||||
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
|
||||
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
|
||||
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif self.dim==3:
|
||||
indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
|
||||
indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
|
||||
indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
|
||||
indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
|
||||
indzd = (self.gridFz[:,2]==min(self.gridFz[:,2]))
|
||||
indzu = (self.gridFz[:,2]==max(self.gridFz[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
@property
|
||||
def cellBoundaryInd(self):
|
||||
"""
|
||||
Find indices of boundary faces in each direction
|
||||
"""
|
||||
if self.dim==1:
|
||||
indxd = (self.gridCC==min(self.gridCC))
|
||||
indxu = (self.gridCC==max(self.gridCC))
|
||||
return indxd, indxu
|
||||
elif self.dim==2:
|
||||
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
|
||||
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
|
||||
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
|
||||
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
|
||||
return indxd, indxu, indyd, indyu
|
||||
elif self.dim==3:
|
||||
indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
|
||||
indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
|
||||
indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
|
||||
indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
|
||||
indzd = (self.gridCC[:,2]==min(self.gridCC[:,2]))
|
||||
indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
print('Welcome to tensor mesh!')
|
||||
|
||||
|
||||
@@ -74,7 +74,7 @@ class TensorView(object):
|
||||
if I.size != np.prod(self.vnEz): ex, ey, I = self.r(I,'E','E','M')
|
||||
elif imageType[0] == 'E':
|
||||
plotAll = len(imageType) == 1
|
||||
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
|
||||
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
|
||||
fig = plt.figure(figNum)
|
||||
# Determine the subplot number: 131, 121
|
||||
numPlots = 130 if plotAll else len(imageType)/2*10+100
|
||||
@@ -92,10 +92,11 @@ class TensorView(object):
|
||||
ax_z = plt.subplot(numPlots+pltNum)
|
||||
self.plotImage(ez, imageType='Ez', ax=ax_z, **options)
|
||||
pltNum +=1
|
||||
if showIt: plt.show()
|
||||
return
|
||||
elif imageType[0] == 'F':
|
||||
plotAll = len(imageType) == 1
|
||||
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
|
||||
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
|
||||
fig = plt.figure(figNum)
|
||||
# Determine the subplot number: 131, 121
|
||||
numPlots = 130 if plotAll else len(imageType)/2*10+100
|
||||
@@ -113,6 +114,7 @@ class TensorView(object):
|
||||
ax_z = plt.subplot(numPlots+pltNum)
|
||||
self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
|
||||
pltNum +=1
|
||||
if showIt: plt.show()
|
||||
return
|
||||
else:
|
||||
raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'")
|
||||
|
||||
@@ -93,9 +93,9 @@ class OrderTest(unittest.TestCase):
|
||||
h3 = np.ones(nc)/nc
|
||||
h = [h1, h2, h3]
|
||||
elif 'random' in self._meshType:
|
||||
h1 = np.random.rand(nc)
|
||||
h2 = np.random.rand(nc)
|
||||
h3 = np.random.rand(nc)
|
||||
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
|
||||
else:
|
||||
raise Exception('Unexpected meshType')
|
||||
@@ -111,10 +111,12 @@ class OrderTest(unittest.TestCase):
|
||||
kwrd = 'rotate'
|
||||
else:
|
||||
raise Exception('Unexpected meshType')
|
||||
if self.meshDimension == 2:
|
||||
if self.meshDimension == 1:
|
||||
raise Exception('Lom not supported for 1D')
|
||||
elif self.meshDimension == 2:
|
||||
X, Y = Utils.exampleLomGird([nc, nc], kwrd)
|
||||
self.M = LogicallyOrthogonalMesh([X, Y])
|
||||
if self.meshDimension == 3:
|
||||
elif self.meshDimension == 3:
|
||||
X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd)
|
||||
self.M = LogicallyOrthogonalMesh([X, Y, Z])
|
||||
return 1./nc
|
||||
|
||||
@@ -0,0 +1,166 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
from TestUtils import OrderTest
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG import Utils, Solver
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh']
|
||||
|
||||
class Test1D_InhomogeneousDirichlet(OrderTest):
|
||||
name = "1D - Dirichlet"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 1
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.cos(np.pi*x)
|
||||
j_fun = lambda x: -np.pi*np.sin(np.pi*x)
|
||||
q_fun = lambda x: -(np.pi**2)*np.cos(np.pi*x)
|
||||
|
||||
xc_anal = phi(self.M.gridCC)
|
||||
q_anal = q_fun(self.M.gridCC)
|
||||
j_anal = j_fun(self.M.gridFx)
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
# vec = self.M.vectorNx
|
||||
vec = self.M.vectorCCx
|
||||
bc = phi(vec[[0,-1]])
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'dirichlet']])
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol)
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_anal + P*bc)
|
||||
q = D*j
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = D*McI*G
|
||||
rhs = q_anal - D*McI*P*bc
|
||||
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((j-j_anal), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-q_anal), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
xc = Solver(A).solve(rhs)
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
elif self.myTest == 'xcJ':
|
||||
xc = Solver(A).solve(rhs)
|
||||
j = McI*(G*xc + P*bc)
|
||||
err = np.linalg.norm((j-j_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderX(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Inverse"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class Test2D_InhomogeneousDirichlet(OrderTest):
|
||||
name = "2D - Dirichlet"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
j_funX = lambda x: -np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
j_funY = lambda x: -np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
|
||||
q_fun = lambda x: -2*(np.pi**2)*phi(x)
|
||||
|
||||
xc_anal = phi(self.M.gridCC)
|
||||
q_anal = q_fun(self.M.gridCC)
|
||||
jX_anal = j_funX(self.M.gridFx)
|
||||
jY_anal = j_funY(self.M.gridFy)
|
||||
j_anal = np.r_[jX_anal,jY_anal]
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
# fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
# gBFx = self.M.gridFx[(fxm|fxp),:]
|
||||
# gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
fxm,fxp,fym,fyp = self.M.cellBoundaryInd
|
||||
gBFx = self.M.gridCC[(fxm|fxp),:]
|
||||
gBFy = self.M.gridCC[(fym|fyp),:]
|
||||
|
||||
bc = phi(np.r_[gBFx,gBFy])
|
||||
|
||||
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF('dirichlet')
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol)
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_anal + P*bc)
|
||||
q = D*j
|
||||
|
||||
# self.M.plotImage(j, 'FxFy', showIt=True)
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = D*McI*G
|
||||
rhs = q_anal - D*McI*P*bc
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((j-j_anal), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-q_anal), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
xc = Solver(A).solve(rhs)
|
||||
err = np.linalg.norm((xc-xc_anal), np.inf)
|
||||
elif self.myTest == 'xcJ':
|
||||
xc = Solver(A).solve(rhs)
|
||||
j = McI*(G*xc + P*bc)
|
||||
err = np.linalg.norm((j-j_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderX(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Inverse"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user