From dcd4fbf97311762376c3fe2698cb980580d0f086 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sun, 24 Apr 2016 15:23:14 -0700 Subject: [PATCH] Implemented mixed B.C. to CC problem. Fix bugs in get fuction getxBCyBC_CC --- SimPEG/EM/Static/DC/BoundaryUtils.py | 36 +++++----- SimPEG/EM/Static/DC/ProblemDC.py | 84 ++++++++++++++++++++---- SimPEG/EM/Static/DC/SrcDC.py | 2 +- tests/em/static/test_DC.py | 9 +-- tests/mesh/test_Mixed_boundaryPoisson.py | 63 +++++++++++------- 5 files changed, 135 insertions(+), 59 deletions(-) diff --git a/SimPEG/EM/Static/DC/BoundaryUtils.py b/SimPEG/EM/Static/DC/BoundaryUtils.py index 5ec3282c..3967eb46 100644 --- a/SimPEG/EM/Static/DC/BoundaryUtils.py +++ b/SimPEG/EM/Static/DC/BoundaryUtils.py @@ -28,7 +28,8 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_xm, h_xp = mesh.hx[0], mesh.hx[-1] a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -47,19 +48,19 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): raise Exception("Lenght of list, alpha should be 4") - fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd fxm,fxp,fym,fyp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -94,23 +95,24 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): elif mesh.dim == 3: #3D if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): raise Exception("Lenght of list, alpha should be 6") - fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] - h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] - alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4] + alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] - h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 7ec32919..9f05d876 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -115,12 +115,7 @@ class Problem3D_CC(BaseDCProblem): def __init__(self, mesh, **kwargs): BaseDCProblem.__init__(self, mesh, **kwargs) - - def setBC(self): - self.Div = V * self.mesh.faceDiv - P_BC, B = self.mesh.getBCProjWF_simple() - M = B*self.mesh.aveCC2F - Grad = Div.T - P_BC*Utils.sdiag(y_BC)*M + self.setBC() def getA(self): """ @@ -131,11 +126,11 @@ class Problem3D_CC(BaseDCProblem): """ - V = self.Vol - D = V * self.mesh.faceDiv + D = self.Div + G = self.Grad # TODO: this won't work for full anisotropy MfRhoI = self.MfRhoI - A = D * MfRhoI * D.T + A = D * MfRhoI * G # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: @@ -144,19 +139,19 @@ class Problem3D_CC(BaseDCProblem): def getADeriv(self, u, v, adjoint= False): - V = self.Vol - D = V * self.mesh.faceDiv + D = self.Div + G = self.Grad MfRhoIDeriv = self.MfRhoIDeriv if adjoint: # if self._makeASymmetric is True: # v = V * v - return(MfRhoIDeriv( D.T * u ).T) * ( D.T * v) + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: # return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) ) - return D * (MfRhoIDeriv( D.T * u ) * v) + return D * (MfRhoIDeriv( G * u ) * v) def getRHS(self): """ @@ -182,6 +177,69 @@ class Problem3D_CC(BaseDCProblem): # return qDeriv return Zero() + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + class Problem3D_N(BaseDCProblem): diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 3ed6067d..1e64835e 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -15,7 +15,7 @@ class BaseSrc(SimPEG.Survey.BaseSrc): raise NotImplementedError def evalDeriv(self, prob): - Zero() + return Zero() class Dipole(BaseSrc): diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py index 99b384aa..227d4614 100644 --- a/tests/em/static/test_DC.py +++ b/tests/em/static/test_DC.py @@ -8,7 +8,7 @@ class DCProblemTestsCC(unittest.TestCase): def setUp(self): aSpacing=2.5 - nElecs=10 + nElecs=5 surveySize = nElecs*aSpacing - aSpacing cs = surveySize/nElecs/4 @@ -16,7 +16,7 @@ class DCProblemTestsCC(unittest.TestCase): mesh = Mesh.TensorMesh([ [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], [(cs,3, -1.3),(cs,3,1.3)], - # [(cs,5, -1.3),(cs,10)] + # [(cs,5, -1.3),(cs,10)] ],'CN') srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) @@ -44,7 +44,7 @@ class DCProblemTestsCC(unittest.TestCase): def test_misfit(self): derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) def test_adjoint(self): @@ -60,7 +60,7 @@ class DCProblemTestsCC(unittest.TestCase): def test_dataObj(self): derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] - passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) class DCProblemTestsN(unittest.TestCase): @@ -122,5 +122,6 @@ class DCProblemTestsN(unittest.TestCase): derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) + if __name__ == '__main__': unittest.main() diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py index d5dcb6fa..3aa1dbcd 100644 --- a/tests/mesh/test_Mixed_boundaryPoisson.py +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -6,10 +6,23 @@ from SimPEG import * MESHTYPES = ['uniformTensorMesh'] -def getxBCyBC(mesh, alpha, beta, gamma): +def getxBCyBC_CC(mesh, alpha, beta, gamma): # def getxBCyBC(mesh, alpha, beta, gamma): """ + This is a subfunction generating mixed-boundary condition: + .. math:: + + \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q + + \rho \vec{j} = -\nabla \phi \phi + + \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega + + xBC = f_1(\alpha, \beta, \gamma) + yBC = f(\alpha, \beta, \gamma) + + Computes xBC and yBC for cell-centered discretizations """ if mesh.dim == 1: #1D if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): @@ -21,7 +34,8 @@ def getxBCyBC(mesh, alpha, beta, gamma): alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_xm, h_xp = mesh.hx[0], mesh.hx[-1] a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -40,19 +54,19 @@ def getxBCyBC(mesh, alpha, beta, gamma): if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): raise Exception("Lenght of list, alpha should be 4") - fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd fxm,fxp,fym,fyp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -87,23 +101,24 @@ def getxBCyBC(mesh, alpha, beta, gamma): elif mesh.dim == 3: #3D if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): raise Exception("Lenght of list, alpha should be 6") - fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd - nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() - h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] - h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] - h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] - alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] - alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4] + alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5] - h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] - h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] - h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) @@ -182,7 +197,7 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): phi_bc = phi_fun(vecN[[0,-1]]) phi_deriv_bc = phi_deriv(vecN[[0,-1]]) gamma = alpha*phi_bc + beta*phi_deriv_bc - x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) sigma = np.ones(self.M.nC) @@ -265,7 +280,7 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest): beta = [beta_xm, beta_xp, beta_ym, beta_yp] gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] - x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) sigma = np.ones(self.M.nC) @@ -335,8 +350,8 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) - alpha_zm, alpha_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) - beta_zm, beta_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) + alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) @@ -345,7 +360,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) - phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funY(gBFzm), phideriv_funY(gBFzp) + phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp) gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) @@ -359,7 +374,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] - x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) sigma = np.ones(self.M.nC)