From 20e80ed983e2992ab369505541b5f752d0ed7a28 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 21 Jul 2016 11:55:14 -0700 Subject: [PATCH] use `@property` decorator in DiffOperators.py --- SimPEG/Mesh/DiffOperators.py | 725 +++++++++++++++++++---------------- 1 file changed, 394 insertions(+), 331 deletions(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index b4a4b1ba..e755372a 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -18,13 +18,15 @@ def checkBC(bc): for bc_i in bc: assert type(bc_i) is str, "each bc must be a string" - assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" + assert bc_i in ['dirichlet', 'neumann'], ("each bc must be either," + "'dirichlet' or 'neumann'") return bc def ddxCellGrad(n, bc): """ - Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + Create 1D derivative operator from cell-centers to nodes this means we + go from n to n+1 For Cell-Centered **Dirichlet**, use a ghost point:: @@ -52,7 +54,8 @@ def ddxCellGrad(n, bc): """ bc = checkBC(bc) - D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") + D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, + format="csr") # Set the first side if(bc[0] == 'dirichlet'): D[0, 0] = 2 @@ -65,10 +68,11 @@ def ddxCellGrad(n, bc): D[-1, -1] = 0 return D + def ddxCellGradBC(n, bc): """ - - Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + Create 1D derivative operator from cell-centers to nodes this means we + go from n to n+1 For Cell-Centered **Dirichlet**, use a ghost point:: @@ -99,7 +103,7 @@ def ddxCellGradBC(n, bc): """ bc = checkBC(bc) - ij = (np.array([0, n]),np.array([0, 1])) + ij = (np.array([0, n]), np.array([0, 1])) vals = np.zeros(2) # Set the first side @@ -112,7 +116,7 @@ def ddxCellGradBC(n, bc): vals[1] = 2 elif(bc[1] == 'neumann'): vals[1] = 0 - D = sp.csr_matrix((vals, ij), shape=(n+1,2)) + D = sp.csr_matrix((vals, ij), shape=(n+1, 2)) return D @@ -121,175 +125,166 @@ class DiffOperators(object): Class creates the differential operators that you need! """ def __init__(self): - raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') + raise Exception('DiffOperators is a base class providing differential' + 'operators on meshes and cannot run on its own.' + 'Inherit to your favorite Mesh class.') - def faceDiv(): - doc = "Construct divergence operator (face-stg to cell-centres)." - - def fget(self): - if(self._faceDiv is None): - # The number of cell centers in each direction - n = self.vnC - # Compute faceDivergence operator on faces - if(self.dim == 1): - D = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]), ddx(n[0])) - D2 = sp.kron(ddx(n[1]), speye(n[0])) - D = sp.hstack((D1, D2), format="csr") - elif(self.dim == 3): - D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) - D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) - D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) - D = sp.hstack((D1, D2, D3), format="csr") - # Compute areas of cell faces & volumes - S = self.area - V = self.vol - self._faceDiv = sdiag(1/V)*D*sdiag(S) - - return self._faceDiv - return locals() - _faceDiv = None - faceDiv = property(**faceDiv()) - - def faceDivx(): - doc = "Construct divergence operator in the x component (face-stg to cell-centres)." - - def fget(self): - if(self._faceDivx is None): - # The number of cell centers in each direction - n = self.vnC - # Compute faceDivergence operator on faces - if(self.dim == 1): - D1 = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]), ddx(n[0])) - elif(self.dim == 3): - D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) - # Compute areas of cell faces & volumes - S = self.r(self.area, 'F', 'Fx', 'V') - V = self.vol - self._faceDivx = sdiag(1/V)*D1*sdiag(S) - - return self._faceDivx - return locals() - _faceDivx = None - faceDivx = property(**faceDivx()) - - def faceDivy(): - doc = "Construct divergence operator in the y component (face-stg to cell-centres)." - - def fget(self): - if(self.dim < 2): return None - if(self._faceDivy is None): - # The number of cell centers in each direction - n = self.vnC - # Compute faceDivergence operator on faces - if(self.dim == 2): - D2 = sp.kron(ddx(n[1]), speye(n[0])) - elif(self.dim == 3): - D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) - # Compute areas of cell faces & volumes - S = self.r(self.area, 'F', 'Fy', 'V') - V = self.vol - self._faceDivy = sdiag(1/V)*D2*sdiag(S) - - return self._faceDivy - return locals() - _faceDivy = None - faceDivy = property(**faceDivy()) - - def faceDivz(): - doc = "Construct divergence operator in the z component (face-stg to cell-centres)." - - def fget(self): - if(self.dim < 3): return None - if(self._faceDivz is None): - # The number of cell centers in each direction - n = self.vnC - # Compute faceDivergence operator on faces + @property + def faceDiv(self): + """ + Construct divergence operator (face-stg to cell-centres). + """ + if getattr(self, '_faceDiv', None) is None: + n = self.vnC + # Compute faceDivergence operator on faces + if(self.dim == 1): + D = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0])) + D = sp.hstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) + D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) - # Compute areas of cell faces & volumes - S = self.r(self.area, 'F', 'Fz', 'V') - V = self.vol - self._faceDivz = sdiag(1/V)*D3*sdiag(S) + D = sp.hstack((D1, D2, D3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.vol + self._faceDiv = sdiag(1/V)*D*sdiag(S) + return self._faceDiv - return self._faceDivz - return locals() - _faceDivz = None - faceDivz = property(**faceDivz()) + @property + def faceDivx(self): + """ + Construct divergence operator in the x component (face-stg to + cell-centres). + """ + if getattr(self, '_faceDivx', None) is None: + # The number of cell centers in each direction + n = self.vnC + # Compute faceDivergence operator on faces + if(self.dim == 1): + D1 = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]), ddx(n[0])) + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F', 'Fx', 'V') + V = self.vol + self._faceDivx = sdiag(1/V)*D1*sdiag(S) - def nodalGrad(): - doc = "Construct gradient operator (nodes to edges)." + return self._faceDivx - def fget(self): - if(self._nodalGrad is None): - # The number of cell centers in each direction - n = self.vnC - # Compute divergence operator on faces - if(self.dim == 1): - G = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]+1), ddx(n[0])) - D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) - G = sp.vstack((D1, D2), format="csr") - elif(self.dim == 3): - D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) - D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) - D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) - G = sp.vstack((D1, D2, D3), format="csr") - # Compute lengths of cell edges - L = self.edge - self._nodalGrad = sdiag(1/L)*G - return self._nodalGrad - return locals() - _nodalGrad = None - nodalGrad = property(**nodalGrad()) + @property + def faceDivy(self): + if(self.dim < 2): + return None + if getattr(self, '_faceDivy', None) is None: + # The number of cell centers in each direction + n = self.vnC + # Compute faceDivergence operator on faces + if(self.dim == 2): + D2 = sp.kron(ddx(n[1]), speye(n[0])) + elif(self.dim == 3): + D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F', 'Fy', 'V') + V = self.vol + self._faceDivy = sdiag(1/V)*D2*sdiag(S) + return self._faceDivy - def nodalLaplacian(): - doc = "Construct laplacian operator (nodes to edges)." + @property + def faceDivz(self): + """ + Construct divergence operator in the z component (face-stg to + cell-centres). + """ + if(self.dim < 3): + return None + if(self._faceDivz is None): + # The number of cell centers in each direction + n = self.vnC + # Compute faceDivergence operator on faces + D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F', 'Fz', 'V') + V = self.vol + self._faceDivz = sdiag(1/V)*D3*sdiag(S) + return self._faceDivz - def fget(self): - if(self._nodalLaplacian is None): - print 'Warning: Laplacian has not been tested rigorously.' - # The number of cell centers in each direction - n = self.vnC - # Compute divergence operator on faces - if(self.dim == 1): - D1 = sdiag(1./self.hx) * ddx(mesh.nCx) - L = - D1.T*D1 - elif(self.dim == 2): - D1 = sdiag(1./self.hx) * ddx(n[0]) - D2 = sdiag(1./self.hy) * ddx(n[1]) - L1 = sp.kron(speye(n[1]+1), - D1.T * D1) - L2 = sp.kron(- D2.T * D2, speye(n[0]+1)) - L = L1 + L2 - elif(self.dim == 3): - D1 = sdiag(1./self.hx) * ddx(n[0]) - D2 = sdiag(1./self.hy) * ddx(n[1]) - D3 = sdiag(1./self.hz) * ddx(n[2]) - L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1) - L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1)) - L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1)) - L = L1 + L2 + L3 - self._nodalLaplacian = L - return self._nodalLaplacian - return locals() - _nodalLaplacian = None - nodalLaplacian = property(**nodalLaplacian()) + @property + def nodalGrad(self): + """ + Construct gradient operator (nodes to edges). + """ + if getattr(self, '_nodalGrad', None) is None: + # The number of cell centers in each direction + n = self.vnC + # Compute divergence operator on faces + if(self.dim == 1): + G = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]+1), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) + G = sp.vstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) + D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) + D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) + G = sp.vstack((D1, D2, D3), format="csr") + # Compute lengths of cell edges + L = self.edge + self._nodalGrad = sdiag(1/L)*G + return self._nodalGrad + + @property + def nodalLaplacian(self): + """ + Construct laplacian operator (nodes to edges). + """ + if getattr(self, '_nodalLaplacian', None) is None: + print 'Warning: Laplacian has not been tested rigorously.' + # The number of cell centers in each direction + n = self.vnC + # Compute divergence operator on faces + if(self.dim == 1): + D1 = sdiag(1./self.hx) * ddx(mesh.nCx) + L = - D1.T*D1 + elif(self.dim == 2): + D1 = sdiag(1./self.hx) * ddx(n[0]) + D2 = sdiag(1./self.hy) * ddx(n[1]) + L1 = sp.kron(speye(n[1]+1), - D1.T * D1) + L2 = sp.kron(- D2.T * D2, speye(n[0]+1)) + L = L1 + L2 + elif(self.dim == 3): + D1 = sdiag(1./self.hx) * ddx(n[0]) + D2 = sdiag(1./self.hy) * ddx(n[1]) + D3 = sdiag(1./self.hz) * ddx(n[2]) + L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1) + L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1)) + L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1)) + L = L1 + L2 + L3 + self._nodalLaplacian = L + return self._nodalLaplacian def setCellGradBC(self, BC): """ - Function that sets the boundary conditions for cell-centred derivative operators. + Function that sets the boundary conditions for cell-centred derivative + operators. Examples:: + # Neumann in all directions + BC = 'neumann' - 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 + # 3D, Dirichlet in y Neumann else + BC = ['neumann', 'dirichlet', 'neumann'] + # 3D, Neumann in x on bottom of domain, Dirichlet else + BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] """ + if(type(BC) is str): BC = [BC]*self.dim if(type(BC) is list): @@ -323,47 +318,69 @@ class DiffOperators(object): G = sp.vstack((G1, G2, G3), format="csr") return G - def cellGrad(): - doc = "The cell centered Gradient, takes you to cell faces." + @property + def cellGrad(self): + """ + The cell centered Gradient, takes you to cell faces. + """ + if(self._cellGrad is None): + G = self._cellGradStencil() + S = self.area # Compute areas of cell faces & volumes + V = self.aveCC2F*self.vol # Average volume between adjacent cells + self._cellGrad = sdiag(S/V)*G + return self._cellGrad - def fget(self): - if(self._cellGrad is None): - G = self._cellGradStencil() - # Compute areas of cell faces & volumes - S = self.area - V = self.aveCC2F*self.vol # Average volume between adjacent cells - self._cellGrad = sdiag(S/V)*G - return self._cellGrad - return locals() - _cellGrad = None - cellGrad = property(**cellGrad()) + @property + def cellGradBC(self): + """ + The cell centered Gradient boundary condition matrix + """ + if getattr(self, '_cellGradBC', None) is None: + BC = self.setCellGradBC(self._cellGradBC_list) + n = self.vnC + if(self.dim == 1): + G = ddxCellGradBC(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0])) + G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0])) + G = sp.block_diag((G1, G2), format="csr") + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0])) + G = sp.block_diag((G1, G2, G3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.aveCC2F*self.vol # Average volume between adjacent cells + self._cellGradBC = sdiag(S/V)*G + return self._cellGradBC - def cellGradBC(): - doc = "The cell centered Gradient boundary condition matrix" + # def cellGradBC(): + # doc = "The cell centered Gradient boundary condition matrix" - def fget(self): - if(self._cellGradBC is None): - BC = self.setCellGradBC(self._cellGradBC_list) - n = self.vnC - if(self.dim == 1): - G = ddxCellGradBC(n[0], BC[0]) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0])) - G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0])) - G = sp.block_diag((G1, G2), format="csr") - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0])) - G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0])) - G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0])) - G = sp.block_diag((G1, G2, G3), format="csr") - # Compute areas of cell faces & volumes - S = self.area - V = self.aveCC2F*self.vol # Average volume between adjacent cells - self._cellGradBC = sdiag(S/V)*G - return self._cellGradBC - return locals() - _cellGradBC = None - cellGradBC = property(**cellGradBC()) + # def fget(self): + # if(self._cellGradBC is None): + # BC = self.setCellGradBC(self._cellGradBC_list) + # n = self.vnC + # if(self.dim == 1): + # G = ddxCellGradBC(n[0], BC[0]) + # elif(self.dim == 2): + # G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0])) + # G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0])) + # G = sp.block_diag((G1, G2), format="csr") + # elif(self.dim == 3): + # G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0])) + # G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0])) + # G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0])) + # G = sp.block_diag((G1, G2, G3), format="csr") + # # Compute areas of cell faces & volumes + # S = self.area + # V = self.aveCC2F*self.vol # Average volume between adjacent cells + # self._cellGradBC = sdiag(S/V)*G + # return self._cellGradBC + # return locals() + # _cellGradBC = None + # cellGradBC = property(**cellGradBC()) def _cellGradxStencil(self): BC = ['neumann', 'neumann'] @@ -376,20 +393,19 @@ class DiffOperators(object): G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) return G1 - - def cellGradx(): - doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." - - def fget(self): - if getattr(self, '_cellGradx', None) is None: - G1 = self._cellGradxStencil() - # Compute areas of cell faces & volumes - V = self.aveCC2F*self.vol - L = self.r(self.area/V, 'F','Fx', 'V') - self._cellGradx = sdiag(L)*G1 - return self._cellGradx - return locals() - cellGradx = property(**cellGradx()) + @property + def cellGradx(self): + """ + Cell centered Gradient in the x dimension. Has neumann boundary + conditions. + """ + if getattr(self, '_cellGradx', None) is None: + G1 = self._cellGradxStencil() + # Compute areas of cell faces & volumes + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fx', 'V') + self._cellGradx = sdiag(L)*G1 + return self._cellGradx def _cellGradyStencil(self): if self.dim < 2: return None @@ -401,19 +417,17 @@ class DiffOperators(object): G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) return G2 - def cellGrady(): - doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." - def fget(self): - if self.dim < 2: return None - if getattr(self, '_cellGrady', None) is None: - G2 = self._cellGradyStencil() - # Compute areas of cell faces & volumes - V = self.aveCC2F*self.vol - L = self.r(self.area/V, 'F','Fy', 'V') - self._cellGrady = sdiag(L)*G2 - return self._cellGrady - return locals() - cellGrady = property(**cellGrady()) + @property + def cellGrady(self): + if self.dim < 2: + return None + if getattr(self, '_cellGrady', None) is None: + G2 = self._cellGradyStencil() + # Compute areas of cell faces & volumes + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F', 'Fy', 'V') + self._cellGrady = sdiag(L)*G2 + return self._cellGrady def _cellGradzStencil(self): if self.dim < 3: return None @@ -422,66 +436,61 @@ class DiffOperators(object): G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) return G3 - def cellGradz(): - doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." - def fget(self): - if self.dim < 3: return None - if getattr(self, '_cellGradz', None) is None: - G3 = self._cellGradzStencil() - # Compute areas of cell faces & volumes - V = self.aveCC2F*self.vol - L = self.r(self.area/V, 'F','Fz', 'V') - self._cellGradz = sdiag(L)*G3 - return self._cellGradz - return locals() - cellGradz = property(**cellGradz()) + @property + def cellGradz(self): + """ + Cell centered Gradient in the x dimension. Has neumann boundary + conditions. + """ + if self.dim < 3: + return None + if getattr(self, '_cellGradz', None) is None: + G3 = self._cellGradzStencil() + # Compute areas of cell faces & volumes + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F', 'Fz', 'V') + self._cellGradz = sdiag(L)*G3 + return self._cellGradz - def edgeCurl(): - doc = "Construct the 3D curl operator." + @property + def edgeCurl(self): + """ + Construct the 3D curl operator. + """ + if getattr(self, '_edgeCurl', None) is None: + assert self.dim > 1, "Edge Curl only programed for 2 or 3D." - def fget(self): - if(self._edgeCurl is None): - assert self.dim > 1, "Edge Curl only programed for 2 or 3D." - # The number of cell centers in each direction - n = self.vnC + n = self.vnC # The number of cell centers in each direction + L = self.edge # Compute lengths of cell edges + S = self.area # Compute areas of cell faces - # Compute lengths of cell edges - L = self.edge + # Compute divergence operator on faces + if self.dim == 2: - # Compute areas of cell faces - S = self.area + D21 = sp.kron(ddx(n[1]), speye(n[0])) + D12 = sp.kron(speye(n[1]), ddx(n[0])) + C = sp.hstack((-D21, D12), format="csr") + self._edgeCurl = C*sdiag(1/S) - # Compute divergence operator on faces - if self.dim == 2: + elif self.dim == 3: - D21 = sp.kron(ddx(n[1]), speye(n[0])) - D12 = sp.kron(speye(n[1]), ddx(n[0])) - C = sp.hstack((-D21, D12), format="csr") - self._edgeCurl = C*sdiag(1/S) + D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1)) + D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1)) + D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0])) + D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0])) + D21 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0])) + D12 = kron3(speye(n[2]+1), speye(n[1]), ddx(n[0])) - elif self.dim == 3: + O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) + O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) + O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1)) - D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1)) - D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0])) - D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0])) - D21 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0])) - D12 = kron3(speye(n[2]+1), speye(n[1]), ddx(n[0])) + C = sp.vstack((sp.hstack((O1, -D32, D23)), + sp.hstack((D31, O2, -D13)), + sp.hstack((-D21, D12, O3))), format="csr") - O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) - O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) - O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - - C = sp.vstack((sp.hstack((O1, -D32, D23)), - sp.hstack((D31, O2, -D13)), - sp.hstack((-D21, D12, O3))), format="csr") - - self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) - - return self._edgeCurl - return locals() - _edgeCurl = None - edgeCurl = property(**edgeCurl()) + self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) + return self._edgeCurl def getBCProjWF(self, BC, discretization='CC'): """ @@ -489,16 +498,19 @@ class DiffOperators(object): The weak form boundary condition projection matrices. Examples:: + # Neumann in all directions + BC = 'neumann' - 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 + # 3D, Dirichlet in y Neumann else + BC = ['neumann', 'dirichlet', 'neumann'] + # 3D, Neumann in x on bottom of domain, Dirichlet else + BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] """ if discretization is not 'CC': - raise NotImplementedError('Boundary conditions only implemented for CC discretization.') + raise NotImplementedError('Boundary conditions only implemented' + 'for CC discretization.') if(type(BC) is str): BC = [BC for _ in self.vnC] # Repeat the str self.dim times @@ -510,35 +522,34 @@ class DiffOperators(object): 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] + 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)) + 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:,:] + P = P[1:, :] if(bc[1] == 'neumann'): - P = P[:-1,:] + P = P[:-1, :] return P def projNeumannOut(n, bc): bc = checkBC(bc) - ij = ([0, 1],[0, n]) + 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)) + return sp.csr_matrix((vals, ij), shape=(2, n+1)) n = self.vnC indF = self.faceBoundaryInd @@ -550,6 +561,7 @@ class DiffOperators(object): 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])) @@ -564,12 +576,14 @@ class DiffOperators(object): 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])] + 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])) @@ -586,36 +600,36 @@ class DiffOperators(object): def getBCProjWF_simple(self, discretization='CC'): """ - The weak form boundary condition projection matrices when mixed boundary condition is used - - """ if discretization is not 'CC': - raise NotImplementedError('Boundary conditions only implemented for CC discretization.') + raise NotImplementedError('Boundary conditions only implemented' + 'for CC discretization.') def projBC(n): - ij = ([0,n], [0,1]) - vals = [0,0] + ij = ([0, n], [0, 1]) + vals = [0, 0] vals[0] = 1 vals[1] = 1 - return sp.csr_matrix((vals, ij), shape=(n+1,2)) + return sp.csr_matrix((vals, ij), shape=(n+1, 2)) def projDirichlet(n, bc): bc = checkBC(bc) - ij = ([0,n], [0,1]) - vals = [0,0] + 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)) + return sp.csr_matrix((vals, ij), shape=(n+1, 2)) - BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']] + BC = [['dirichlet', 'dirichlet'], ['dirichlet', 'dirichlet'], + ['dirichlet', 'dirichlet']] n = self.vnC indF = self.faceBoundaryInd + if(self.dim == 1): Pbc = projDirichlet(n[0], BC[0]) B = projBC(n[0]) @@ -653,9 +667,11 @@ class DiffOperators(object): if(self.dim == 1): return self.aveFx2CC elif(self.dim == 2): - return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr") + return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), + format="csr") elif(self.dim == 3): - return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr") + return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, + self.aveFz2CC), format="csr") @property def aveF2CCV(self): @@ -665,11 +681,16 @@ class DiffOperators(object): elif(self.dim == 2): return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr") elif(self.dim == 3): - return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr") + return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), + format="csr") @property def aveFx2CC(self): - "Construct the averaging operator on cell faces in the x direction to cell centers." + """ + Construct the averaging operator on cell faces in the x direction to + cell centers. + """ + if getattr(self, '_aveFx2CC', None) is None: n = self.vnC if(self.dim == 1): @@ -682,8 +703,12 @@ class DiffOperators(object): @property def aveFy2CC(self): - "Construct the averaging operator on cell faces in the y direction to cell centers." - if self.dim < 2: return None + """ + Construct the averaging operator on cell faces in the y direction to + cell centers. + """ + if self.dim < 2: + return None if getattr(self, '_aveFy2CC', None) is None: n = self.vnC if(self.dim == 2): @@ -694,7 +719,10 @@ class DiffOperators(object): @property def aveFz2CC(self): - "Construct the averaging operator on cell faces in the z direction to cell centers." + """ + Construct the averaging operator on cell faces in the z direction to + cell centers. + """ if self.dim < 3: return None if getattr(self, '_aveFz2CC', None) is None: n = self.vnC @@ -711,12 +739,18 @@ class DiffOperators(object): if(self.dim == 1): self._aveCC2F = avExtrap(n[0]) elif(self.dim == 2): - self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])), - sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr") + self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), + avExtrap(n[0])), + sp.kron(avExtrap(n[1]), + speye(n[0]))), format="csr") elif(self.dim == 3): - self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])), - kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])), - kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr") + self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), + avExtrap(n[0])), + kron3(speye(n[2]), avExtrap(n[1]), + speye(n[0])), + kron3(avExtrap(n[2]), speye(n[1]), + speye(n[0]))), + format="csr") return self._aveCC2F @property @@ -727,7 +761,8 @@ class DiffOperators(object): elif(self.dim == 2): return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr") elif(self.dim == 3): - return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr") + return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, + self.aveEz2CC), format="csr") @property def aveE2CCV(self): @@ -737,11 +772,15 @@ class DiffOperators(object): elif(self.dim == 2): return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr") elif(self.dim == 3): - return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr") + return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), + format="csr") @property def aveEx2CC(self): - "Construct the averaging operator on cell edges in the x direction to cell centers." + """ + Construct the averaging operator on cell edges in the x direction to + cell centers. + """ if getattr(self, '_aveEx2CC', None) is None: # The number of cell centers in each direction n = self.vnC @@ -755,8 +794,12 @@ class DiffOperators(object): @property def aveEy2CC(self): - "Construct the averaging operator on cell edges in the y direction to cell centers." - if self.dim < 2: return None + """ + Construct the averaging operator on cell edges in the y direction to + cell centers. + """ + if self.dim < 2: + return None if getattr(self, '_aveEy2CC', None) is None: # The number of cell centers in each direction n = self.vnC @@ -768,8 +811,12 @@ class DiffOperators(object): @property def aveEz2CC(self): - "Construct the averaging operator on cell edges in the z direction to cell centers." - if self.dim < 3: return None + """ + Construct the averaging operator on cell edges in the z direction to + cell centers. + """ + if self.dim < 3: + return None if getattr(self, '_aveEz2CC', None) is None: # The number of cell centers in each direction n = self.vnC @@ -793,7 +840,10 @@ class DiffOperators(object): @property def aveN2E(self): - "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate." + """ + Construct the averaging operator on cell nodes to cell edges, keeping + each dimension separate. + """ if getattr(self, '_aveN2E', None) is None: # The number of cell centers in each direction @@ -802,16 +852,24 @@ class DiffOperators(object): self._aveN2E = av(n[0]) elif(self.dim == 2): self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])), - sp.kron(av(n[1]), speye(n[0]+1))), format="csr") + sp.kron(av(n[1]), speye(n[0]+1))), + format="csr") elif(self.dim == 3): - self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])), - kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)), - kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr") + self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), + av(n[0])), + kron3(speye(n[2]+1), av(n[1]), + speye(n[0]+1)), + kron3(av(n[2]), speye(n[1]+1), + speye(n[0]+1))), + format="csr") return self._aveN2E @property def aveN2F(self): - "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate." + """ + Construct the averaging operator on cell nodes to cell faces, keeping + each dimension separate. + """ if getattr(self, '_aveN2F', None) is None: # The number of cell centers in each direction n = self.vnC @@ -819,9 +877,14 @@ class DiffOperators(object): self._aveN2F = av(n[0]) elif(self.dim == 2): self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)), - sp.kron(speye(n[1]+1), av(n[0]))), format="csr") + sp.kron(speye(n[1]+1), av(n[0]))), + format="csr") elif(self.dim == 3): - self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)), - kron3(av(n[2]), speye(n[1]+1), av(n[0])), - kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") + self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), + speye(n[0]+1)), + kron3(av(n[2]), speye(n[1]+1), + av(n[0])), + kron3(speye(n[2]+1), av(n[1]), + av(n[0]))), + format="csr") return self._aveN2F