From fbdde2ac245a5b002c67064065472f3d28bb557e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Feb 2014 12:00:50 -0800 Subject: [PATCH] boundaryCondition initial work. --- SimPEG/Mesh/DiffOperators.py | 98 ++++++++++++++++ SimPEG/Mesh/TensorMesh.py | 50 ++++++++ SimPEG/Mesh/TensorView.py | 6 +- SimPEG/Tests/TestUtils.py | 12 +- SimPEG/Tests/test_boundaryPoisson.py | 166 +++++++++++++++++++++++++++ 5 files changed, 325 insertions(+), 7 deletions(-) create mode 100644 SimPEG/Tests/test_boundaryPoisson.py diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 1712964e..8ca498d7 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -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 diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 12253bb5..8f83e793 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -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!') diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index e4dd9243..597bf989 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -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'") diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index d84408d4..b46fbf10 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -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 diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py new file mode 100644 index 00000000..fd6ede5d --- /dev/null +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -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()