diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index a6fb21f7..43d6c31f 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -101,7 +101,7 @@ class InnerProducts(object): If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: - P = [P000, P001, P010, P011, P100, P101, P110, P111] + P = [P000, P100, P010, P110, P001, P101, P011, P111] Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: @@ -205,7 +205,7 @@ class InnerProducts(object): If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: - P = [P000, P001, P010, P011, P100, P101, P110, P111] + P = [P000, P100, P010, P110, P001, P101, P011, P111] Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index dec87423..ebd089ef 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1,4 +1,6 @@ from SimPEG import np, sp, Utils, Solver +from BaseMesh import BaseMesh +from InnerProducts import InnerProducts import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import matplotlib.colors as colors @@ -582,15 +584,16 @@ class TreeCell(TreeObject): self.mesh.cells.remove(self) - @property - def faceIndex(self): + def faceIndex(self, theseFaces='all', addDirection=True): #TODO: preallocate I, J, V = np.empty(0,dtype=float), np.empty(0,dtype=float), np.empty(0,dtype=float) for face in self.faces: + thisFace = 'all' == theseFaces or face in theseFaces + if not thisFace: continue j = self.faces[face].index i = j*0+self.num v = j*0+1 - if 'm' in face: + if addDirection and 'm' in face: v *= -1 I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] return I, J, V @@ -598,7 +601,6 @@ class TreeCell(TreeObject): @property def vol(self): return self.sz.prod() - def viz(self, ax, color='none', text=False): if not self.isleaf: return x0, sz = self.x0, self.sz @@ -615,8 +617,11 @@ class TreeCell(TreeObject): if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) -class TreeMesh(object): +class TreeMesh(InnerProducts, BaseMesh): """TreeMesh""" + + _meshType = 'TREE' + def __init__(self, h_in, x0=None): assert type(h_in) is list, 'h_in must be a list' h = range(len(h_in)) @@ -635,7 +640,9 @@ class TreeMesh(object): assert type(x0) in [list, np.ndarray], 'x0 must be a numpy array or a list' x0 = np.array(x0, dtype=float) assert len(x0) == self.dim, 'x0 must have the same dimensions as the mesh' - self.x0 = x0 + + # TODO: this has a lot of stuff which doesn't work for this style of mesh... + BaseMesh.__init__(self, np.array([x.size for x in h]), x0) # set the sets for holding the cells, nodes, faces, and edges self.cells = set() @@ -846,7 +853,7 @@ class TreeMesh(object): # TODO: Preallocate! I, J, V = np.empty(0), np.empty(0), np.empty(0) for cell in self.sortedCells: - i, j, v = cell.faceIndex + i, j, v = cell.faceIndex('all') I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] VOL = self.vol @@ -855,6 +862,43 @@ class TreeMesh(object): self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S) return self._faceDiv + def _getFacePxx(self, xFace, yFace): + self.number() + I, J, V = np.empty(0), np.empty(0), np.empty(0) + for cell in self.sortedCells: + i, j, v = cell.faceIndex('fX'+xFace) + I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] + xP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + + I, J, V = np.empty(0), np.empty(0), np.empty(0) + for cell in self.sortedCells: + i, j, v = cell.faceIndex('fY'+yFace) + I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] + yP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + return sp.vstack((xP, yP)) + + def _getFacePxxx(self, xFace, yFace, zFace): + self.number() + I, J, V = np.empty(0), np.empty(0), np.empty(0) + for cell in self.sortedCells: + i, j, v = cell.faceIndex('fX'+xFace) + I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] + xP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + + I, J, V = np.empty(0), np.empty(0), np.empty(0) + for cell in self.sortedCells: + i, j, v = cell.faceIndex('fY'+yFace) + I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] + yP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + + I, J, V = np.empty(0), np.empty(0), np.empty(0) + for cell in self.sortedCells: + i, j, v = cell.faceIndex('fZ'+zFace) + I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v] + zP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + + return sp.vstack((xP, yP, zP)) + def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, plotE=False, plotEx=False, plotEy=False, plotEz=False, showIt=False): axOpts = {'projection':'3d'} if self.dim == 3 else {} if ax is None: ax = plt.subplot(111, **axOpts) diff --git a/SimPEG/Tests/test_TreeMesh.py b/SimPEG/Tests/test_TreeMesh.py index 4abdfd05..10d1dba6 100644 --- a/SimPEG/Tests/test_TreeMesh.py +++ b/SimPEG/Tests/test_TreeMesh.py @@ -489,13 +489,15 @@ class SimpleOctreeOperatorTests(unittest.TestCase): h1 = np.random.rand(5) h2 = np.random.rand(7) h3 = np.random.rand(3) - self.tM = TensorMesh([h1,h2,1]) - self.oM = TreeMesh([h1,h2,1]) + self.tM = TensorMesh([h1,h2,h3]) + self.oM = TreeMesh([h1,h2,h3]) def test_faceDiv(self): - print (self.tM.faceDiv - self.oM.faceDiv) self.assertTrue((self.tM.faceDiv - self.oM.faceDiv).toarray().sum() == 0) + def test_InnerProducts(self): + self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() == 0) + if __name__ == '__main__': unittest.main()