From 5afea4e2b60119bf59ad6a5c49afe00baa39e626 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 11 Feb 2015 16:36:57 -0800 Subject: [PATCH] update grids. Think refining is working properly --- SimPEG/Mesh/NewTreeMesh.py | 155 ++++++++++++++++++++++++++----- SimPEG/Tests/test_NewTreeMesh.py | 114 ++++++++++++++++++++--- 2 files changed, 229 insertions(+), 40 deletions(-) diff --git a/SimPEG/Mesh/NewTreeMesh.py b/SimPEG/Mesh/NewTreeMesh.py index 96e525c0..462f0c74 100644 --- a/SimPEG/Mesh/NewTreeMesh.py +++ b/SimPEG/Mesh/NewTreeMesh.py @@ -232,6 +232,10 @@ class TreeCell(TreeIndexer): def n6(self): return self.fZp.n2 @property def n7(self): return self.fZp.n3 + @property + def center(self): + return (self.n0.vec + self.n1.vec + self.n2.vec + self.n3.vec + + self.n4.vec + self.n5.vec + self.n6.vec + self.n7.vec)/8.0 class TreeMesh(BaseMesh): @@ -533,8 +537,11 @@ class TreeMesh(BaseMesh): @property def gridCC(self): - F = TreeFace(self, 'active') - return F.sort(F.center) + if self.dim == 2: + F = TreeFace(self, 'active') + return F.sort(F.center) + C = TreeCell(self, 'active') + return C.sort(C.center) @property def gridEx(self): @@ -660,23 +667,23 @@ class TreeMesh(BaseMesh): nEi, nE = self._push('_edges', nE) # Add four new faces - Fs = np.zeros((4,8)) - Fs[:, ACTIVE] = 1 - Fs[:, PARENT] = index - Fs[:, FDIR] = f[FDIR] + nF = np.zeros((4,8)) + nF[:, ACTIVE] = 1 + nF[:, PARENT] = index + nF[:, FDIR] = f[FDIR] - fInds = [FEDGE0,FEDGE1,FEDGE2,FEDGE3] - Fs[0, fInds] = [E0i[0], nEi[0], E2i[0], nEi[2]] - Fs[1, fInds] = [E0i[1], nEi[1], nEi[2], E3i[0]] - Fs[2, fInds] = [nEi[0], E1i[0], E2i[1], nEi[3]] - Fs[3, fInds] = [nEi[1], E1i[1], nEi[3], E3i[1]] + feInds = [FEDGE0,FEDGE1,FEDGE2,FEDGE3] + nF[0, feInds] = [E0i[0], nEi[0], E2i[0], nEi[2]] + nF[1, feInds] = [E0i[1], nEi[1], nEi[2], E3i[0]] + nF[2, feInds] = [nEi[0], E1i[0], E2i[1], nEi[3]] + nF[3, feInds] = [nEi[1], E1i[1], nEi[3], E3i[1]] - return self._push('_faces', Fs) + return self._push('_faces', nF) def refineCell(self, index): c = self._cells[index,:] - if f[ACTIVE] == 0: + if c[ACTIVE] == 0: # search for the children up to one level deep subInds = np.argwhere(self._cells[:,PARENT] == index).flatten() return subInds, self._cells[subInds,:] @@ -684,16 +691,114 @@ class TreeMesh(BaseMesh): self._cells[index, ACTIVE] = 0 # Refine the outer faces - F0i, F0 = self.refineFace(c[CFACE0]) - F1i, F1 = self.refineFace(c[CFACE1]) - F2i, F2 = self.refineFace(c[CFACE2]) - F3i, F3 = self.refineFace(c[CFACE3]) - F4i, F4 = self.refineFace(c[CFACE4]) - F5i, F5 = self.refineFace(c[CFACE5]) + fXm, rfXm = self.refineFace(c[CFACE0]) + fXp, rfXp = self.refineFace(c[CFACE1]) + fYm, rfYm = self.refineFace(c[CFACE2]) + fYp, rfYp = self.refineFace(c[CFACE3]) + fZm, rfZm = self.refineFace(c[CFACE4]) + fZp, rfZp = self.refineFace(c[CFACE5]) - nodeNums = self._edges[f[[FEDGE0, FEDGE1]],:][:,[ENODE0, ENODE1]] + nodes = TreeCell(self, index).fXm.nodes + TreeCell(self, index).fXp.nodes + nodeNums = [n.index for n in nodes] newNode, node = self.addNode(nodeNums) + nE = np.zeros((6,6)) + nE[:, ACTIVE] = 1 + nE[:, PARENT] = -1 + nE[:, EDIR] = [0, 0, 1, 1, 2, 2] + + nE[0, ENODE0] = TreeFace(self, fXm[0]).n3.index + nE[0, ENODE1] = newNode + nE[1, ENODE0] = newNode + nE[1, ENODE1] = TreeFace(self, fXp[0]).n3.index + + nE[2, ENODE0] = TreeFace(self, fYm[0]).n3.index + nE[2, ENODE1] = newNode + nE[3, ENODE0] = newNode + nE[3, ENODE1] = TreeFace(self, fYp[0]).n3.index + + nE[4, ENODE0] = TreeFace(self, fZm[0]).n3.index + nE[4, ENODE1] = newNode + nE[5, ENODE0] = newNode + nE[5, ENODE1] = TreeFace(self, fZp[0]).n3.index + + nEi, nE = self._push('_edges', nE) + + nF = np.zeros((12,8)) + nF[:, ACTIVE] = 1 + nF[:, PARENT] = -1 + nF[:, FDIR] = [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2] + + feInds = [FEDGE0,FEDGE1,FEDGE2,FEDGE3] + nF[0, feInds] = [rfZm[0, FEDGE3], nEi[2], rfYm[0, FEDGE3], nEi[4]] + nF[1, feInds] = [rfZm[2, FEDGE3], nEi[3], nEi[4], rfYp[0, FEDGE3]] + nF[2, feInds] = [nEi[2], rfZp[0, FEDGE3], rfYm[2, FEDGE3], nEi[5]] + nF[3, feInds] = [nEi[3], rfZp[2, FEDGE3], nEi[5], rfYp[2, FEDGE3]] + + nF[4, feInds] = [rfZm[0, FEDGE1], nEi[0], rfXm[0, FEDGE3], nEi[4]] + nF[5, feInds] = [rfZm[1, FEDGE1], nEi[1], nEi[4], rfXp[0, FEDGE3]] + nF[6, feInds] = [nEi[0], rfZp[0, FEDGE1], rfXm[2, FEDGE3], nEi[5]] + nF[7, feInds] = [nEi[1], rfZp[1, FEDGE1], nEi[5], rfXp[2, FEDGE3]] + + nF[8, feInds] = [rfYm[0, FEDGE1], nEi[0], rfXm[0, FEDGE1], nEi[2]] + nF[9, feInds] = [rfYm[1, FEDGE1], nEi[1], nEi[2], rfXp[0, FEDGE1]] + nF[10,feInds] = [nEi[0], rfYp[0, FEDGE1], rfXm[2, FEDGE1], nEi[3]] + nF[11,feInds] = [nEi[1], rfYp[1, FEDGE1], nEi[3], rfXp[2, FEDGE1]] + + nFi, nF = self._push('_faces', nF) + + nC = np.zeros((8,9)) + nC[:, ACTIVE] = 1 + nC[:, PARENT] = index + + # .----------------.----------------. + # /| /| /| + # / | / | / | + # / | 6 / | 7 / | + # / | / | / | + # .----------------.----+-----------. | + # /| . ---------/|----.----------/|----. + # fZp / | /| / | /| / | /| + # | / | / | 4 / | / | 5 / | / | + # 6 ------eX3------ 7 / | / | / | / | / | / | + # /| | / | . -------------- .----------------. |/ | + # /eZ2 . / eZ3 | . ---+------|----.----+------|----. | + # eY2 | fYp eY3 | | /| .______|___/|____.______|___/|____. + # / | / fXp| | / | / 2 | / | / 3 | / | / + # 4 ------eX2----- 5 | | / | / | / | / | / | / + # |fXm 2 -----eX1--|---- 3 z . ---+---------- . ---+---------- . | / + # eZ0 / | eY1 ^ y | |/ | |/ | |/ + # | eY0 . fYm eZ1 / | / | . ----------|----.-----------|----. + # | / | | / | / | / 0 | / 1 | / + # 0 ------eX0------1 o----> x | / | / | / + # | | / | / | / + # fZm . -------------- . -------------- . + # + # + # fX fY fZ + # 2___________3 2___________3 2___________3 + # | e1 | | e1 | | e1 | + # | | | | | | + # e2 | x | e3 z e2 | x | e3 z e2 | x | e3 y + # | | ^ | | ^ | | ^ + # |___________| |___> y |___________| |___> x |___________| |___> x + # 0 e0 1 0 e0 1 0 e0 1 + + + cfInds = [CFACE0,CFACE1,CFACE2,CFACE3,CFACE4,CFACE5] + nC[0, cfInds] = [fXm[0], nFi[0], fYm[0], nFi[4], fZm[0], nFi[ 8]] + nC[1, cfInds] = [nFi[0], fXp[0], fYm[1], nFi[5], fZm[1], nFi[ 9]] + nC[2, cfInds] = [fXm[1], nFi[1], nFi[4], fYp[0], fZm[2], nFi[10]] + nC[3, cfInds] = [nFi[1], fXp[1], nFi[5], fYp[1], fZm[3], nFi[11]] + nC[4, cfInds] = [fXm[2], nFi[2], fYm[2], nFi[6], nFi[ 8], fZp[0]] + nC[5, cfInds] = [nFi[2], fXp[2], fYm[3], nFi[7], nFi[ 9], fZp[1]] + nC[6, cfInds] = [fXm[3], nFi[3], nFi[6], fYp[2], nFi[10], fZp[2]] + nC[7, cfInds] = [nFi[3], fXp[3], nFi[7], fYp[3], nFi[11], fZp[3]] + + return self._push('_cells', nC) + + + def _index(self, attr, index): index = [index] if np.isscalar(index) else list(index) @@ -794,7 +899,7 @@ if __name__ == '__main__': print tM._nodes[:,NUM] print tM._edges[:,NUM] - print TreeFace(tM,[0]).e2.n0.x + print TreeFace(tM,0).e3.n1.index @@ -810,10 +915,10 @@ if __name__ == '__main__': # print tM._edges[:,[0,1,3, 4,5 ]] - plt.subplot(211) - plt.spy(tM.faceDiv) - tM.plotGrid(ax=plt.subplot(212)) + # plt.subplot(211) + # plt.spy(tM.faceDiv) + # tM.plotGrid(ax=plt.subplot(212)) # plt.figure(2) # plt.plot(SortByX0(tM.gridCC),'b.') - plt.show() + # plt.show() diff --git a/SimPEG/Tests/test_NewTreeMesh.py b/SimPEG/Tests/test_NewTreeMesh.py index a2a8235f..03cdaba6 100644 --- a/SimPEG/Tests/test_NewTreeMesh.py +++ b/SimPEG/Tests/test_NewTreeMesh.py @@ -108,9 +108,9 @@ class TestOcTreeObjects(unittest.TestCase): self.M = TreeMesh([2,1,1]) self.M.number() - # self.Mr = TreeMesh([2,1,1]) - # self.Mr.children[0,0,0].refine() - # self.Mr.number() + self.Mr = TreeMesh([2,1,1]) + self.Mr.refineCell(0) + self.Mr.number() def test_counts(self): self.assertTrue(self.M.nC == 2) @@ -124,22 +124,106 @@ class TestOcTreeObjects(unittest.TestCase): self.assertTrue(self.M.nE == 20) self.assertTrue(self.M.nN == 12) - # self.assertTrue(self.Mr.nC == 9) - # self.assertTrue(self.Mr.nFx == 13) - # self.assertTrue(self.Mr.nFy == 14) - # self.assertTrue(self.Mr.nFz == 14) - # self.assertTrue(self.Mr.nF == 41) + self.assertTrue(self.Mr.nC == 9) + self.assertTrue(self.Mr.nFx == 13) + self.assertTrue(self.Mr.nFy == 14) + self.assertTrue(self.Mr.nFz == 14) + self.assertTrue(self.Mr.nF == 41) + + self.assertTrue(self.Mr.nN == 31) + self.assertTrue(self.Mr.nEx == 22) + self.assertTrue(self.Mr.nEy == 20) + self.assertTrue(self.Mr.nEz == 20) - # for cell in self.Mr.sortedCells: - # for e in cell.edgeDict: - # self.assertTrue(cell.edgeDict[e].edgeType==e[1].lower()) + def test_gridCC(self): + x = np.r_[0.25,0.75] + y = np.r_[0.5,0.5] + z = np.r_[0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridCC).flatten()) == 0) - # self.assertTrue(self.Mr.nN == 31) - # self.assertTrue(self.Mr.nEx == 22) - # self.assertTrue(self.Mr.nEy == 20) - # self.assertTrue(self.Mr.nEz == 20) + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375] + y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75] + z = np.r_[0.25,0.25,0.5,0.25,0.25,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridCC).flatten()) == 0) + def test_gridN(self): + x = np.r_[0,0.5,1,0,0.5,1,0,0.5,1,0,0.5,1] + y = np.r_[0,0,0,1,1,1,0,0,0,1,1,1.] + z = np.r_[0,0,0,0,0,0,1,1,1,1,1,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridN).flatten()) == 0) + + x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1] + y = np.r_[0,0,0,0,0.5,0.5,0.5,1,1,1,1,0,0,0,0.5,0.5,0.5,1,1,1,0,0,0,0,0.5,0.5,0.5,1,1,1,1] + z = np.r_[0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridN).flatten()) == 0) + + def test_gridFx(self): + x = np.r_[0.0,0.5,1.0] + y = np.r_[0.5,0.5,0.5] + z = np.r_[0.5,0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFx).flatten()) == 0) + + x = np.r_[0.0,0.25,0.5,1.0,0.0,0.25,0.5,0.0,0.25,0.5,0.0,0.25,0.5] + y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75] + z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.75,0.75,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFx).flatten()) == 0) + + def test_gridFy(self): + x = np.r_[0.25,0.75,0.25,0.75] + y = np.r_[0,0,1.,1.] + z = np.r_[0.5,0.5,0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFy).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375] + y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1] + z = np.r_[0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFy).flatten()) == 0) + + def test_gridFz(self): + x = np.r_[0.25,0.75,0.25,0.75] + y = np.r_[0.5,0.5,0.5,0.5] + z = np.r_[0,0,1.,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFz).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375] + y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75,0.25,0.25,0.5,0.75,0.75] + z = np.r_[0,0,0,0,0,0.5,0.5,0.5,0.5,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFz).flatten()) == 0) + + + def test_gridEx(self): + x = np.r_[0.25,0.75,0.25,0.75,0.25,0.75,0.25,0.75] + y = np.r_[0,0,1.,1.,0,0,1.,1.] + z = np.r_[0,0,0,0,1.,1.,1.,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEx).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75] + y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1,0,0,0,0.5,0.5,1,1,1] + z = np.r_[0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEx).flatten()) == 0) + + def test_gridEy(self): + x = np.r_[0,0.5,1,0,0.5,1] + y = np.r_[0.5,0.5,0.5,0.5,0.5,0.5] + z = np.r_[0,0,0,1.,1.,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEy).flatten()) == 0) + + x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5] + y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75,0.25,0.25,0.25,0.5,0.75,0.75,0.75] + z = np.r_[0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEy).flatten()) == 0) + + def test_gridEz(self): + x = np.r_[0,0.5,1,0,0.5,1] + y = np.r_[0,0,0,1.,1.,1.] + z = np.r_[0.5,0.5,0.5,0.5,0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEz).flatten()) == 0) + + x = np.r_[0,0.25,0.5,1,0 ,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0 ,0.25,0.5,0 ,0.25,0.5] + y = np.r_[0,0 ,0 ,0,0.5,0.5 ,0.5,1,1 ,1 ,1,0,0 ,0 ,0.5,0.5 ,0.5,1 ,1 ,1 ] + z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEz).flatten()) == 0) if __name__ == '__main__': unittest.main()