From f25d5bd524f0391e1f711e2c91648714ac1099b5 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 19 Dec 2013 14:25:34 -0800 Subject: [PATCH] Testing and Face numbering (shared faces are shared!) --- SimPEG/__init__.py | 2 +- SimPEG/mesh/QuadTreeMesh.py | 157 +++++++++++++++++++++------------- SimPEG/mesh/__init__.py | 1 + SimPEG/tests/test_QuadTree.py | 21 +++++ 4 files changed, 120 insertions(+), 61 deletions(-) create mode 100644 SimPEG/tests/test_QuadTree.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index f421bc82..70f29c46 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -6,7 +6,7 @@ import mesh import data import forward import inverse -import visualize +# import visualize import examples import scipy.version as _v diff --git a/SimPEG/mesh/QuadTreeMesh.py b/SimPEG/mesh/QuadTreeMesh.py index 5d5c6f41..934147d9 100644 --- a/SimPEG/mesh/QuadTreeMesh.py +++ b/SimPEG/mesh/QuadTreeMesh.py @@ -2,16 +2,45 @@ import numpy as np import matplotlib.pyplot as plt + +def SortByX0(): + eps = 1e-7 + def mycmp(c1,c2): + if np.abs(c1.x0[1] - c2.x0[1]) < eps: + return c1.x0[0] - c2.x0[0] + return c1.x0[1] - c2.x0[1] + class K(object): + def __init__(self, obj, *args): + self.obj = obj + def __lt__(self, other): + return mycmp(self.obj, other.obj) < 0 + def __gt__(self, other): + return mycmp(self.obj, other.obj) > 0 + def __eq__(self, other): + return mycmp(self.obj, other.obj) == 0 + def __le__(self, other): + return mycmp(self.obj, other.obj) <= 0 + def __ge__(self, other): + return mycmp(self.obj, other.obj) >= 0 + def __ne__(self, other): + return mycmp(self.obj, other.obj) != 0 + return K + + class TreeFace(object): """docstring for TreeFace""" def __init__(self, mesh, x0=[0,0], faceType=None, dim=2, sz=1, depth=0, parent=None): self.mesh = mesh self.children = None + self.numFace = None self.x0 = np.array(x0,dtype=float) self.faceType = faceType self.sz = sz self.dim = dim + self.depth = depth mesh.faces.add(self) + if faceType is 'x': self.mesh.faceX.add(self) + elif faceType is 'y': self.mesh.faceY.add(self) self.tangent = np.zeros(dim) self.tangent[1 if faceType is 'x' else 0] = 1 self.normal = np.zeros(dim) @@ -25,30 +54,40 @@ class TreeFace(object): self.children = np.empty(2,dtype=TreeFace) # Create refined x0's x0r_0 = self.x0 - x0r_1 = self.x0+self.tangent*self.sz/2 - self.children[0] = TreeFace(self.mesh, x0=x0r_0, faceType=self.faceType, dim=self.dim, sz=self.sz/2, depth=self.depth+1,parent=self) - self.children[1] = TreeFace(self.mesh, x0=x0r_1, faceType=self.faceType, dim=self.dim, sz=self.sz/2, depth=self.depth+1,parent=self) + x0r_1 = self.x0+0.5*self.tangent*self.sz + self.children[0] = TreeFace(self.mesh, x0=x0r_0, faceType=self.faceType, dim=self.dim, sz=0.5*self.sz, depth=self.depth+1,parent=self) + self.children[1] = TreeFace(self.mesh, x0=x0r_1, faceType=self.faceType, dim=self.dim, sz=0.5*self.sz, depth=self.depth+1,parent=self) self.mesh.faces.remove(self) + if self.faceType is 'x': + self.mesh.faceX.remove(self) + elif self.faceType is 'y': + self.mesh.faceY.remove(self) def viz(self, ax): + if not self.isleaf: return ax.plot(np.r_[self.x0[0],self.x0[0]+self.tangent[0]*self.sz], np.r_[self.x0[1], self.x0[1]+self.tangent[1]*self.sz],'rx-') + if self.faceType is 'y': + ax.text(self.x0[0]+0.5*self.tangent[0]*self.sz, self.x0[1]+0.5*self.tangent[1]*self.sz,self.numFace) + class TreeNode(object): """docstring for TreeNode""" def __init__(self, mesh, x0=[0,0], dim=2, depth=0, sz=[1,1], parent=None, fXm=None, fXp=None, fYm=None, fYp=None): - self.fXm = fXm if fXm is not None else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='x', dim=dim, sz=sz[1], depth=depth, parent=parent) - self.fXp = fXp if fXp is not None else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] ], faceType='x', dim=dim, sz=sz[1], depth=depth, parent=parent) - self.fYm = fYm if fYm is not None else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='y', dim=dim, sz=sz[0], depth=depth, parent=parent) - self.fYp = fYp if fYp is not None else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1]], faceType='y', dim=dim, sz=sz[0], depth=depth, parent=parent) + fXm = fXm if fXm is not None else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='x', dim=dim, sz=sz[1], depth=depth, parent=parent) + fXp = fXp if fXp is not None else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] ], faceType='x', dim=dim, sz=sz[1], depth=depth, parent=parent) + fYm = fYm if fYm is not None else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='y', dim=dim, sz=sz[0], depth=depth, parent=parent) + fYp = fYp if fYp is not None else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1]], faceType='y', dim=dim, sz=sz[0], depth=depth, parent=parent) + + self.faces = {"fXm":fXm, "fXp":fXp, "fYm":fYm, "fYp":fYp} self.mesh = mesh - self.x0 = np.array(x0,dtype=float) + self.x0 = np.array(x0, dtype=float) self.dim = dim self.depth = depth - self.sz = np.array(sz,dtype=float) + self.sz = np.array(sz, dtype=float) self.parent = parent self.children = None self.numCell = None @@ -72,17 +111,33 @@ class TreeNode(object): def refine(self): if not self.isleaf: return + self.children = np.empty((2,2),dtype=TreeNode) x0, sz = self.x0, self.sz - corners = [np.r_[x0[0] , x0[1] ], - np.r_[x0[0]+sz[0]/2, x0[1] ], - np.r_[x0[0] , x0[1]+sz[1]/2], - np.r_[x0[0]+sz[0]/2, x0[1]+sz[1]/2]] - children = [TreeNode(self.mesh, x0=corners[i],dim=self.dim, depth=self.depth+1, sz=sz/2, parent=self) for i in range(2*self.dim)] + for faceName in self.faces: + self.faces[faceName].refine() + + i, j = 0, 0 + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1]] + fXm, fXp, fYm, fYp = self.faces['fXm'].children[0], None, self.faces['fYm'].children[0], None + self.children[i,j] = TreeNode(self.mesh, x0=x0r,dim=self.dim, depth=self.depth+1, sz=0.5*sz, parent=self, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp) + + i, j = 1, 0 + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1]] + fXm, fXp, fYm, fYp = self.children[0,0].faces['fXp'], self.faces['fXp'].children[0], self.faces['fYm'].children[1], None + self.children[i,j] = TreeNode(self.mesh, x0=x0r,dim=self.dim, depth=self.depth+1, sz=0.5*sz, parent=self, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp) + + i, j = 0, 1 + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1]] + fXm, fXp, fYm, fYp = self.faces['fXm'].children[1], None, self.children[0,0].faces['fYp'], self.faces['fYp'].children[0] + self.children[i,j] = TreeNode(self.mesh, x0=x0r,dim=self.dim, depth=self.depth+1, sz=0.5*sz, parent=self, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp) + + i, j = 1, 1 + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1]] + fXm, fXp, fYm, fYp = self.children[0,1].faces['fXp'], self.faces['fXp'].children[1], self.children[1,0].faces['fYp'], self.faces['fYp'].children[1] + self.children[i,j] = TreeNode(self.mesh, x0=x0r,dim=self.dim, depth=self.depth+1, sz=0.5*sz, parent=self, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp) + self.mesh.cells.remove(self) - self.mesh.cells = self.mesh.cells.union(set(children)) - self.children = np.array(children,dtype=TreeNode).reshape([2,2],order='F') - def countCell(self, COUNTER, column, maxColumns): @@ -115,7 +170,6 @@ class TreeNode(object): np.r_[x0[0] , x0[1]+sz[1]], np.r_[x0[0] , x0[1] ]].T ax.plot(corners[:,0],corners[:,1], 'b') - self.fXm.viz(ax) if self.numCell is not None: ax.text(x0[0]+sz[0]/2,x0[1]+sz[1]/2,'%d'%self.numCell) else: @@ -128,13 +182,15 @@ class QuadTreeMesh(object): """docstring for QuadTreeMesh""" def __init__(self, cells, sz): self.faces = set() + self.faceX = set() + self.faceY = set() self.cells = set() self.children = np.empty(cells,dtype=TreeNode) for i in range(cells[0]): for j in range(cells[1]): - fXm = None if i is 0 else self.children[i-1][j].fXp - fYm = None if j is 0 else self.children[i][j-1].fYp + fXm = None if i is 0 else self.children[i-1][j].faces['fXp'] + fYm = None if j is 0 else self.children[i][j-1].faces['fYp'] self.children[i][j] = TreeNode(self, x0=[i*sz[0],j*sz[1]], dim=2, depth=0, sz=sz, fXm=fXm, fYm=fYm) @@ -144,13 +200,21 @@ class QuadTreeMesh(object): def number(self): - COUNTER = 0 - for col in range(self.children.shape[1]): - coldepth = np.max([node.branchdepth for node in self.children[:,col]]) - maxColumns = 2**coldepth - for COL in range(maxColumns): - for row in range(self.children.shape[0]): - COUNTER = self.children[row,col].countCell(COUNTER, COL, maxColumns) + # COUNTER = 0 + # for col in range(self.children.shape[1]): + # coldepth = np.max([node.branchdepth for node in self.children[:,col]]) + # maxColumns = 2**coldepth + # for COL in range(maxColumns): + # for row in range(self.children.shape[0]): + # COUNTER = self.children[row,col].countCell(COUNTER, COL, maxColumns) + + sortCells = sorted(M.cells,key=SortByX0()) + sortFaceX = sorted(M.faceX,key=SortByX0()) + sortFaceY = sorted(M.faceY,key=SortByX0()) + nFx = len(sortFaceX) + for i, sc in enumerate(sortCells): sc.numCell = i + for i, sfx in enumerate(sortFaceX): sfx.numFace = i + for i, sfy in enumerate(sortFaceY): sfy.numFace = i + nFx def viz(self): ax = plt.subplot(111) @@ -159,45 +223,18 @@ class QuadTreeMesh(object): if __name__ == '__main__': - M = QuadTreeMesh([4,6],[1,2]) - M.children[0,0].refine() - M.children[1,1].refine() - for ii in range(3): + M = QuadTreeMesh([3,2],[1,2]) + for ii in range(1): M.children[ii,ii].refine() - M.children[ii,ii].children[0,0].refine() + # M.children[ii,ii].children[0,0].refine() - eps = 1e-7 - def mycmp(c1,c2): - if np.abs(c1.x0[1] - c2.x0[1]) < eps: - return c1.x0[0] - c2.x0[0] - return c1.x0[1] - c2.x0[1] - def cmp_to_key(mycmp): - 'Convert a cmp= function into a key= function' - class K(object): - def __init__(self, obj, *args): - self.obj = obj - def __lt__(self, other): - return mycmp(self.obj, other.obj) < 0 - def __gt__(self, other): - return mycmp(self.obj, other.obj) > 0 - def __eq__(self, other): - return mycmp(self.obj, other.obj) == 0 - def __le__(self, other): - return mycmp(self.obj, other.obj) <= 0 - def __ge__(self, other): - return mycmp(self.obj, other.obj) >= 0 - def __ne__(self, other): - return mycmp(self.obj, other.obj) != 0 - return K - # M.number() - sortCells = sorted(M.cells,key=cmp_to_key(mycmp)) - for i, sc in enumerate(sortCells): - sc.numCell = i + M.number() - print M.children[0,0].fXp is M.children[1,0].fXm + print M.children[0,0].faces['fXp'] is M.children[1,0].faces['fXm'] + print len(M.faces) M.viz() # plt.gca().invert_yaxis() plt.show() diff --git a/SimPEG/mesh/__init__.py b/SimPEG/mesh/__init__.py index 3b8e1eef..66ee88c8 100644 --- a/SimPEG/mesh/__init__.py +++ b/SimPEG/mesh/__init__.py @@ -1,5 +1,6 @@ from Cyl1DMesh import Cyl1DMesh from TensorMesh import TensorMesh +from QuadTreeMesh import QuadTreeMesh from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh from BaseMesh import BaseMesh from TensorView import TensorView diff --git a/SimPEG/tests/test_QuadTree.py b/SimPEG/tests/test_QuadTree.py new file mode 100644 index 00000000..c085c7ae --- /dev/null +++ b/SimPEG/tests/test_QuadTree.py @@ -0,0 +1,21 @@ +from SimPEG import mesh, np +import unittest + + + +class TestCheckDerivative(unittest.TestCase): + + def setUp(self): + M = mesh.QuadTreeMesh([3,2],[1,2]) + for ii in range(1): + M.children[ii,ii].refine() + self.M = M + + def test_MeshSizes(self): + self.assertTrue(len(self.M.faces)==25) + self.assertTrue(len(self.M.cells)==9) + + + +if __name__ == '__main__': + unittest.main()