mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-13 17:45:30 +08:00
Edge Bug
This commit is contained in:
Rowan Cockett
+33
-30
@@ -227,7 +227,10 @@ class Tree(object):
|
||||
P += SortGrid(self.gridFx, offset=self.nEx)
|
||||
return sp.identity(self.nE).tocsr()[P,:]
|
||||
if self.dim == 3:
|
||||
raise Exception()
|
||||
P = SortGrid(self.gridEx)
|
||||
P += SortGrid(self.gridEy, offset=self.nEx)
|
||||
P += SortGrid(self.gridEz, offset=self.nEx+self.nEy)
|
||||
return sp.identity(self.nE).tocsr()[P,:]
|
||||
|
||||
def _index(self, pointer):
|
||||
assert len(pointer) is self.dim+1
|
||||
@@ -705,7 +708,7 @@ class Tree(object):
|
||||
gridEy.append( [n[0], n[1] + h[1]/2.0, n[2]] )
|
||||
edgeEy.append( h[1] )
|
||||
self._gridEy = np.array(gridEy)
|
||||
self._edgeEyFull = np.array(edgeEx)
|
||||
self._edgeEyFull = np.array(edgeEy)
|
||||
|
||||
gridEz = []
|
||||
edgeEz = []
|
||||
@@ -717,7 +720,7 @@ class Tree(object):
|
||||
gridEz.append( [n[0], n[1], n[2] + h[2]/2.0] )
|
||||
edgeEz.append( h[2] )
|
||||
self._gridEz = np.array(gridEz)
|
||||
self._edgeEzFull = np.array(edgeEx)
|
||||
self._edgeEzFull = np.array(edgeEz)
|
||||
|
||||
self.__dirtyEdges__ = False
|
||||
|
||||
@@ -1123,14 +1126,14 @@ class Tree(object):
|
||||
ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r:')
|
||||
ax.plot(self.gridCC[[0,-1],0], self.gridCC[[0,-1],1], 'ro')
|
||||
if nodes:
|
||||
ax.plot(self.gridN[:,0], self.gridN[:,1], 'ms')
|
||||
ax.plot(self.gridN[self._hangingN.keys(),0], self.gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m')
|
||||
ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms')
|
||||
ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m')
|
||||
if facesX:
|
||||
ax.plot(self.gridFx[self._hangingFx.keys(),0], self.gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.plot(self.gridFx[:,0], self.gridFx[:,1], 'g>')
|
||||
ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>')
|
||||
if facesY:
|
||||
ax.plot(self.gridFy[self._hangingFy.keys(),0], self.gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.plot(self.gridFy[:,0], self.gridFy[:,1], 'g^')
|
||||
ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
|
||||
ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^')
|
||||
elif self.dim == 3:
|
||||
if cells:
|
||||
ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.', zs=self.gridCC[:,2])
|
||||
@@ -1139,61 +1142,61 @@ class Tree(object):
|
||||
ax.plot(self.gridCC[[0,-1],0], self.gridCC[[0,-1],1], 'ro', zs=self.gridCC[[0,-1],2])
|
||||
|
||||
if nodes:
|
||||
ax.plot(self.gridN[:,0], self.gridN[:,1], 'ms', zs=self.gridN[:,2])
|
||||
ax.plot(self.gridN[self._hangingN.keys(),0], self.gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m', zs=self.gridN[self._hangingN.keys(),2])
|
||||
ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms', zs=self._gridN[:,2])
|
||||
ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m', zs=self._gridN[self._hangingN.keys(),2])
|
||||
for key in self._hangingN.keys():
|
||||
for hf in self._hangingN[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridN[ind,0], self.gridN[ind,1], 'm:', zs=self.gridN[ind,2])
|
||||
ax.plot(self._gridN[ind,0], self._gridN[ind,1], 'm:', zs=self._gridN[ind,2])
|
||||
|
||||
if facesX:
|
||||
ax.plot(self.gridFx[:,0], self.gridFx[:,1], 'g>', zs=self.gridFx[:,2])
|
||||
ax.plot(self.gridFx[self._hangingFx.keys(),0], self.gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g', zs=self.gridFx[self._hangingFx.keys(),2])
|
||||
ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>', zs=self._gridFx[:,2])
|
||||
ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g', zs=self._gridFx[self._hangingFx.keys(),2])
|
||||
for key in self._hangingFx.keys():
|
||||
for hf in self._hangingFx[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridFx[ind,0], self.gridFx[ind,1], 'g:', zs=self.gridFx[ind,2])
|
||||
ax.plot(self._gridFx[ind,0], self._gridFx[ind,1], 'g:', zs=self._gridFx[ind,2])
|
||||
|
||||
if facesY:
|
||||
ax.plot(self.gridFy[:,0], self.gridFy[:,1], 'g^', zs=self.gridFy[:,2])
|
||||
ax.plot(self.gridFy[self._hangingFy.keys(),0], self.gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g', zs=self.gridFy[self._hangingFy.keys(),2])
|
||||
ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^', zs=self._gridFy[:,2])
|
||||
ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g', zs=self._gridFy[self._hangingFy.keys(),2])
|
||||
for key in self._hangingFy.keys():
|
||||
for hf in self._hangingFy[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridFy[ind,0], self.gridFy[ind,1], 'g:', zs=self.gridFy[ind,2])
|
||||
ax.plot(self._gridFy[ind,0], self._gridFy[ind,1], 'g:', zs=self._gridFy[ind,2])
|
||||
|
||||
if facesZ:
|
||||
ax.plot(self.gridFz[:,0], self.gridFz[:,1], 'g^', zs=self.gridFz[:,2])
|
||||
ax.plot(self.gridFz[self._hangingFz.keys(),0], self.gridFz[self._hangingFz.keys(),1], 'gs', ms=10, mfc='none', mec='g', zs=self.gridFz[self._hangingFz.keys(),2])
|
||||
ax.plot(self._gridFz[:,0], self._gridFz[:,1], 'g^', zs=self._gridFz[:,2])
|
||||
ax.plot(self._gridFz[self._hangingFz.keys(),0], self._gridFz[self._hangingFz.keys(),1], 'gs', ms=10, mfc='none', mec='g', zs=self._gridFz[self._hangingFz.keys(),2])
|
||||
for key in self._hangingFz.keys():
|
||||
for hf in self._hangingFz[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridFz[ind,0], self.gridFz[ind,1], 'g:', zs=self.gridFz[ind,2])
|
||||
ax.plot(self._gridFz[ind,0], self._gridFz[ind,1], 'g:', zs=self._gridFz[ind,2])
|
||||
|
||||
if edgesX:
|
||||
ax.plot(self.gridEx[:,0], self.gridEx[:,1], 'k>', zs=self.gridEx[:,2])
|
||||
ax.plot(self.gridEx[self._hangingEx.keys(),0], self.gridEx[self._hangingEx.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self.gridEx[self._hangingEx.keys(),2])
|
||||
ax.plot(self._gridEx[:,0], self._gridEx[:,1], 'k>', zs=self._gridEx[:,2])
|
||||
ax.plot(self._gridEx[self._hangingEx.keys(),0], self._gridEx[self._hangingEx.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self._gridEx[self._hangingEx.keys(),2])
|
||||
for key in self._hangingEx.keys():
|
||||
for hf in self._hangingEx[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridEx[ind,0], self.gridEx[ind,1], 'k:', zs=self.gridEx[ind,2])
|
||||
ax.plot(self._gridEx[ind,0], self._gridEx[ind,1], 'k:', zs=self._gridEx[ind,2])
|
||||
|
||||
|
||||
if edgesY:
|
||||
ax.plot(self.gridEy[:,0], self.gridEy[:,1], 'k<', zs=self.gridEy[:,2])
|
||||
ax.plot(self.gridEy[self._hangingEy.keys(),0], self.gridEy[self._hangingEy.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self.gridEy[self._hangingEy.keys(),2])
|
||||
ax.plot(self._gridEy[:,0], self._gridEy[:,1], 'k<', zs=self._gridEy[:,2])
|
||||
ax.plot(self._gridEy[self._hangingEy.keys(),0], self._gridEy[self._hangingEy.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self._gridEy[self._hangingEy.keys(),2])
|
||||
for key in self._hangingEy.keys():
|
||||
for hf in self._hangingEy[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridEy[ind,0], self.gridEy[ind,1], 'k:', zs=self.gridEy[ind,2])
|
||||
ax.plot(self._gridEy[ind,0], self._gridEy[ind,1], 'k:', zs=self._gridEy[ind,2])
|
||||
|
||||
if edgesZ:
|
||||
ax.plot(self.gridEz[:,0], self.gridEz[:,1], 'k^', zs=self.gridEz[:,2])
|
||||
ax.plot(self.gridEz[self._hangingEz.keys(),0], self.gridEz[self._hangingEz.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self.gridEz[self._hangingEz.keys(),2])
|
||||
ax.plot(self._gridEz[:,0], self._gridEz[:,1], 'k^', zs=self._gridEz[:,2])
|
||||
ax.plot(self._gridEz[self._hangingEz.keys(),0], self._gridEz[self._hangingEz.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self._gridEz[self._hangingEz.keys(),2])
|
||||
for key in self._hangingEz.keys():
|
||||
for hf in self._hangingEz[key]:
|
||||
ind = [key, hf[0]]
|
||||
ax.plot(self.gridEz[ind,0], self.gridEz[ind,1], 'k:', zs=self.gridEz[ind,2])
|
||||
ax.plot(self._gridEz[ind,0], self._gridEz[ind,1], 'k:', zs=self._gridEz[ind,2])
|
||||
|
||||
|
||||
# ax.axis('equal')
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -51,5 +51,5 @@ setup(
|
||||
platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],
|
||||
use_2to3 = False,
|
||||
include_dirs=[np.get_include()],
|
||||
ext_modules = cythonize('SimPEG/Utils/interputils_cython.pyx')
|
||||
ext_modules = cythonize(['SimPEG/Utils/interputils_cython.pyx', 'SimPEG/Mesh/TreeUtils.pyx'])
|
||||
)
|
||||
|
||||
@@ -12,44 +12,41 @@ class TestSimpleQuadTree(unittest.TestCase):
|
||||
|
||||
def test_counts(self):
|
||||
|
||||
T = Tree([8,8])
|
||||
T._refineCell([0,0,0])
|
||||
T._refineCell([4,4,1])
|
||||
T._refineCell([0,0,1])
|
||||
T._refineCell([2,2,2])
|
||||
|
||||
T.number()
|
||||
T.plotGrid(showIt=True)
|
||||
assert sorted(T._treeInds) == [2, 34, 66, 99, 107, 115, 123, 129, 257, 386, 418, 450, 482]
|
||||
assert len(T._hangingFacesX) == 7
|
||||
assert T.nFx == 18
|
||||
assert T.vol.sum() == 1.0
|
||||
M = Tree([8,8])
|
||||
M._refineCell([0,0,0])
|
||||
M._refineCell([0,0,1])
|
||||
M.number()
|
||||
# M.plotGrid(showIt=True)
|
||||
# assert sorted(M._cells) == [2, 34, 66, 99, 107, 115, 123, 129, 257, 386, 418, 450, 482]
|
||||
assert M.nhFx == 2
|
||||
assert M.nFx == 9
|
||||
assert M.vol.sum() == 1.0
|
||||
|
||||
|
||||
def test_connectivity(self):
|
||||
T = Tree([8,8])
|
||||
T._refineCell([0,0,0])
|
||||
T._refineCell([4,4,1])
|
||||
T._refineCell([0,0,1])
|
||||
T._refineCell([2,2,2])
|
||||
T.number()
|
||||
assert T._getNextCell([4,0,1]) is None
|
||||
assert T._getNextCell([0,4,1]) == [T._index([4,4,2]), T._index([4,6,2])]
|
||||
assert T._getNextCell([0,2,2]) == [T._index([2,2,3]), T._index([2,3,3])]
|
||||
assert T._getNextCell([4,4,2]) == T._index([6,4,2])
|
||||
assert T._getNextCell([6,4,2]) is None
|
||||
assert T._getNextCell([2,0,2]) == T._index([4,0,1])
|
||||
assert T._getNextCell([4,0,1], positive=False) == [T._index([2,0,2]), [T._index([3,2,3]), T._index([3,3,3])]]
|
||||
assert T._getNextCell([3,3,3]) == T._index([4,0,1])
|
||||
assert T._getNextCell([3,2,3]) == T._index([4,0,1])
|
||||
assert T._getNextCell([2,2,3]) == T._index([3,2,3])
|
||||
assert T._getNextCell([3,2,3], positive=False) == T._index([2,2,3])
|
||||
# def test_connectivity(self):
|
||||
# T = Tree([8,8])
|
||||
# T._refineCell([0,0,0])
|
||||
# T._refineCell([4,4,1])
|
||||
# T._refineCell([0,0,1])
|
||||
# T._refineCell([2,2,2])
|
||||
# T.number()
|
||||
# assert T._getNextCell([4,0,1]) is None
|
||||
# assert T._getNextCell([0,4,1]) == [T._index([4,4,2]), T._index([4,6,2])]
|
||||
# assert T._getNextCell([0,2,2]) == [T._index([2,2,3]), T._index([2,3,3])]
|
||||
# assert T._getNextCell([4,4,2]) == T._index([6,4,2])
|
||||
# assert T._getNextCell([6,4,2]) is None
|
||||
# assert T._getNextCell([2,0,2]) == T._index([4,0,1])
|
||||
# assert T._getNextCell([4,0,1], positive=False) == [T._index([2,0,2]), [T._index([3,2,3]), T._index([3,3,3])]]
|
||||
# assert T._getNextCell([3,3,3]) == T._index([4,0,1])
|
||||
# assert T._getNextCell([3,2,3]) == T._index([4,0,1])
|
||||
# assert T._getNextCell([2,2,3]) == T._index([3,2,3])
|
||||
# assert T._getNextCell([3,2,3], positive=False) == T._index([2,2,3])
|
||||
|
||||
|
||||
assert T._getNextCell([0,0,2], direction=1) == T._index([0,2,2])
|
||||
assert T._getNextCell([0,2,2], direction=1, positive=False) == T._index([0,0,2])
|
||||
assert T._getNextCell([0,2,2], direction=1) == T._index([0,4,1])
|
||||
assert T._getNextCell([0,4,1], direction=1, positive=False) == [T._index([0,2,2]), [T._index([2,3,3]), T._index([3,3,3])]]
|
||||
# assert T._getNextCell([0,0,2], direction=1) == T._index([0,2,2])
|
||||
# assert T._getNextCell([0,2,2], direction=1, positive=False) == T._index([0,0,2])
|
||||
# assert T._getNextCell([0,2,2], direction=1) == T._index([0,4,1])
|
||||
# assert T._getNextCell([0,4,1], direction=1, positive=False) == [T._index([0,2,2]), [T._index([2,3,3]), T._index([3,3,3])]]
|
||||
|
||||
|
||||
class TestOperatorsQuadTree(unittest.TestCase):
|
||||
@@ -81,29 +78,29 @@ class TestOperatorsQuadTree(unittest.TestCase):
|
||||
|
||||
class TestOperatorsOcTree(unittest.TestCase):
|
||||
|
||||
def test_counts(self):
|
||||
def test_faceDiv(self):
|
||||
|
||||
hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
|
||||
T = Tree([hx, hy, hz], levels=2)
|
||||
T.refine(lambda xc:2)
|
||||
# T.plotGrid(showIt=True)
|
||||
M = Mesh.TensorMesh([hx, hy, hz])
|
||||
assert M.nC == T.nC
|
||||
assert M.nF == T.nF
|
||||
assert M.nFx == T.nFx
|
||||
assert M.nFy == T.nFy
|
||||
# assert M.nE == T.nE
|
||||
# assert M.nEx == T.nEx
|
||||
# assert M.nEy == T.nEy
|
||||
assert np.allclose(M.area, T.permuteF*T.area)
|
||||
# assert np.allclose(M.edge, T.permuteE*T.edge)
|
||||
assert np.allclose(M.vol, T.permuteCC*T.vol)
|
||||
M = Tree([hx, hy, hz], levels=2)
|
||||
M.refine(lambda xc:2)
|
||||
# M.plotGrid(showIt=True)
|
||||
Mr = Mesh.TensorMesh([hx, hy, hz])
|
||||
assert M.nC == Mr.nC
|
||||
assert M.nF == Mr.nF
|
||||
assert M.nFx == Mr.nFx
|
||||
assert M.nFy == Mr.nFy
|
||||
assert M.nE == Mr.nE
|
||||
assert M.nEx == Mr.nEx
|
||||
assert M.nEy == Mr.nEy
|
||||
assert np.allclose(Mr.area, M.permuteF*M.area)
|
||||
assert np.allclose(Mr.edge, M.permuteE*M.edge)
|
||||
assert np.allclose(Mr.vol, M.permuteCC*M.vol)
|
||||
|
||||
# plt.subplot(211).spy(M.faceDiv)
|
||||
# plt.subplot(212).spy(T.permuteCC.T*T.faceDiv*T.permuteF)
|
||||
# plt.subplot(211).spy(Mr.faceDiv)
|
||||
# plt.subplot(212).spy(M.permuteCC.T*M.faceDiv*M.permuteF)
|
||||
# plt.show()
|
||||
|
||||
assert (M.faceDiv - T.permuteCC*T.faceDiv*T.permuteF.T).nnz == 0
|
||||
assert (Mr.faceDiv - M.permuteCC*M.faceDiv*M.permuteF.T).nnz == 0
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user