diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index e46bed20..00bf0259 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -571,35 +571,58 @@ class DiffOperators(object): @property def aveF2CC(self): "Construct the averaging operator on cell faces to cell centers." - if getattr(self, '_aveF2CC', None) is None: - n = self.vnC - if(self.dim == 1): - self._aveF2CC = av(n[0]) - elif(self.dim == 2): - self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") - elif(self.dim == 3): - self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") - return self._aveF2CC - + if(self.dim == 1): + return self.aveFx2CC + elif(self.dim == 2): + 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") @property def aveF2CCV(self): "Construct the averaging operator on cell faces to cell centers." - if getattr(self, '_aveF2CCV', None) is None: + if(self.dim == 1): + return self.aveFx2CC + 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") + + @property + def aveFx2CC(self): + "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): - self._aveF2CCV = av(n[0]) + self._aveFx2CC = av(n[0]) elif(self.dim == 2): - self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") + self._aveFx2CC = sp.kron(speye(n[1]), av(n[0])) elif(self.dim == 3): - self._aveF2CCV = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") - return self._aveF2CCV + self._aveFx2CC = kron3(speye(n[2]), speye(n[1]), av(n[0])) + return self._aveFx2CC + + @property + def aveFy2CC(self): + "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): + self._aveFy2CC = sp.kron(av(n[1]), speye(n[0])) + elif(self.dim == 3): + self._aveFy2CC = kron3(speye(n[2]), av(n[1]), speye(n[0])) + return self._aveFy2CC + + @property + def aveFz2CC(self): + "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 + if(self.dim == 3): + self._aveFz2CC = kron3(av(n[2]), speye(n[1]), speye(n[0])) + return self._aveFz2CC + @property def aveCC2F(self): @@ -620,36 +643,60 @@ class DiffOperators(object): @property def aveE2CC(self): "Construct the averaging operator on cell edges to cell centers." - if getattr(self, '_aveE2CC', None) is None: - # The number of cell centers in each direction - n = self.vnC - if(self.dim == 1): - self._aveE2CC = speye(n[0]) - elif(self.dim == 2): - self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveE2CC + if(self.dim == 1): + return self.aveEx2CC + 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") @property def aveE2CCV(self): "Construct the averaging operator on cell edges to cell centers." - if getattr(self, '_aveE2CCV', None) is None: + if(self.dim == 1): + return self.aveEx2CC + 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") + + @property + def aveEx2CC(self): + "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 if(self.dim == 1): - raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') + self._aveEx2CC = speye(n[0]) elif(self.dim == 2): - self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") + self._aveEx2CC = sp.kron(av(n[1]), speye(n[0])) elif(self.dim == 3): - self._aveE2CCV = sp.block_diag((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveE2CCV + self._aveEx2CC = kron3(av(n[2]), av(n[1]), speye(n[0])) + return self._aveEx2CC + + @property + def aveEy2CC(self): + "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 + if(self.dim == 2): + self._aveEy2CC = sp.kron(speye(n[1]), av(n[0])) + elif(self.dim == 3): + self._aveEy2CC = kron3(av(n[2]), speye(n[1]), av(n[0])) + return self._aveEy2CC + + @property + def aveEz2CC(self): + "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 + if(self.dim == 3): + self._aveEz2CC = kron3(speye(n[2]), av(n[1]), av(n[0])) + return self._aveEz2CC @property def aveN2CC(self):