mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-22 13:00:21 +08:00
Testing and Face numbering (shared faces are shared!)
This commit is contained in:
rowanc1
+1
-1
@@ -6,7 +6,7 @@ import mesh
|
||||
import data
|
||||
import forward
|
||||
import inverse
|
||||
import visualize
|
||||
# import visualize
|
||||
import examples
|
||||
|
||||
import scipy.version as _v
|
||||
|
||||
+97
-60
@@ -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()
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user