diff --git a/SimPEG/Mesh/PointerTree.py b/SimPEG/Mesh/PointerTree.py index 339e3008..1db06c1a 100644 --- a/SimPEG/Mesh/PointerTree.py +++ b/SimPEG/Mesh/PointerTree.py @@ -1,51 +1,10 @@ from SimPEG import np, sp, Utils, Solver import matplotlib.pyplot as plt import matplotlib +import TreeUtils +import time -class ZCurve(object): - """ - The Z-order curve is generated by interleaving the bits of an offset. - - See: - - https://github.com/cortesi/scurve - Aldo Cortesi - - """ - def __init__(self, dimension, bits): - """ - dimension: Number of dimensions - bits: The number of bits per co-ordinate. Total number of points is - 2**(bits*dimension). - """ - self.dimension, self.bits = dimension, bits - - def bitrange(self, x, width, start, end): - """ - Extract a bit range as an integer. - (start, end) is inclusive lower bound, exclusive upper bound. - """ - return x >> (width-end) & ((2**(end-start))-1) - - def index(self, p): - p.reverse() - idx = 0 - iwidth = self.bits * self.dimension - for i in range(iwidth): - bitoff = self.bits-(i/self.dimension)-1 - poff = self.dimension-(i%self.dimension)-1 - b = self.bitrange(p[poff], self.bits, bitoff, bitoff+1) << i - idx |= b - return idx - - def point(self, idx): - p = [0]*self.dimension - iwidth = self.bits * self.dimension - for i in range(iwidth): - b = self.bitrange(idx, iwidth, i, i+1) << (iwidth-i-1)/self.dimension - p[i%self.dimension] |= b - p.reverse() - return p +MAX_BITS = 20 def SortGrid(grid, offset=0): """ @@ -112,7 +71,6 @@ class Tree(object): self.__dirty__ = True #: The numbering is dirty! - self._z = ZCurve(self.dim, 20) self._treeInds = set() self._treeInds.add(0) @@ -216,8 +174,6 @@ class Tree(object): raise Exception() def _structureChange(self): - if self.__dirty__: return - deleteThese = ['__sortedInds', '_gridCC', '_gridN', '_gridFx', '_gridFy', '_gridFz', '_gridEx', '_gridEy', '_gridEz', '_area', '_edge', '_vol'] for p in deleteThese: if hasattr(self, p): delattr(self, p) @@ -226,39 +182,42 @@ class Tree(object): def _index(self, pointer): assert len(pointer) is self.dim+1 assert pointer[-1] <= self.levels - x = self._z.index([p for p in pointer[:-1]]) # copy - return (x << self._levelBits) + pointer[-1] + return TreeUtils.index(self.dim, MAX_BITS, self._levelBits, pointer[:-1], pointer[-1]) def _pointer(self, index): assert type(index) in [int, long] - n = index & (2**self._levelBits-1) - p = self._z.point(index >> self._levelBits) - return p + [n] + return TreeUtils.point(self.dim, MAX_BITS, self._levelBits, index) def __contains__(self, v): if type(v) in [int, long]: return v in self._treeInds return self._index(v) in self._treeInds - def refine(self, function=None, recursive=True, cells=None, balance=True): + def refine(self, function=None, recursive=True, cells=None, balance=True, _inRecursion=False): + + if not _inRecursion: + self._structureChange() + print 'Refining Mesh' cells = cells if cells is not None else sorted(self._treeInds) recurse = [] + tic = time.time() for cell in cells: p = self._pointer(cell) do = function(self._cellC(cell)) > p[-1] if do: recurse += self._refineCell(cell) - if recursive and len(recurse) > 0: - recurse += self.refine(function=function, recursive=True, cells=recurse, balance=balance) + print ' ', time.time() - tic - if balance: - recurse += self.balance(recursive=True) + if recursive and len(recurse) > 0: + recurse += self.refine(function=function, recursive=True, cells=recurse, balance=balance, _inRecursion=True) + + if balance and not _inRecursion: + self.balance() return recurse def _refineCell(self, pointer): - self._structureChange() pointer = self._asPointer(pointer) ind = self._asIndex(pointer) assert ind in self @@ -283,10 +242,13 @@ class Tree(object): self._treeInds.remove(ind) return added - def _corsenCell(self, pointer): + def corsen(self, function=None): self._structureChange() raise Exception('Not yet implemented') + def _corsenCell(self, pointer): + raise Exception('Not yet implemented') + def _asPointer(self, ind): if type(ind) in [int, long]: return self._pointer(ind) @@ -393,10 +355,24 @@ class Tree(object): return self._getNextCell(self._parentPointer(pointer), direction=direction, positive=positive) - def balance(self, recursive=True): + def balance(self, recursive=True, cells=None, _inRecursion=False): + + tic = time.time() + if not _inRecursion: + self._structureChange() + print 'Balancing Mesh:' + + cells = cells if cells is not None else sorted(self._treeInds) + + # calcDepth = lambda i: lambda A: i if type(A) is not list else max(map(calcDepth(i+1), A)) + # flatten = lambda A: A if calcDepth(0)(A) == 1 else flatten([_ for __ in A for _ in (__ if type(__) is list else [__])]) + + recurse = set() + + for cell in cells: + p = self._asPointer(cell) + if p[-1] == self.levels: continue - recurse = [] - for cell in sorted(self._treeInds): cs = range(6) cs[0] = self._getNextCell(cell, direction=0, positive=False) cs[1] = self._getNextCell(cell, direction=0, positive=True) @@ -410,14 +386,20 @@ class Tree(object): for c in cs if c is not None ]) - # print cs, do - if do: - recurse += self._refineCell(cell) + # depth = calcDepth(0)(cs) + # print depth, depth > 2, do, [jj for jj in flatten(cs) if jj is not None] + # recurse += [jj for jj in flatten(cs) if jj is not None] + if do and cell in self: + newCells = self._refineCell(cell) + recurse.update([_ for _ in cs if type(_) in [int, long]]) # only add the bigger ones! + recurse.update(newCells) + + print ' ', len(cells), time.time() - tic if recursive and len(recurse) > 0: - # print 'here' - recurse += self.balance() - return recurse + self.balance(cells=sorted(recurse), _inRecursion=True) + + @property def gridCC(self): @@ -961,7 +943,7 @@ class Tree(object): facesX=False, facesY=False, facesZ=False, edgesX=False, edgesY=False, edgesZ=False): - self.number() + # self.number() axOpts = {'projection':'3d'} if self.dim == 3 else {} if ax is None: @@ -1070,7 +1052,7 @@ class Tree(object): ax.plot(self.gridEz[ind,0], self.gridEz[ind,1], 'k:', zs=self.gridEz[ind,2]) - ax.axis('equal') + # ax.axis('equal') if showIt:plt.show() @@ -1079,36 +1061,47 @@ if __name__ == '__main__': def function(xc): - r = xc - np.r_[0.5*128,0.5*128] + r = xc - np.array([0.5*128]*len(xc)) dist = np.sqrt(r.dot(r)) # if dist < 0.05: # return 5 if dist < 0.1*128: - return 5 + return 7 if dist < 0.3*128: - return 3 - # if dist < 1.0*128: - # return 2 + return 5 + if dist < 1.0*128: + return 2 else: return 0 - # T = Tree([[(1,8)],[(1,8)],[(1,8)]],levels=3) + T = Tree([[(1,128)],[(1,128)],[(1,128)]],levels=7) # T = Tree([[(1,16)],[(1,16)]],levels=4) - T = Tree([[(1,128)],[(1,128)]],levels=7) - T.refine(lambda xc:2, balance=False) - T.refine(function) + # T = Tree([[(1,128)],[(1,128)]],levels=7) + # T.refine(lambda xc:6, balance=False) + # T._index([0,0,0]) + # T._pointer(0) + + + tic = time.time() + T.refine(function)#, balance=False) + print time.time() - tic + print T.nC + + asdf # T._refineCell([8,0,1]) # T._refineCell([8,0,2]) # T._refineCell([12,0,2]) # T._refineCell([8,4,2]) # T._refineCell([6,0,3]) # T._refineCell([8,8,1]) - # T._refineCell([4,4,4,1]) + # T._refineCell([0,0,0,1]) + ax = plt.subplot(211) - ax.spy(T.faceDiv) + # ax.spy(T.faceDiv) + T.plotGrid(ax=ax) # R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i) # print R @@ -1120,7 +1113,7 @@ if __name__ == '__main__': # ax.spy(T.faceDiv[:,:T.nFx] * R) - + T.balance() T.plotGrid(ax=ax) # cx = T._getNextCell([0,0,1],direction=0,positive=True) @@ -1145,7 +1138,6 @@ if __name__ == '__main__': - # print T.nN plt.show() diff --git a/SimPEG/Mesh/TreeUtils.pyx b/SimPEG/Mesh/TreeUtils.pyx new file mode 100644 index 00000000..d9bf5238 --- /dev/null +++ b/SimPEG/Mesh/TreeUtils.pyx @@ -0,0 +1,85 @@ +# from __future__ import division +# import numpy as np +# cimport numpy as np +# from libcpp.vector cimport vector + + +""" + The Z-order curve is generated by interleaving the bits of an offset. + + See also: + + https://github.com/cortesi/scurve + Aldo Cortesi + +""" + +def bitrange(long x, int width, int start, int end): + """ + Extract a bit range as an integer. + (start, end) is inclusive lower bound, exclusive upper bound. + """ + return x >> (width-end) & ((2**(end-start))-1) + +def index(int dimension, int bits, int levelBits, list p, int level): + cdef long idx = 0 + cdef int iwidth + cdef int i + cdef long b + cdef int bitoff + + p = [_ for _ in p] + + p.reverse() + iwidth = bits * dimension + for i in range(iwidth): + bitoff = bits-(i/dimension)-1 + poff = dimension-(i%dimension)-1 + b = bitrange(p[poff], bits, bitoff, bitoff+1) << i + idx |= b + + return (idx << levelBits) + level + +def point(int dimension, int bits, int levelBits, long idx): + cdef list p + cdef int iwidth + cdef int i, n + cdef long b + + n = idx & (2**levelBits-1) + idx = idx >> levelBits + + p = [0]*dimension + iwidth = bits * dimension + for i in range(iwidth): + b = bitrange(idx, iwidth, i, i+1) << (iwidth-i-1)/dimension + p[i%dimension] |= b + p.reverse() + return p + [n] + + +# def _refineCell(int dimension, int bits, self, pointer): +# self._structureChange() +# pointer = self._asPointer(pointer) +# ind = self._asIndex(pointer) +# assert ind in self +# h = self._levelWidth(pointer[-1])/2 # halfWidth +# nL = pointer[-1] + 1 # new level +# add = lambda p:p[0]+p[1] +# added = [] +# def addCell(p): +# i = self._index(p+[nL]) +# self._treeInds.add(i) +# added.append(i) + +# addCell(map(add, zip(pointer[:-1], [0,0,0]))) +# addCell(map(add, zip(pointer[:-1], [h,0,0]))) +# addCell(map(add, zip(pointer[:-1], [0,h,0]))) +# addCell(map(add, zip(pointer[:-1], [h,h,0]))) +# if self.dim == 3: +# addCell(map(add, zip(pointer[:-1], [0,0,h]))) +# addCell(map(add, zip(pointer[:-1], [h,0,h]))) +# addCell(map(add, zip(pointer[:-1], [0,h,h]))) +# addCell(map(add, zip(pointer[:-1], [h,h,h]))) +# self._treeInds.remove(ind) +# return added