From 26f9b83b76c81e9a000286f844b840f54d7c7f51 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 18 Oct 2013 14:21:12 -0700 Subject: [PATCH 1/6] Added general getMass method to reduce code duplication. --- SimPEG/DiffOperators.py | 43 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 3b5d8333..666bd98a 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -286,6 +286,49 @@ class DiffOperators(object): _nodalVectorAve = None nodalVectorAve = property(**nodalVectorAve()) + def getMass(self, loc='e', materialProp=None, inv=False): + """ Produces mass matricies. + + Kwargs: + loc (str): 'e' - Average to edges + 'f' faces + materialProp: property to be averaged (see below) + inv (bool): True returns matrix inverse + + Returns: + scipy.sparse.csr.csr_matrix + + materialProp can be: + None -> takes materialProp = 1 (default) + float -> a constant value for entire domain + numpy.ndarray -> if materialProp.size == self.nC + 3D property model + if materialProp.size = self.nCz + 1D (layered eath) property model + """ + if materialProp is None: + materialProp = np.ones(self.nC) + elif type(materialProp) is float: + materialProp = np.ones(self.nC)*materialProp + elif materialProp.shape == (self.nCz,): + materialProp = materialProp.repeat(self.nCx*self.nCy) + materialProp = mkvc(materialProp) + assert materialProp.shape == (self.nC,), "materialProp incorrect shape" + + if loc=='e': + Av = self.edgeAve + elif loc=='f': + Av = self.faceAve + else: + raise ValueError('Invalid loc') + + diag = Av.T * (self.vol * mkvc(materialProp)) + + if inv: + diag = 1/diag + + return sdiag(diag) + def getEdgeMass(self, materialProp=None): """mass matrix for products of edge functions w'*M(materialProp)*e""" if(materialProp is None): From 1f770dfd63482dd21c7385b1949fc16ddbc0c360 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 18 Oct 2013 14:23:57 -0700 Subject: [PATCH 2/6] Updated getEdgeMass and getFaceMass to use the getMass method. --- SimPEG/DiffOperators.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 666bd98a..86fd35f7 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -329,19 +329,13 @@ class DiffOperators(object): return sdiag(diag) - def getEdgeMass(self, materialProp=None): + def getEdgeMass(self, materialProp=None, inv=False): """mass matrix for products of edge functions w'*M(materialProp)*e""" - if(materialProp is None): - materialProp = np.ones(self.nC) - Av = self.edgeAve - return sdiag(Av.T * (self.vol * mkvc(materialProp))) + return self.getMass(loc='e', materialProp=materialProp, inv=inv) - def getFaceMass(self, materialProp=None): + def getFaceMass(self, materialProp=None, inv=False): """mass matrix for products of face functions w'*M(materialProp)*f""" - if(materialProp is None): - materialProp = np.ones(self.nC) - Av = self.faceAve - return sdiag(Av.T * (self.vol * mkvc(materialProp))) + return self.getMass(loc='f', materialProp=materialProp, inv=inv) def getFaceMassDeriv(self): Av = self.faceAve From 0bebcd48add4b73bb1acf7a01114187bb297d24b Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 18 Oct 2013 14:32:33 -0700 Subject: [PATCH 3/6] Changed kwarg order for easier use. --- SimPEG/DiffOperators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 86fd35f7..b2a7a6a5 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -286,7 +286,7 @@ class DiffOperators(object): _nodalVectorAve = None nodalVectorAve = property(**nodalVectorAve()) - def getMass(self, loc='e', materialProp=None, inv=False): + def getMass(self, materialProp=None, loc='e', inv=False): """ Produces mass matricies. Kwargs: From 1c3a04f3375386e1645cbc722054f7e084f22828 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Fri, 18 Oct 2013 14:41:26 -0700 Subject: [PATCH 4/6] Changed nomenclature for averaging operators (issue #6) --- SimPEG/DiffOperators.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index b2a7a6a5..88756ad8 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -203,46 +203,46 @@ class DiffOperators(object): _edgeCurl = None edgeCurl = property(**edgeCurl()) - def faceAve(): + def aveF2CC(): doc = "Construct the averaging operator on cell faces to cell centers." def fget(self): - if(self._faceAve is None): + if(self._aveF2CC is None): n = self.n if(self.dim == 1): - self._faceAve = av(n[0]) + self._aveF2CC = av(n[0]) elif(self.dim == 2): - self._faceAve = sp.hstack((sp.kron(speye(n[1]), av(n[0])), + self._aveF2CC = 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._faceAve = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + self._aveF2CC = 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._faceAve + return self._aveF2CC return locals() - _faceAve = None - faceAve = property(**faceAve()) + _aveF2CC = None + aveF2CC = property(**aveF2CC()) - def edgeAve(): + def aveE2CC(): doc = "Construct the averaging operator on cell edges to cell centers." def fget(self): - if(self._edgeAve is None): + if(self._aveE2CC is None): # The number of cell centers in each direction n = self.n if(self.dim == 1): raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') elif(self.dim == 2): - self._edgeAve = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + self._aveE2CC = 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._edgeAve = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + self._aveE2CC = 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._edgeAve + return self._aveE2CC return locals() - _edgeAve = None - edgeAve = property(**edgeAve()) + _aveE2CC = None + aveE2CC = property(**aveE2CC()) def nodalAve(): doc = "Construct the averaging operator on cell nodes to cell centers." @@ -316,9 +316,9 @@ class DiffOperators(object): assert materialProp.shape == (self.nC,), "materialProp incorrect shape" if loc=='e': - Av = self.edgeAve + Av = self.aveE2CC elif loc=='f': - Av = self.faceAve + Av = self.aveF2CC else: raise ValueError('Invalid loc') @@ -338,5 +338,5 @@ class DiffOperators(object): return self.getMass(loc='f', materialProp=materialProp, inv=inv) def getFaceMassDeriv(self): - Av = self.faceAve + Av = self.aveF2CC return Av.T * sdiag(self.vol) From 04dfc83e1a74a2d1c51e21070944fe15476783fc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 18 Oct 2013 15:27:02 -0700 Subject: [PATCH 5/6] documentation and a few more renaming. --- SimPEG/DiffOperators.py | 58 +++++++++++++++---------------- SimPEG/LogicallyOrthogonalMesh.py | 2 +- 2 files changed, 29 insertions(+), 31 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 88756ad8..f26252d9 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -244,61 +244,59 @@ class DiffOperators(object): _aveE2CC = None aveE2CC = property(**aveE2CC()) - def nodalAve(): + def aveN2CC(): doc = "Construct the averaging operator on cell nodes to cell centers." def fget(self): - if(self._nodalAve is None): + if(self._aveN2CC is None): # The number of cell centers in each direction n = self.n if(self.dim == 1): - self._nodalAve = av(n[0]) + self._aveN2CC = av(n[0]) elif(self.dim == 2): - self._nodalAve = sp.hstack((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") + self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])), + sp.kron(av(n[1]), av(n[0]))), format="csr") elif(self.dim == 3): - self._nodalAve = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._nodalAve + self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") + return self._aveN2CC return locals() - _nodalAve = None - nodalAve = property(**nodalAve()) + _aveN2CC = None + aveN2CC = property(**aveN2CC()) - def nodalVectorAve(): + def aveN2CCv(): doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." def fget(self): - if(self._nodalVectorAve is None): + if(self._aveN2CCv is None): # The number of cell centers in each direction n = self.n if(self.dim == 1): - self._nodalVectorAve = av(n[0]) + self._aveN2CCv = av(n[0]) elif(self.dim == 2): - self._nodalVectorAve = sp.block_diag((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") + self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])), + sp.kron(av(n[1]), av(n[0]))), format="csr") elif(self.dim == 3): - self._nodalVectorAve = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._nodalVectorAve + self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0])), + kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") + return self._aveN2CCv return locals() - _nodalVectorAve = None - nodalVectorAve = property(**nodalVectorAve()) + _aveN2CCv = None + aveN2CCv = property(**aveN2CCv()) def getMass(self, materialProp=None, loc='e', inv=False): """ Produces mass matricies. - Kwargs: - loc (str): 'e' - Average to edges - 'f' faces - materialProp: property to be averaged (see below) - inv (bool): True returns matrix inverse + :param str loc: Average to location: 'e'-edges, 'f'-faces + :param None,float,numpy.ndarray materialProp: property to be averaged (see below) + :param bool inv: True returns matrix inverse + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the mass matrix - Returns: - scipy.sparse.csr.csr_matrix + materialProp can be:: - materialProp can be: None -> takes materialProp = 1 (default) float -> a constant value for entire domain numpy.ndarray -> if materialProp.size == self.nC diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index 07c37255..5d2b7ffa 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -45,7 +45,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): def fget(self): if self._gridCC is None: - ccV = (self.nodalVectorAve*mkvc(self.gridN)) + ccV = (self.aveN2CCv*mkvc(self.gridN)) self._gridCC = ccV.reshape((-1, self.dim), order='F') return self._gridCC return locals() From e56e5cbf05986bd5b8507bde3a5f3bce26511b73 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 18 Oct 2013 17:04:23 -0700 Subject: [PATCH 6/6] Removed inversion option. --- SimPEG/DiffOperators.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index f26252d9..110fe9bd 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -286,12 +286,11 @@ class DiffOperators(object): _aveN2CCv = None aveN2CCv = property(**aveN2CCv()) - def getMass(self, materialProp=None, loc='e', inv=False): + def getMass(self, materialProp=None, loc='e'): """ Produces mass matricies. :param str loc: Average to location: 'e'-edges, 'f'-faces :param None,float,numpy.ndarray materialProp: property to be averaged (see below) - :param bool inv: True returns matrix inverse :rtype: scipy.sparse.csr.csr_matrix :return: M, the mass matrix @@ -322,18 +321,15 @@ class DiffOperators(object): diag = Av.T * (self.vol * mkvc(materialProp)) - if inv: - diag = 1/diag - return sdiag(diag) - def getEdgeMass(self, materialProp=None, inv=False): + def getEdgeMass(self, materialProp=None): """mass matrix for products of edge functions w'*M(materialProp)*e""" - return self.getMass(loc='e', materialProp=materialProp, inv=inv) + return self.getMass(loc='e', materialProp=materialProp) - def getFaceMass(self, materialProp=None, inv=False): + def getFaceMass(self, materialProp=None): """mass matrix for products of face functions w'*M(materialProp)*f""" - return self.getMass(loc='f', materialProp=materialProp, inv=inv) + return self.getMass(loc='f', materialProp=materialProp) def getFaceMassDeriv(self): Av = self.aveF2CC