diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index c44afebe..f797c1c5 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -3,39 +3,10 @@ import scipy.sparse as sp from scipy.constants import pi from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap from TensorMesh import BaseTensorMesh +from InnerProducts import InnerProducts - # def getNearest(self, loc, locType): - # """ Returns the index of the closest face or edge to a given location - - # :param numpy.ndarray loc: Test point - # :param str locType: Type of location desired (see below) - # :rtype: int - # :return: ind: - - # locType can be:: - - # 'fz' -> location of nearest z-face - # 'fr' -> location of nearest r-face - # 'et' -> location of nearest edge - # """ - - # if locType=='et': - # dr = self.gridN[:,0] - loc[0] - # dz = self.gridN[:,1] - loc[1] - # elif locType=='fz': - # dr = self.gridFz[:,0] - loc[0] - # dz = self.gridFz[:,1] - loc[1] - # elif locType=='fr': - # dr = self.gridFr[:,0] - loc[0] - # dz = self.gridFr[:,1] - loc[1] - # else: - # raise ValueError('Invalid locType') - # R = np.sqrt(dr**2 + dz**2) - # ind = np.argmin(R) - # return ind - -class CylMesh(BaseTensorMesh): +class CylMesh(BaseTensorMesh, InnerProducts): """ CylMesh is a mesh class for cylindrical problems """ @@ -49,6 +20,10 @@ class CylMesh(BaseTensorMesh): assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" + @property + def isSymmetric(self): + return self.nCy == 1 + @property def nNx(self): """ @@ -57,7 +32,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNx """ - if self.nCy == 1: return self.nCx + if self.isSymmetric: return self.nCx return self.nCx + 1 @property @@ -68,7 +43,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNy """ - if self.nCy == 1: return 0 + if self.isSymmetric: return 0 return self.nCy @property @@ -89,7 +64,7 @@ class CylMesh(BaseTensorMesh): :rtype: numpy.array (dim, ) :return: vnEy or None if dim < 2 """ - nNx = self.nNx if self.nCy == 1 else self.nNx - 1 + nNx = self.nNx if self.isSymmetric else self.nNx - 1 return np.r_[nNx, self.nCy, self.nNz] @property @@ -100,7 +75,7 @@ class CylMesh(BaseTensorMesh): :rtype: numpy.array (dim, ) :return: vnEz or None if nCy > 1 """ - if self.nCy == 1: + if self.isSymmetric: return np.r_[self.nNx, self.nNy, self.nCz] else: return None @@ -113,7 +88,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nEz """ - if self.nCy == 1: + if self.isSymmetric: return self.vnEz.prod() return (np.r_[self.nNx-1, self.nNy, self.nCz]).prod() + self.nCz @@ -130,14 +105,14 @@ class CylMesh(BaseTensorMesh): @property def vectorNx(self): """Nodal grid vector (1D) in the x direction.""" - if self.nCy == 1: + if self.isSymmetric: return self.hx.cumsum() return np.r_[0, self.hx].cumsum() @property def vectorNy(self): """Nodal grid vector (1D) in the y direction.""" - if self.nCy == 1: + if self.isSymmetric: # There aren't really any nodes, but all the grids need # somewhere to live, why not zero?! return np.r_[0] @@ -147,7 +122,7 @@ class CylMesh(BaseTensorMesh): def edge(self): """Edge lengths""" if getattr(self, '_edge', None) is None: - if self.nCy == 1: + if self.isSymmetric: self._edge = 2*pi*self.gridN[:,0] else: raise NotImplementedError('edges not yet implemented for 3D cyl mesh') @@ -186,7 +161,7 @@ class CylMesh(BaseTensorMesh): # Compute faceDivergence operator on faces D1 = self.faceDivx D3 = self.faceDivz - if self.nCy == 1: + if self.isSymmetric: D = sp.hstack((D1, D3), format="csr") elif self.nCy > 1: D2 = self.faceDivy @@ -237,7 +212,7 @@ class CylMesh(BaseTensorMesh): def nodalGrad(self): """Construct gradient operator (nodes to edges).""" # Nodal grad does not make sense for cylindrically symmetric mesh. - if self.nCy == 1: return None + if self.isSymmetric: return None raise NotImplementedError('nodalGrad not yet implemented') @property @@ -263,178 +238,49 @@ class CylMesh(BaseTensorMesh): self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0) return self._edgeCurl + # @property + # def aveE2CC(self): + # """Averaging operator from cell edges to cell centers""" + # if getattr(self, '_aveE2CC', None) is None: + # if self.isSymmetric: + # az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') + # ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') + # ar[0,0] = 1 + # self._aveE2CC = (0.5)*sp.kron(az, ar).T + # else: + # raise NotImplementedError('wrapping in the averaging is not yet implemented') + # return self._aveE2CC + @property def aveE2CC(self): - """Averaging operator from cell edges to cell centres""" + "Construct the averaging operator on cell edges to cell centers." if getattr(self, '_aveE2CC', None) is None: - az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') - ar[0,0] = 1 - self._aveE2CC = sp.kron(az, ar).T + # The number of cell centers in each direction + n = self.vnC + if self.isSymmetric: + avR = av(n[0])[:,1:] + avR[0,0] = 1. + self._aveE2CC = (0.5)*sp.kron(av(n[2]), avR, format="csr") + else: + raise NotImplementedError('wrapping in the averaging is not yet implemented') + # 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 + @property def aveF2CC(self): - """Averaging operator from cell faces to cell centres""" + "Construct the averaging operator on cell faces to cell centers." if getattr(self, '_aveF2CC', None) is None: - az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') - ar[0,0] = 1 - Afr = sp.kron(sp.eye(self.nCz),ar) - Afz = sp.kron(az,sp.eye(self.nCx)) - self._aveF2CC = sp.vstack((Afr,Afz)).T + n = self.vnC + if self.isSymmetric: + avR = av(n[0])[:,1:] + avR[0,0] = 1. + self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[2]), avR), sp.kron(av(n[2]), speye(n[0]))), format="csr") + else: + raise NotImplementedError('wrapping in the averaging is not yet implemented') + # 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 - - def getFaceMassDeriv(self): - Av = self.aveF2CC - return Av.T * sdiag(self.vol) - - def getEdgeMassDeriv(self): - Av = self.aveE2CC - return Av.T * sdiag(self.vol) - - - #################################################### - # Methods - #################################################### - - - def getMass(self, materialProp=None, loc='e'): - """ Produces mass matricies. - - :param None,float,numpy.ndarray materialProp: property to be averaged (see below) - :param str loc: Average to location: 'e'-edges, 'f'-faces - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the mass 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) - materialProp = mkvc(materialProp) - assert materialProp.shape == (self.nC,), "materialProp incorrect shape" - - if loc=='e': - Av = self.aveE2CC - elif loc=='f': - Av = self.aveF2CC - else: - raise ValueError('Invalid loc') - - diag = Av.T * (self.vol * mkvc(materialProp)) - - return sdiag(diag) - - def getEdgeMass(self, materialProp=None): - """mass matrix for products of edge functions w'*M(materialProp)*e""" - return self.getMass(loc='e', materialProp=materialProp) - - def getFaceMass(self, materialProp=None): - """mass matrix for products of face functions w'*M(materialProp)*f""" - return self.getMass(loc='f', materialProp=materialProp) - - # def getInterpolationMat(self, loc, locType='fz'): - # """ Produces intrpolation matrix - - # :param numpy.ndarray loc: Location of points to interpolate to - # :param str locType: What to interpolate (see below) - # :rtype: scipy.sparse.csr.csr_matrix - # :return: M, the intrpolation matrix - - # locType can be:: - - # 'fz' -> z-component of field defined on faces - # 'fr' -> r-component of field defined on faces - # 'et' -> theta-component of field defined on edges - # """ - - # loc = np.atleast_2d(loc) - - # assert np.all(loc[:,0]<=self.vectorNx.max()) & \ - # np.all(loc[:,1]>=self.vectorNz.min()) & \ - # np.all(loc[:,1]<=self.vectorNz.max()), \ - # "Points outside of mesh" - - - # if locType=='fz': - # Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float) - - # for i, iloc in enumerate(loc): - # # Point is on a z-interface - # if np.any(np.abs(self.vectorNz-iloc[1])<0.001): - # dFz = self.gridFz-iloc #Distance to z faces - # dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left... - # indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one - # if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - # zFL = self.gridFz[indL,:] - # zFLL = self.gridFz[indL-1,:] - # Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0]) - # Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0]) - # else: - # zFL = self.gridFz[indL,:] - # zFR = self.gridFz[indL+1,:] - # Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0]) - # Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0]) - # # Point is in a cell - # else: - # dFz = self.gridFz-iloc - # dFz[dFz>0] = np.inf - # dFz = np.sum(dFz**2, axis=1) - - # indBL = np.argmin(dFz) # Face below and to the left - # indAL = indBL + self.nCx # Face above and to the left - - # zF_BL = self.gridFz[indBL,:] - # zF_AL = self.gridFz[indAL,:] - - # dzB = iloc[1] - zF_BL[1] # z-distance to face below - # dzA = zF_AL[1] - iloc[1] # z-distance to face above - - # if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - # zF_BLL = self.gridFz[indBL-1,:] - # zF_ALL = self.gridFz[indAL-1,:] - - # DZ = zF_AL[1] - zF_BL[1] - # DR = zF_AL[0] - zF_ALL[0] - - # drL = iloc[0] - zF_AL[0] - # drLL = iloc[0] - zF_ALL[0] - - # Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR) - # Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR) - # Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR) - # Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR) - # else: - # indBR = indBL+1 # Face below and to the right - # indAR = indAL + 1 # Face above and to the right - # zF_BR = self.gridFz[indBR,:] - - # drL = iloc[0] - zF_BL[0] # r-distance to face on left - # drR = zF_BR[0] - iloc[0] # r-distance to face on right - - # drz = (drL + drR)*(dzB + dzA) - # Q[i,indBL+self.nFr] = drR*dzA/drz - # Q[i,indBR+self.nFr] = drL*dzA/drz - # Q[i,indAL+self.nFr] = drR*dzB/drz - # Q[i,indAR+self.nFr] = drL*dzB/drz - - # elif locType=='fr': - # raise NotImplementedError('locType==fr') - # elif locType=='et': - # raise NotImplementedError('locType==et') - # else: - # raise ValueError('Invalid locType') - # return Q.tocsr() - - diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 15d8ed2b..f8d448d4 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -623,11 +623,11 @@ class DiffOperators(object): raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') 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") + 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") + 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 @property @@ -696,57 +696,3 @@ class DiffOperators(object): kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") return self._aveN2F - # --------------- Methods --------------------- - - 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) - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the mass 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.aveE2CC - elif loc=='f': - Av = self.aveF2CC - else: - raise ValueError('Invalid loc') - - diag = Av.T * (self.vol * mkvc(materialProp)) - - return sdiag(diag) - - def getEdgeMass(self, materialProp=None): - """mass matrix for products of edge functions w'*M(materialProp)*e""" - return self.getMass(loc='e', materialProp=materialProp) - - def getFaceMass(self, materialProp=None): - """mass matrix for products of face functions w'*M(materialProp)*f""" - return self.getMass(loc='f', materialProp=materialProp) - - def getFaceMassDeriv(self): - Av = self.aveF2CC - return Av.T * sdiag(self.vol) - - def getEdgeMassDeriv(self): - Av = self.aveE2CC - return Av.T * sdiag(self.vol) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index cc3bcbb1..4d1b16a2 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -239,6 +239,8 @@ class BaseTensorMesh(BaseRectangularMesh): nF_nE = self.vnF if 'F' in locType else self.vnE components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) + # remove any zero blocks (hstack complains) + components = [comp for comp in components if comp.shape[1] > 0] Q = sp.hstack(components) elif locType in ['CC', 'N']: Q = Utils.interpmat(loc, *self.getTensor(locType)) @@ -251,6 +253,116 @@ class BaseTensorMesh(BaseRectangularMesh): return Q.tocsr() + def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getEdgeInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str AvType: 'E' or 'F' + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None: + materialProperty = np.ones(self.nC) + + if invertProperty: + materialProperty = 1./materialProperty + + if Utils.isScalar(materialProperty): + materialProperty = materialProperty*np.ones(self.nC) + + if materialProperty.size == self.nC: + Av = getattr(self, 'ave'+AvType+'2CC') + Vprop = self.vol * Utils.mkvc(materialProperty) + return self.dim * Utils.sdiag(Av.T * Vprop) + if materialProperty.size == self.nC*self.dim: + Av = getattr(self, 'ave'+AvType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) + + + def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v) + + + def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v) + + + def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None): + """ + :param str AvType: 'E' or 'F' + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None: + return 0.0 + + if Utils.isScalar(materialProperty): + Av = getattr(self, 'ave'+AvType+'2CC') + V = Utils.sdiag(self.vol) + ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) + if v is None: + return self.dim * Av.T * V * ones + return Utils.sdiag(v) * self.dim * Av.T * V * ones + + if materialProperty.size == self.nC: + Av = getattr(self, 'ave'+AvType+'2CC') + V = Utils.sdiag(self.vol) + if v is None: + return self.dim * Av.T * V + return Utils.sdiag(v) * self.dim * Av.T * V + + if materialProperty.size == self.nC*self.dim: # anisotropic + Av = getattr(self, 'ave'+AvType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + if v is None: + return Av.T * V + return Utils.sdiag(v) * Av.T * V + + class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): """ @@ -449,111 +561,7 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) return indxd, indxu, indyd, indyu, indzd, indzu - def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor materialProperty. - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :param bool invertProperty: inverts the material property - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty) - - - def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False): - """ - Fast version of getEdgeInnerProduct. - This does not handle the case of a full tensor materialProperty. - - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :param bool invertProperty: inverts the material property - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nE, nE) - """ - return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty) - - - def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor materialProperty. - - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param str AvType: 'E' or 'F' - :param bool returnP: returns the projection matrices - :param bool invertProperty: inverts the material property - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - if materialProperty is None: - materialProperty = np.ones(self.nC) - - if invertProperty: - materialProperty = 1./materialProperty - - if Utils.isScalar(materialProperty): - materialProperty = materialProperty*np.ones(self.nC) - - if materialProperty.size == self.nC: - Av = getattr(self, 'ave'+AvType+'2CC') - Vprop = self.vol * Utils.mkvc(materialProperty) - return self.dim * Utils.sdiag(Av.T * Vprop) - if materialProperty.size == self.nC*self.dim: - Av = getattr(self, 'ave'+AvType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) - - - def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None): - """ - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v) - - - def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None): - """ - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nE, nE) - """ - return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v) - - - def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None): - """ - :param str AvType: 'E' or 'F' - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - if materialProperty is None: - return None - if Utils.isScalar(materialProperty): - Av = getattr(self, 'ave'+AvType+'2CC') - V = Utils.sdiag(self.vol) - ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) - if v is None: - return self.dim * Av.T * V * ones - return Utils.sdiag(v) * self.dim * Av.T * V * ones - if materialProperty.size == self.nC: - Av = getattr(self, 'ave'+AvType+'2CC') - V = Utils.sdiag(self.vol) - if v is None: - return self.dim * Av.T * V - return Utils.sdiag(v) * self.dim * Av.T * V - if materialProperty.size == self.nC*self.dim: # anisotropic - Av = getattr(self, 'ave'+AvType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - if v is None: - return Av.T * V - return Utils.sdiag(v) * Av.T * V if __name__ == '__main__': diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index cc6d6338..85a87bdc 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -337,7 +337,7 @@ if __name__ == '__main__': D = M.faceDiv G = M.cellGrad - Msig = M.getFaceMass() + Msig = M.getFaceInnerProduct() A = D*Msig*G A[0,0] *= 10 # remove the constant null space from the matrix diff --git a/SimPEG/Tests/test_Solver.py b/SimPEG/Tests/test_Solver.py index f8a8983f..76177990 100644 --- a/SimPEG/Tests/test_Solver.py +++ b/SimPEG/Tests/test_Solver.py @@ -22,7 +22,7 @@ class TestSolver(unittest.TestCase): D = M.faceDiv G = M.cellGrad - Msig = M.getFaceMass() + Msig = M.getFaceInnerProduct() A = D*Msig*G A[0,0] *= 10 # remove the constant null space from the matrix diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index f1c6e423..5eaef428 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -206,6 +206,62 @@ class TestEdgeCurl2D(OrderTest): def test_order(self): self.orderTest() + +# class TestInnerProducts2D(OrderTest): +# """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" + +# meshTypes = MESHTYPES +# meshDimension = 2 +# meshSizes = [4, 8, 16, 32, 64, 128] + +# def getError(self): + +# funR = lambda r, t, z: np.cos(2.0*np.pi*z) +# funT = lambda r, t, z: 0*t +# funZ = lambda r, t, z: np.sin(2.0*np.pi*r) + +# call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) + +# sigma1 = lambda r, t, z: z+1 +# sigma2 = lambda r, t, z: r*z+50 +# sigma3 = lambda r, t, z: 3+t*r +# sigma4 = lambda r, t, z: 0.1*r*t*z +# sigma5 = lambda r, t, z: 0.2*z*r*t +# sigma6 = lambda r, t, z: 0.1*t + +# Gc = self.M.gridCC +# if self.sigmaTest == 1: +# sigma = np.c_[call(sigma1, Gc)] +# analytic = 144877./360 # Found using sympy. z=5 +# elif self.sigmaTest == 2: +# sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)] +# analytic = 189959./120 # Found using sympy. z=5 +# elif self.sigmaTest == 3: +# sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] +# analytic = 781427./360 # Found using sympy. z=5 + +# if self.location == 'edges': +# E = call(funT, self.M.gridEy) +# A = self.M.getEdgeInnerProduct(sigma) +# numeric = E.T.dot(A.dot(E)) +# elif self.location == 'faces': +# Fr = call(funR, self.M.gridFx) +# Fz = call(funZ, self.M.gridFz) +# A = self.M.getFaceInnerProduct(sigma) +# F = np.r_[Fr,Fz] +# numeric = F.T.dot(A.dot(F)) + +# print numeric +# err = np.abs(numeric - analytic) +# return err + +# def test_order1_faces(self): +# self.name = "2D Face Inner Product - Isotropic" +# self.location = 'faces' +# self.sigmaTest = 1 +# self.orderTest() + + class TestCyl3DMesh(unittest.TestCase): def setUp(self):