mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-04 21:04:57 +08:00
Switch to Edges (few bugs in the refinement code..)
This commit is contained in:
rowanc1
+176
-73
@@ -86,14 +86,44 @@ class TreeEdge(TreeObject):
|
||||
elif edgeType is 'y': mesh.edgesY.add(self)
|
||||
elif edgeType is 'z': mesh.edgesZ.add(self)
|
||||
|
||||
self.node0 = node0
|
||||
self.node1 = node1
|
||||
self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=self.x0)
|
||||
self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent*self.sz[0])
|
||||
|
||||
|
||||
def refine(self):
|
||||
if not self.isleaf: return
|
||||
self.mesh.isNumbered = False
|
||||
|
||||
self.children = np.empty(2,dtype=TreeFace)
|
||||
# Create refined x0's
|
||||
x0r_0 = self.x0
|
||||
x0r_1 = self.x0+0.5*self.tangent*self.sz
|
||||
self.children[0] = TreeEdge(self.mesh, x0=x0r_0, edgeType=self.edgeType, sz=0.5*self.sz, depth=self.depth+1, parent=self, node0=self.node0)
|
||||
self.children[1] = TreeEdge(self.mesh, x0=x0r_1, edgeType=self.edgeType, sz=0.5*self.sz, depth=self.depth+1, parent=self, node0=self.children[0].node1, node1=self.node1)
|
||||
self.mesh.edges.remove(self)
|
||||
if self.edgeType is 'x':
|
||||
self.mesh.edgesX.remove(self)
|
||||
elif self.edgeType is 'y':
|
||||
self.mesh.edgesY.remove(self)
|
||||
elif self.edgeType is 'z':
|
||||
self.mesh.edgesZ.remove(self)
|
||||
|
||||
@property
|
||||
def tangent(self):
|
||||
if self.edgeType is 'x': return np.r_[1.,0,0]
|
||||
elif self.edgeType is 'y': return np.r_[0,1.,0]
|
||||
elif self.edgeType is 'z': return np.r_[0,0,1.]
|
||||
|
||||
def plotGrid(self, ax, text=False):
|
||||
line = np.c_[self.node0.x0, self.node1.x0].T
|
||||
ax.plot(line[:,0], line[:,1],'r-', zs=line[:,2])
|
||||
|
||||
|
||||
class TreeFace(TreeObject):
|
||||
"""docstring for TreeFace"""
|
||||
def __init__(self, mesh, x0=[0,0], faceType=None, sz=[1,], depth=0,
|
||||
node0=None, node1=None, node2=None, node3=None,
|
||||
node0=None, node1=None,
|
||||
edge0=None, edge1=None, edge2=None, edge3=None,
|
||||
parent=None):
|
||||
TreeObject.__init__(self, mesh, parent)
|
||||
|
||||
@@ -106,12 +136,46 @@ class TreeFace(TreeObject):
|
||||
if faceType is 'x': mesh.facesX.add(self)
|
||||
elif faceType is 'y': mesh.facesY.add(self)
|
||||
elif faceType is 'z': mesh.facesZ.add(self)
|
||||
# Add the nodes:
|
||||
self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=self.x0)
|
||||
self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent0*self.sz[0])
|
||||
if self.dim == 2:
|
||||
# Add the nodes:
|
||||
self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=self.x0)
|
||||
self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent0*self.sz[0])
|
||||
if self.dim == 3:
|
||||
self.node2 = node2 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent1*self.sz[1])
|
||||
self.node3 = node3 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent0*self.sz[0] + self.tangent1*self.sz[1])
|
||||
#TODO: Change this to edges
|
||||
|
||||
#
|
||||
# 2___________3
|
||||
# | e1 |
|
||||
# | |
|
||||
# e2| x |e3 t1
|
||||
# | | ^
|
||||
# |___________| |___> t0
|
||||
# 0 e0 1
|
||||
#
|
||||
|
||||
N = {}
|
||||
N["n0"] = getattr(edge0, 'node0', None) or getattr(edge2, 'node0', None)
|
||||
N["n1"] = getattr(edge0, 'node1', None) or getattr(edge3, 'node0', None)
|
||||
N["n2"] = getattr(edge1, 'node0', None) or getattr(edge2, 'node1', None)
|
||||
N["n3"] = getattr(edge1, 'node1', None) or getattr(edge3, 'node1', None)
|
||||
|
||||
eType = ['x', 'y'] if self.faceType == 'z' else ['x', 'z'] if self.faceType == 'y' else ['y', 'z']
|
||||
|
||||
e0 = edge0 if isinstance(edge0,TreeEdge) else TreeEdge(mesh, x0=self.x0, edgeType=eType[0], sz=np.r_[sz[0]], depth=depth, parent=parent, node0=N['n0'], node1=N['n1'])
|
||||
N["n0"], N["n1"] = e0.node0, e0.node1
|
||||
|
||||
e1 = edge1 if isinstance(edge1,TreeEdge) else TreeEdge(mesh, x0=self.x0 + self.tangent1*self.sz[1], edgeType=eType[0], sz=np.r_[sz[0]], depth=depth, parent=parent, node0=N['n2'], node1=N['n3'])
|
||||
N["n2"], N["n3"] = e1.node0, e1.node1
|
||||
|
||||
e2 = edge2 if isinstance(edge2,TreeEdge) else TreeEdge(mesh, x0=self.x0, edgeType=eType[1], sz=np.r_[sz[1]], depth=depth, parent=parent, node0=N['n0'], node1=N['n2'])
|
||||
N["n0"], N["n2"] = e2.node0, e2.node1
|
||||
|
||||
e3 = edge3 if isinstance(edge3,TreeEdge) else TreeEdge(mesh, x0=self.x0 + self.tangent0*self.sz[0], edgeType=eType[1], sz=np.r_[sz[1]], depth=depth, parent=parent, node0=N['n1'], node1=N['n3'])
|
||||
N["n1"], N["n3"] = e3.node0, e3.node1
|
||||
|
||||
self.nodes = N
|
||||
self.edges = {'e0':e0, 'e1':e1, 'e2':e2, 'e3':e3}
|
||||
|
||||
|
||||
@property
|
||||
def tangent0(self):
|
||||
@@ -169,31 +233,48 @@ class TreeFace(TreeObject):
|
||||
|
||||
def _refine3D(self):
|
||||
self.children = np.empty((2,2),dtype=TreeFace)
|
||||
# Create refined x0's
|
||||
x0r_0 = self.x0
|
||||
x0r_1 = self.x0+0.5*self.tangent0*self.sz[0]
|
||||
x0r_2 = self.x0+0.5*self.tangent1*self.sz[1]
|
||||
x0r_3 = self.x0+0.5*self.tangent0*self.sz[0]+0.5*self.tangent1*self.sz[1]
|
||||
|
||||
for edgeName in self.edges:
|
||||
self.edges[edgeName].refine()
|
||||
|
||||
#
|
||||
# 2_______________3 _______________
|
||||
# | | | | |
|
||||
# ^ | | | (0,1) | (1,1) |
|
||||
# | | | | | |
|
||||
# | | x | ---> |-------+-------|
|
||||
# t1 | | | | |
|
||||
# | e1--> | | | |
|
||||
# ^ | | ^ | (0,1) | (1,1) |
|
||||
# | | | | | | |
|
||||
# | | x | | ---> |-------+-------|
|
||||
# e2 | | e3 | | |
|
||||
# | | | (0,0) | (1,0) |
|
||||
# |_______________| |_______|_______|
|
||||
# 0 t0--> 1
|
||||
# 0 e0--> 1
|
||||
|
||||
c00 = TreeFace(self.mesh, x0=x0r_0, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, parent=self, node0=self.node0)
|
||||
c01 = TreeFace(self.mesh, x0=x0r_1, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, parent=self, node0=c00.node2, node1=c00.node3, node2=self.node2)
|
||||
c10 = TreeFace(self.mesh, x0=x0r_2, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, parent=self, node0=c00.node1, node1=self.node1, node2=c00.node3)
|
||||
c11 = TreeFace(self.mesh, x0=x0r_3, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, parent=self, node0=c00.node3, node1=c10.node3, node2=c01.node3, node3=self.node3)
|
||||
|
||||
C = np.empty((2,2), dtype=TreeFace)
|
||||
C[0,0], C[0,1], C[1,0], C[1,1] = c00, c01, c10, c11
|
||||
self.children = C
|
||||
order = [{'c':[0,0],
|
||||
'e0': ('p', 'e0', [0]), 'e1': 'new' ,
|
||||
'e2': ('p', 'e2', [0]), 'e3': 'new' },
|
||||
{'c':[1,0],
|
||||
'e0': ('p', 'e0', [1]), 'e1': 'new' ,
|
||||
'e2': ('c', 'e3', [0,0]), 'e3': ('p', 'e3', [0])},
|
||||
{'c':[0,1],
|
||||
'e0': ('c', 'e1', [0,0]), 'e1': ('p', 'e1', [0]),
|
||||
'e2': ('p', 'e2', [1]), 'e3': 'new' },
|
||||
{'c':[1,1],
|
||||
'e0': ('c', 'e1', [1,0]), 'e1': ('p', 'e1', [1]),
|
||||
'e2': ('c', 'e3', [0,1]), 'e3': ('p', 'e3', [1])}]
|
||||
|
||||
def getEdge(pointer):
|
||||
return
|
||||
if pointer is 'new': return
|
||||
if pointer[0] == 'p':
|
||||
return self.edges[pointer[1]].children[pointer[2][0]]
|
||||
if pointer[0] == 'c':
|
||||
return self.children[pointer[2][0],pointer[2][1]].edges[pointer[1]]
|
||||
|
||||
for O in order:
|
||||
i, j = O['c']
|
||||
x0r = 0.5*i*self.tangent0*self.sz[0] + 0.5*j*self.tangent1*self.sz[1]
|
||||
e0, e1, e2, e3 = getEdge(O['e0']), getEdge(O['e1']), getEdge(O['e2']), getEdge(O['e3'])
|
||||
self.children[i,j] = TreeFace(self.mesh, x0=x0r, faceType=self.faceType, depth=self.depth+1, sz=0.5*self.sz, parent=self, edge0=e0, edge1=e1, edge2=e2, edge3=e3)
|
||||
|
||||
self.mesh.faces.remove(self)
|
||||
if self.faceType is 'x':
|
||||
@@ -210,8 +291,6 @@ class TreeFace(TreeObject):
|
||||
ax.plot(line[:,0], line[:,1],'r-')
|
||||
if text: ax.text(self.center[0], self.center[1],self.num)
|
||||
elif self.dim == 3:
|
||||
line = np.c_[self.node0.x0, self.node1.x0, self.node3.x0, self.node2.x0, self.node0.x0].T
|
||||
ax.plot(line[:,0], line[:,1],'r-', zs=line[:,2])
|
||||
if text: ax.text(self.center[0], self.center[1], self.center[2], self.num)
|
||||
|
||||
@property
|
||||
@@ -270,65 +349,76 @@ class TreeCell(TreeObject):
|
||||
elif self.dim == 3:
|
||||
# fZp
|
||||
# |
|
||||
# 6 --------------- 7
|
||||
# 6 ------eX3------ 7
|
||||
# /| | / |
|
||||
# / | . / |
|
||||
# / | fYp / |
|
||||
# /eZ2 . / eZ3
|
||||
# eY2 | fYp eY3 |
|
||||
# / | / fXp|
|
||||
# 4 -------------- 5 |
|
||||
# |fXm 2 ----------|---- 3 z
|
||||
# | / | / ^ y
|
||||
# | / fYm . | / | /
|
||||
# | / | | / | /
|
||||
# 0 -------------- 1 o----> x
|
||||
# |
|
||||
# fZm
|
||||
# 4 ------eX2----- 5 |
|
||||
# |fXm 2 -----eX1--|---- 3 z
|
||||
# eZ0 / | eY1 ^ y
|
||||
# | eY0 . fYm eZ1 / | /
|
||||
# | / | | / | /
|
||||
# 0 ------eX0------1 o----> x
|
||||
# |
|
||||
# fZm
|
||||
#
|
||||
#
|
||||
# fX fY fZ
|
||||
# 2___________3 2___________3 2___________3
|
||||
# | eYp | | eXp | | eXp |
|
||||
# | e1 | | e1 | | e1 |
|
||||
# | | | | | |
|
||||
# eZm| x |eZp z eZm| x |eZp z eYm| x |eYp y
|
||||
# e2 | x | e3 z e2 | x | e3 z e2 | x | e3 y
|
||||
# | | ^ | | ^ | | ^
|
||||
# |___________| |___> y |___________| |___> x |___________| |___> x
|
||||
# 0 eYm 1 0 eXm 1 0 eXm 1
|
||||
# 0 e0 1 0 e0 1 0 e0 1
|
||||
#
|
||||
# Mapping Nodes: numOnFace > numOnCell
|
||||
#
|
||||
# fXm 0>0, 1>2, 2>4, 3>6 fYm 0>0, 1>1, 2>4, 3>5 fZm 0>0, 1>1, 2>2, 3>3
|
||||
# fXp 0>1, 1>3, 2>5, 3>7 fYp 0>2, 1>3, 2>6, 3>7 fZp 0>4, 1>5, 2>6, 3>7
|
||||
|
||||
N = {}
|
||||
N["n0"] = getattr(fXm, 'node0', None) or getattr(fYm, 'node0', None) or getattr(fZm, 'node0', None)
|
||||
N["n1"] = getattr(fXp, 'node0', None) or getattr(fYm, 'node1', None) or getattr(fZm, 'node1', None)
|
||||
N["n2"] = getattr(fXm, 'node1', None) or getattr(fYp, 'node0', None) or getattr(fZm, 'node2', None)
|
||||
N["n3"] = getattr(fXp, 'node1', None) or getattr(fYp, 'node1', None) or getattr(fZm, 'node3', None)
|
||||
N["n4"] = getattr(fXm, 'node2', None) or getattr(fYm, 'node2', None) or getattr(fZp, 'node0', None)
|
||||
N["n5"] = getattr(fXp, 'node2', None) or getattr(fYm, 'node3', None) or getattr(fZp, 'node1', None)
|
||||
N["n6"] = getattr(fXm, 'node3', None) or getattr(fYp, 'node2', None) or getattr(fZp, 'node2', None)
|
||||
N["n7"] = getattr(fXp, 'node3', None) or getattr(fYp, 'node3', None) or getattr(fZp, 'node3', None)
|
||||
def getEdge(face, key):
|
||||
if face is None: return
|
||||
return face.edges[key]
|
||||
|
||||
E = {}
|
||||
eX0 = getEdge(fYm, 'e0') or getEdge(fZm, 'e0')
|
||||
eX1 = getEdge(fYp, 'e0') or getEdge(fZm, 'e1')
|
||||
eX2 = getEdge(fYm, 'e1') or getEdge(fZp, 'e0')
|
||||
eX3 = getEdge(fYp, 'e1') or getEdge(fZp, 'e1')
|
||||
|
||||
eY0 = getEdge(fXm, 'e0') or getEdge(fZm, 'e2')
|
||||
eY1 = getEdge(fXp, 'e0') or getEdge(fZm, 'e3')
|
||||
eY2 = getEdge(fXm, 'e1') or getEdge(fZp, 'e2')
|
||||
eY3 = getEdge(fXp, 'e1') or getEdge(fZp, 'e3')
|
||||
|
||||
eZ0 = getEdge(fXm, 'e2') or getEdge(fYm, 'e2')
|
||||
eZ1 = getEdge(fXp, 'e2') or getEdge(fYm, 'e3')
|
||||
eZ2 = getEdge(fXm, 'e3') or getEdge(fYp, 'e2')
|
||||
eZ3 = getEdge(fXp, 'e3') or getEdge(fYp, 'e3')
|
||||
|
||||
|
||||
fXm = fXm if isinstance(fXm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, parent=parent, node0=N["n0"], node1=N["n2"], node2=N["n4"], node3=N["n6"])
|
||||
N["n0"], N["n2"], N["n4"], N["n6"] = fXm.node0, fXm.node1, fXm.node2, fXm.node3
|
||||
fXm = fXm if isinstance(fXm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, parent=parent, edge0=eY0, edge1=eY2, edge2=eZ0, edge3=eZ2)
|
||||
eY0, eY2, eZ0, eZ2 = fXm.edges['e0'], fXm.edges['e1'], fXm.edges['e2'], fXm.edges['e3']
|
||||
|
||||
fXp = fXp if isinstance(fXp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, parent=parent, node0=N["n1"], node1=N["n3"], node2=N["n5"], node3=N["n7"])
|
||||
N["n1"], N["n3"], N["n5"], N["n7"] = fXp.node0, fXp.node1, fXp.node2, fXp.node3
|
||||
fXp = fXp if isinstance(fXp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, parent=parent, edge0=eY1, edge1=eY3, edge2=eZ1, edge3=eZ3)
|
||||
eY1, eY3, eZ1, eZ3 = fXp.edges['e0'], fXp.edges['e1'], fXp.edges['e2'], fXp.edges['e3']
|
||||
|
||||
fYm = fYm if isinstance(fYm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, parent=parent, node0=N["n0"], node1=N["n1"], node2=N["n4"], node3=N["n5"])
|
||||
N["n0"], N["n1"], N["n4"], N["n5"] = fYm.node0, fYm.node1, fYm.node2, fYm.node3
|
||||
fYm = fYm if isinstance(fYm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, parent=parent, edge0=eX0, edge1=eX2, edge2=eZ0, edge3=eZ1)
|
||||
eX0, eX2, eZ0, eZ1 = fYm.edges['e0'], fYm.edges['e1'], fYm.edges['e2'], fYm.edges['e3']
|
||||
|
||||
fYp = fYp if isinstance(fYp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1], x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, parent=parent, node0=N["n2"], node1=N["n3"], node2=N["n6"], node3=N["n7"])
|
||||
N["n2"], N["n3"], N["n6"], N["n7"] = fYp.node0, fYp.node1, fYp.node2, fYp.node3
|
||||
fYp = fYp if isinstance(fYp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1], x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, parent=parent, edge0=eX1, edge1=eX3, edge2=eZ2, edge3=eZ3)
|
||||
eX1, eX3, eZ2, eZ3 = fYp.edges['e0'], fYp.edges['e1'], fYp.edges['e2'], fYp.edges['e3']
|
||||
|
||||
fZm = fZm if isinstance(fZm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, parent=parent, node0=N["n0"], node1=N["n1"], node2=N["n2"], node3=N["n3"])
|
||||
N["n0"], N["n1"], N["n2"], N["n3"] = fZm.node0, fZm.node1, fZm.node2, fZm.node3
|
||||
fZm = fZm if isinstance(fZm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, parent=parent, edge0=eX0, edge1=eX1, edge2=eY0, edge3=eY1)
|
||||
eX0, eX1, eY0, eY1 = fZm.edges['e0'], fZm.edges['e1'], fZm.edges['e2'], fZm.edges['e3']
|
||||
|
||||
fZp = fZp if isinstance(fZp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2]+sz[2]], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, parent=parent, node0=N["n4"], node1=N["n5"], node2=N["n6"], node3=N["n7"])
|
||||
N["n4"], N["n5"], N["n6"], N["n7"] = fZp.node0, fZp.node1, fZp.node2, fZp.node3
|
||||
fZp = fZp if isinstance(fZp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2]+sz[2]], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, parent=parent, edge0=eX2, edge1=eX3, edge2=eY2, edge3=eY3)
|
||||
eX2, eX3, eY2, eY3 = fZp.edges['e0'], fZp.edges['e1'], fZp.edges['e2'], fZp.edges['e3']
|
||||
|
||||
self.faces = {"fXm":fXm, "fXp":fXp, "fYm":fYm, "fYp":fYp, "fZm":fZm, "fZp":fZp}
|
||||
self.edges = {"eX0":eX0, "eX1":eX1, "eX2":eX2, "eX3":eX3, "eY0":eY0, "eY1":eY1, "eY2":eY2, "eY3":eY3, "eZ0":eZ0, "eZ1":eZ1, "eZ2":eZ2, "eZ3":eZ3}
|
||||
|
||||
mesh.cells.add(self)
|
||||
|
||||
@@ -350,9 +440,9 @@ class TreeCell(TreeObject):
|
||||
if not do: return
|
||||
|
||||
if self.dim == 2:
|
||||
return self._refine2D()
|
||||
self._refine2D()
|
||||
elif self.dim == 3:
|
||||
return self._refine3D()
|
||||
self._refine3D()
|
||||
|
||||
# pass the refine function to the children
|
||||
if function is not None:
|
||||
@@ -590,8 +680,8 @@ class TreeMesh(object):
|
||||
return np.max([node.branchdepth for node in self.children.flatten('F')])
|
||||
|
||||
def refine(self, function):
|
||||
for node in self.children.flatten():
|
||||
node.refine(function)
|
||||
for cell in self.children.flatten():
|
||||
cell.refine(function)
|
||||
|
||||
def number(self):
|
||||
if self.isNumbered: return
|
||||
@@ -642,22 +732,31 @@ class TreeMesh(object):
|
||||
def nFy(self): return len(self.facesY)
|
||||
|
||||
@property
|
||||
def nFz(self): return len(self.facesZ)
|
||||
def nFz(self): return None if self.dim < 3 else len(self.facesZ)
|
||||
|
||||
@property
|
||||
def nE(self): return len(self.faces)
|
||||
def nE(self):
|
||||
if self.dim == 2:
|
||||
return len(self.faces)
|
||||
elif self.dim == 3:
|
||||
return len(self.edges)
|
||||
|
||||
@property
|
||||
def nEx(self):
|
||||
if self.dim == 2:
|
||||
return len(self.facesY)
|
||||
else: raise NotImplementedError('nEx')
|
||||
elif self.dim == 3:
|
||||
return len(self.edgesX)
|
||||
|
||||
@property
|
||||
def nEy(self):
|
||||
if self.dim == 2:
|
||||
return len(self.facesX)
|
||||
else: raise NotImplementedError('nEy')
|
||||
elif self.dim == 3:
|
||||
return len(self.edgesY)
|
||||
|
||||
@property
|
||||
def nEz(self): return None if self.dim < 3 else len(self.edgesZ)
|
||||
|
||||
@property
|
||||
def gridCC(self):
|
||||
@@ -745,12 +844,16 @@ class TreeMesh(object):
|
||||
self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S)
|
||||
return self._faceDiv
|
||||
|
||||
def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, showIt=False):
|
||||
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)
|
||||
|
||||
if plotC: [c.plotGrid(ax, text=text) for c in self.cells]
|
||||
if plotF: [f.plotGrid(ax, text=text) for f in self.faces]
|
||||
if plotE and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edges]
|
||||
if plotEx and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesX]
|
||||
if plotEy and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesY]
|
||||
if plotEz and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesZ]
|
||||
|
||||
ax.set_xlim((self.x0[0], self.h[0].sum()))
|
||||
ax.set_ylim((self.x0[1], self.h[1].sum()))
|
||||
@@ -797,7 +900,7 @@ if __name__ == '__main__':
|
||||
DIV = M.faceDiv
|
||||
# plt.subplot(211)
|
||||
# plt.spy(DIV)
|
||||
M.plotGrid(ax=plt.subplot(111),text=True)
|
||||
M.plotGrid(ax=plt.subplot(111),text=True,showIt=True)
|
||||
|
||||
q = np.zeros(M.nC)
|
||||
q[208] = -1.0
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
from SimPEG.Mesh.TreeMesh import TreeMesh, TreeFace, TreeCell
|
||||
import numpy as np
|
||||
import unittest
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
class TestOcTreeObjects(unittest.TestCase):
|
||||
|
||||
@@ -9,31 +10,156 @@ class TestOcTreeObjects(unittest.TestCase):
|
||||
self.Mr = TreeMesh([2,1,1])
|
||||
self.Mr.children[0,0,0].refine()
|
||||
self.Mr.number()
|
||||
# self.Mr.plotGrid(showIt=True,plotC=True,plotF=True)
|
||||
|
||||
def test_counts(self):
|
||||
ax = plt.subplot(111,projection='3d')
|
||||
# self.Mr.plotGrid(showIt=False,plotC=True,plotEy=True)
|
||||
|
||||
cell = self.Mr.sortedCells[1]
|
||||
[cell.edges[e].plotGrid(ax) for e in cell.edges]
|
||||
cell.plotGrid(ax)
|
||||
plt.show()
|
||||
|
||||
self.assertTrue(self.M.nC == 2)
|
||||
self.assertTrue(self.M.nFx == 3)
|
||||
self.assertTrue(self.M.nFy == 4)
|
||||
self.assertTrue(self.M.nFz == 4)
|
||||
self.assertTrue(self.M.nF == 11)
|
||||
self.assertTrue(self.M.nEx == 8)
|
||||
self.assertTrue(self.M.nEy == 6)
|
||||
self.assertTrue(self.M.nEz == 6)
|
||||
self.assertTrue(self.M.nE == 20)
|
||||
self.assertTrue(self.M.nN == 12)
|
||||
|
||||
print self.Mr.nN
|
||||
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 == 22)
|
||||
cell = self.Mr.sortedCells[1]
|
||||
self.assertTrue(cell.edges['eX0'].edgeType=='x')
|
||||
self.assertTrue(cell.edges['eX1'].edgeType=='x')
|
||||
self.assertTrue(cell.edges['eX2'].edgeType=='x')
|
||||
self.assertTrue(cell.edges['eX3'].edgeType=='x')
|
||||
self.assertTrue(cell.edges['eY0'].edgeType=='y')
|
||||
self.assertTrue(cell.edges['eY1'].edgeType=='y')
|
||||
self.assertTrue(cell.edges['eY2'].edgeType=='y')
|
||||
self.assertTrue(cell.edges['eY3'].edgeType=='y')
|
||||
self.assertTrue(cell.edges['eZ0'].edgeType=='z')
|
||||
self.assertTrue(cell.edges['eZ1'].edgeType=='z')
|
||||
self.assertTrue(cell.edges['eZ2'].edgeType=='z')
|
||||
self.assertTrue(cell.edges['eZ3'].edgeType=='z')
|
||||
|
||||
# def test_pointersM(self):
|
||||
# c0 = self.M.children[0,0,0]
|
||||
# c0fXm = c0.faces['fXm']
|
||||
# c0fXp = c0.faces['fXp']
|
||||
# c0fYm = c0.faces['fYm']
|
||||
# c0fYp = c0.faces['fYp']
|
||||
print self.Mr.nN
|
||||
# self.assertTrue(self.Mr.nN == 22)
|
||||
|
||||
def test_pointersM(self):
|
||||
c0 = self.M.children[0,0,0]
|
||||
c0fXm = c0.faces['fXm']
|
||||
c0fXp = c0.faces['fXp']
|
||||
c0fYm = c0.faces['fYm']
|
||||
c0fYp = c0.faces['fYp']
|
||||
|
||||
c1 = self.M.children[1,0,0]
|
||||
c1fXm = c1.faces['fXm']
|
||||
c1fXp = c1.faces['fXp']
|
||||
c1fYm = c1.faces['fYm']
|
||||
c1fYp = c1.faces['fYp']
|
||||
|
||||
self.assertTrue(c0fXp is c1fXm)
|
||||
self.assertTrue(c0fXp.edges['e0'] is c1fXm.edges['e0'])
|
||||
self.assertTrue(c0fXp.edges['e1'] is c1fXm.edges['e1'])
|
||||
self.assertTrue(c0fXp.edges['e2'] is c1fXm.edges['e2'])
|
||||
self.assertTrue(c0fXp.edges['e3'] is c1fXm.edges['e3'])
|
||||
self.assertTrue(c0fYp is not c1fYm)
|
||||
self.assertTrue(c0fXm is not c1fXm)
|
||||
|
||||
|
||||
def test_pointersMr(self):
|
||||
c0 = self.Mr.sortedCells[0]
|
||||
c0fXm = c0.faces['fXm']
|
||||
c0fXp = c0.faces['fXp']
|
||||
c0fYm = c0.faces['fYm']
|
||||
c0fYp = c0.faces['fYp']
|
||||
c0fZm = c0.faces['fZm']
|
||||
c0fZp = c0.faces['fZp']
|
||||
self.assertTrue(np.all(c0.center==np.r_[0.125,0.25,0.25]))
|
||||
|
||||
c1 = self.Mr.sortedCells[1]
|
||||
c1fXm = c1.faces['fXm']
|
||||
c1fXp = c1.faces['fXp']
|
||||
c1fYm = c1.faces['fYm']
|
||||
c1fYp = c1.faces['fYp']
|
||||
c1fZm = c1.faces['fZm']
|
||||
c1fZp = c1.faces['fZp']
|
||||
self.assertTrue(np.all(c1.center==np.r_[0.375,0.25,0.25]))
|
||||
|
||||
c2 = self.Mr.sortedCells[2]
|
||||
c2fXm = c2.faces['fXm']
|
||||
c2fXp = c2.faces['fXp']
|
||||
c2fYm = c2.faces['fYm']
|
||||
c2fYp = c2.faces['fYp']
|
||||
c2fZm = c2.faces['fZm']
|
||||
c2fZp = c2.faces['fZp']
|
||||
self.assertTrue(np.all(c2.center==np.r_[0.75,0.5,0.5]))
|
||||
|
||||
c4 = self.Mr.sortedCells[4]
|
||||
c4fXm = c4.faces['fXm']
|
||||
c4fXp = c4.faces['fXp']
|
||||
c4fYm = c4.faces['fYm']
|
||||
c4fYp = c4.faces['fYp']
|
||||
c4fZm = c4.faces['fZm']
|
||||
c4fZp = c4.faces['fZp']
|
||||
self.assertTrue(np.all(c4.center==np.r_[0.375,0.75,0.25]))
|
||||
|
||||
c6 = self.Mr.sortedCells[6]
|
||||
c6fXm = c6.faces['fXm']
|
||||
c6fXp = c6.faces['fXp']
|
||||
c6fYm = c6.faces['fYm']
|
||||
c6fYp = c6.faces['fYp']
|
||||
c6fZm = c6.faces['fZm']
|
||||
c6fZp = c6.faces['fZp']
|
||||
self.assertTrue(np.all(c6.center==np.r_[0.375,0.25,0.75]))
|
||||
|
||||
self.assertTrue(c0fXp is c1fXm)
|
||||
self.assertTrue(c0fYp is not c1fYm)
|
||||
self.assertTrue(c0fXm is not c1fXm)
|
||||
|
||||
self.assertTrue(c1fXp is c2fXm.children[0,0])
|
||||
self.assertTrue(c1fXp.parent is c2fXm)
|
||||
|
||||
self.assertTrue(c1fYp is c4fYm)
|
||||
self.assertTrue(c1fZp is c6fZm)
|
||||
|
||||
self.assertTrue(c6fXp is c2fXm.children[0,1])
|
||||
self.assertTrue(c6fXp.parent is c2fXm)
|
||||
|
||||
self.assertTrue(c4fXp is c2fXm.children[1,0])
|
||||
self.assertTrue(c4fXp.parent is c2fXm)
|
||||
|
||||
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)
|
||||
|
||||
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_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]
|
||||
# print self.Mr.gridFx - np.c_[x,y,z]
|
||||
# self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFx).flatten()) == 0)
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user