diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 9a133e0d..40512b2b 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -502,12 +502,17 @@ class Mesh2MeshTopo(IdentityMap): if self.tree==None: self.tree = cKDTree(zip(self.mesh.gridCC[self.actind,0], self.mesh.gridCC[self.actind,1], self.mesh.gridCC[self.actind,2])) d, inds = self.tree.query(zip(self.mesh2.gridCC[self.actind2,0],self.mesh2.gridCC[self.actind2,1],self.mesh2.gridCC[self.actind2,2]), k=self.nIterpPts) - w = 1./ d**2 - w = Utils.sdiag(1./np.sum(w, axis=1)) * w - I = Utils.mkvc(np.arange(inds.shape[0]).reshape([-1,1]).repeat(6, axis=1)) + # Not sure consideration of the volume ... + # vol = np.zeros((self.actind2.sum(), self.nIterpPts)) + # for i in range(self.nIterpPts): + # vol[:,i] = self.mesh.vol[inds[:,i]] + w = 1. / d**2 + w = Utils.sdiag(1./np.sum(w, axis=1)) * (w) + I = Utils.mkvc(np.arange(inds.shape[0]).reshape([-1,1]).repeat(self.nIterpPts, axis=1)) J = Utils.mkvc(inds) P = sp.coo_matrix( (Utils.mkvc(w),(I, J)), shape=(inds.shape[0], (self.actind).sum()) ) - self.P = P.tocsc() + # self.P = Utils.sdiag(self.mesh2.vol[self.actind2])*P.tocsc() + self.P = P.tocsc() @property def shape(self): diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 7a738cf2..201ff56c 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -320,10 +320,12 @@ class GlobalProblem(BaseProblem): subMesh = self.meshes[ind] # subMap = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping? + # subMap = self.getSubMap(subMesh, ind) if self.PropMap is None: prob = self.SubProblem(subMesh, mapping=subMap * self.mapping, **self.probKwargs) else: + # This will not work with a fancier propmap... prob = self.SubProblem(subMesh, mapping=subMap * self._propMapMapping, **self.probKwargs) survey = self.survey.__class__(srcList=self.survey.srcList[self.groups[ind]], **self.surveyKwargs) @@ -331,6 +333,18 @@ class GlobalProblem(BaseProblem): return prob, survey + # Not sure we need this here ... + + # def getSubMap(self, subMesh, ind): + # """The sub""" + # mesh2mesh = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping? + + # if self.PropMap is None: + # subMap = mesh2mesh * self.mapping + # else: + # subMap = mesh2mesh * self._propMapMapping + + # return subMap if __name__ == '__main__':