from SimPEG import np, sp, Utils, Solver import matplotlib.pyplot as plt import matplotlib 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 def SortGrid(grid, offset=0): """ Sorts a grid by the x0 location. """ eps = 1e-7 def mycmp(c1,c2): c1 = grid[c1-offset] c2 = grid[c2-offset] if c1.size == 2: if np.abs(c1[1] - c2[1]) < eps: return c1[0] - c2[0] return c1[1] - c2[1] elif c1.size == 3: if np.abs(c1[2] - c2[2]) < eps: if np.abs(c1[1] - c2[1]) < eps: return c1[0] - c2[0] return c1[1] - c2[1] return c1[2] - c2[2] 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 sorted(range(offset,grid.shape[0]+offset), key=K) class NotBalancedException(Exception): pass class Tree(object): def __init__(self, h_in, levels=3): assert type(h_in) is list, 'h_in must be a list' assert len(h_in) > 1, "len(h_in) must be greater than 1" h = range(len(h_in)) for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: # This gives you something over the unit cube. h_i = np.ones(int(h_i))/int(h_i) elif type(h_i) is list: h_i = Utils.meshTensor(h_i) assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i) assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) assert len(h_i) == 2**levels, "must make h and levels match" h[i] = h_i[:] # make a copy. self.h = h self._levels = levels self._levelBits = int(np.ceil(np.sqrt(levels)))+1 self.__dirty__ = True #: The numbering is dirty! self._z = ZCurve(self.dim, 20) self._treeInds = set() self._treeInds.add(0) @property def __dirty__(self): return self.__dirtyFaces__ or self.__dirtyEdges__ or self.__dirtyNodes__ or self.__dirtyHanging__ @__dirty__.setter def __dirty__(self, val): assert val is True self.__dirtyFaces__ = True self.__dirtyEdges__ = True self.__dirtyNodes__ = True self.__dirtyHanging__ = True @property def levels(self): return self._levels @property def dim(self): return len(self.h) @property def nC(self): return len(self._treeInds) @property def nN(self): self._numberNodes() return len(self._nodes) @property def nF(self): return self.nFx + self.nFy + (0 if self.dim == 2 else self.nFz) @property def nFx(self): self._numberFaces() return len(self._facesX) @property def nFy(self): self._numberFaces() return len(self._facesY) @property def nFz(self): if self.dim == 2: return None self._numberFaces() return len(self._facesZ) @property def nE(self): return self.nEx + self.nEy + (0 if self.dim == 2 else self.nEz) @property def nEx(self): if self.dim == 2:return self.nFy self._numberEdges() return len(self._edgesX) @property def nEy(self): if self.dim == 2:return self.nFx self._numberEdges() return len(self._edgesY) @property def nEz(self): if self.dim == 2: return None self._numberEdges() return len(self._edgesZ) @property def vol(self): self.number() return self._vol @property def area(self): self.number() return self._area @property def edge(self): self.number() if self.dim == 2: return np.r_[self._area[self.nFx:], self._area[:self.nFx]] @property def _sortedInds(self): if getattr(self, '__sortedInds', None) is None: self.__sortedInds = sorted(self._treeInds) return self.__sortedInds @property def permuteCC(self): #TODO: cache these? P = SortGrid(self.gridCC) return sp.identity(self.nC).tocsr()[P,:] @property def permuteF(self): #TODO: cache these? P = SortGrid(self.gridFx) P += SortGrid(self.gridFy, offset=self.nFx) if self.dim == 3: P += SortGrid(self.gridFz, offset=self.nFx+self.nFy) return sp.identity(self.nF).tocsr()[P,:] @property def permuteE(self): #TODO: cache these? if self.dim == 2: P = SortGrid(self.gridFy) P += SortGrid(self.gridFx, offset=self.nEx) return sp.identity(self.nE).tocsr()[P,:] if self.dim == 3: raise Exception() def _structureChange(self): if self.__dirty__: return deleteThese = ['__sortedInds', '_gridCC', '_gridFx'] for p in deleteThese: if hasattr(self, p): delattr(self, p) self.__dirty__ = True 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] 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] 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): cells = cells if cells is not None else sorted(self._treeInds) recurse = [] 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: self.refine(function=function, recursive=True, cells=recurse) return recurse def _refineCell(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][:self.dim]))) addCell(map(add, zip(pointer[:-1], [h,0,0][:self.dim]))) addCell(map(add, zip(pointer[:-1], [0,h,0][:self.dim]))) addCell(map(add, zip(pointer[:-1], [h,h,0][:self.dim]))) 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 def _corsenCell(self, pointer): self._structureChange() raise Exception('Not yet implemented') def _asPointer(self, ind): if type(ind) in [int, long]: return self._pointer(ind) if type(ind) is list: return ind if isinstance(ind, np.ndarray): return ind.tolist() raise Exception def _asIndex(self, pointer): if type(pointer) in [int, long]: return pointer if type(pointer) is list: return self._index(pointer) raise Exception def _parentPointer(self, pointer): mod = self._levelWidth(pointer[-1]-1) return [p - (p % mod) for p in pointer[:-1]] + [pointer[-1]-1] def _cellN(self, p): p = self._asPointer(p) return [hi[:p[ii]].sum() for ii, hi in enumerate(self.h)] def _cellH(self, p): p = self._asPointer(p) w = self._levelWidth(p[-1]) return [hi[p[ii]:p[ii]+w].sum() for ii, hi in enumerate(self.h)] def _cellC(self, p): return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist() def _levelWidth(self, level): return 2**(self.levels - level) def _isInsideMesh(self, pointer): inside = True for p in pointer[:-1]: inside = inside and p >= 0 and p < 2**self.levels return inside def _getNextCell(self, ind, direction=0, positive=True): """ Returns a None, int, list, or nested list The int is the cell number. """ pointer = self._asPointer(ind) step = (1 if positive else -1) * self._levelWidth(pointer[-1]) nextCell = [p if ii is not direction else p + step for ii, p in enumerate(pointer)] if not self._isInsideMesh(nextCell): return None # it might be the same size as me? if nextCell in self: return self._index(nextCell) # it might be smaller than me? if nextCell[-1] + 1 <= self.levels: # if I am not the smallest. nextCell[-1] += 1 if not positive: nextCell[direction] -= step/2 # Get the closer one if nextCell in self: # there is at least one hw = self._levelWidth(pointer[-1]) / 2 nextCell = np.array([p if ii is not direction else p + (step/2 if positive else 0) for ii, p in enumerate(pointer)]) if self.dim == 2: if direction == 0: children = [0,0,1], [0,hw,1] if direction == 1: children = [0,0,1], [hw,0,1] elif self.dim == 3: if direction == 0: children = [0,0,0,1], [0,hw,0,1], [0,0,hw,1], [0,hw,hw,1] if direction == 1: children = [0,0,0,1], [hw,0,0,1], [0,0,hw,1], [hw,0,hw,1] if direction == 2: children = [0,0,0,1], [hw,0,0,1], [0,hw,0,1], [hw,hw,0,1] nextCells = [] for child in children: nextCells.append(self._getNextCell(nextCell + child, direction=direction,positive=positive)) return nextCells # it might be bigger than me? return self._getNextCell(self._parentPointer(pointer), direction=direction, positive=positive) @property def gridCC(self): if getattr(self, '_gridCC', None) is None: self._gridCC = np.zeros((len(self._treeInds),self.dim)) for ii, ind in enumerate(self._sortedInds): p = self._asPointer(ind) self._gridCC[ii, :] = self._cellC(p) return self._gridCC @property def gridN(self): self._numberNodes() return self._gridN @property def gridFx(self): self._numberFaces() return self._gridFx @property def gridFy(self): self._numberFaces() return self._gridFy @property def gridFz(self): if self.dim < 3: return None self._numberFaces() return self._gridFz @property def gridEx(self): if self.dim == 2: return self.gridFy self._numberEdges() return self._gridEx @property def gridEy(self): if self.dim == 2: return self.gridFx self._numberEdges() return self._gridEy @property def gridEz(self): if self.dim < 3: return None self._numberEdges() return self._gridEz def _onSameLevel(self, i0, i1): p0 = self._asPointer(i0) p1 = self._asPointer(i1) return p0[-1] == p1[-1] def _numberNodes(self, force=False): if not self.__dirtyNodes__ and not force: return self._nodes = set() for ind in self._treeInds: p = self._asPointer(ind) w = self._levelWidth(p[-1]) if self.dim == 2: self._nodes.add(self._index([p[0] , p[1] , p[2]])) self._nodes.add(self._index([p[0] + w, p[1] , p[2]])) self._nodes.add(self._index([p[0] , p[1] + w, p[2]])) self._nodes.add(self._index([p[0] + w, p[1] + w, p[2]])) elif self.dim == 3: self._nodes.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._nodes.add(self._index([p[0] + w, p[1] , p[2] , p[3]])) self._nodes.add(self._index([p[0] , p[1] + w, p[2] , p[3]])) self._nodes.add(self._index([p[0] + w, p[1] + w, p[2] , p[3]])) self._nodes.add(self._index([p[0] , p[1] , p[2] + w, p[3]])) self._nodes.add(self._index([p[0] + w, p[1] , p[2] + w, p[3]])) self._nodes.add(self._index([p[0] , p[1] + w, p[2] + w, p[3]])) self._nodes.add(self._index([p[0] + w, p[1] + w, p[2] + w, p[3]])) gridN = [] self._n2i = dict() for ii, n in enumerate(sorted(self._nodes)): self._n2i[n] = ii gridN.append( self._cellN( self._pointer(n)[:-1] ) ) self._gridN = np.array(gridN) self.__dirtyNodes__ = False def _numberFaces(self, force=False): if not self.__dirtyFaces__ and not force: return self._facesX = set() self._facesY = set() if self.dim == 3: self._facesZ = set() for ind in self._treeInds: p = self._asPointer(ind) w = self._levelWidth(p[-1]) if self.dim == 2: self._facesX.add(self._index([p[0] , p[1] , p[2]])) self._facesX.add(self._index([p[0] + w, p[1] , p[2]])) self._facesY.add(self._index([p[0] , p[1] , p[2]])) self._facesY.add(self._index([p[0] , p[1] + w, p[2]])) elif self.dim == 3: self._facesX.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._facesX.add(self._index([p[0] + w, p[1] , p[2] , p[3]])) self._facesY.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._facesY.add(self._index([p[0] , p[1] + w, p[2] , p[3]])) self._facesZ.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._facesZ.add(self._index([p[0] , p[1] , p[2] + w, p[3]])) gridFx = [] self._fx2i = dict() for ii, fx in enumerate(sorted(self._facesX)): self._fx2i[fx] = ii p = self._pointer(fx) n, h = self._cellN(p), self._cellH(p) if self.dim == 2: gridFx.append( [n[0], n[1] + h[1]/2.0] ) elif self.dim == 3: gridFx.append( [n[0], n[1] + h[1]/2.0, n[2] + h[2]/2.0] ) self._gridFx = np.array(gridFx) gridFy = [] self._fy2i = dict() for ii, fy in enumerate(sorted(self._facesY)): self._fy2i[fy] = ii p = self._pointer(fy) n, h = self._cellN(p), self._cellH(p) if self.dim == 2: gridFy.append( [n[0] + h[0]/2.0, n[1]] ) elif self.dim == 3: gridFy.append( [n[0] + h[0]/2.0, n[1], n[2] + h[2]/2.0] ) self._gridFy = np.array(gridFy) if self.dim == 2: self.__dirtyFaces__ = False return gridFz = [] self._fz2i = dict() for ii, fz in enumerate(sorted(self._facesZ)): self._fz2i[fz] = ii p = self._pointer(fz) n, h = self._cellN(p), self._cellH(p) gridFz.append( [n[0] + h[0]/2.0, n[1] + h[1]/2.0, n[2]] ) self._gridFz = np.array(gridFz) self.__dirtyFaces__ = False def _hanging(self, force=False): if not self.__dirtyHanging__ and not force: return self._numberNodes() self._numberFaces() self._numberEdges() self._hangingN = dict() self._hangingFx = dict() self._hangingFy = dict() if self.dim == 3: self._hangingFz = dict() self._hangingEx = dict() self._hangingEy = dict() self._hangingEz = dict() # Compute from x faces for fx in self._facesX: p = self._pointer(fx) if p[-1] + 1 > self.levels: continue sl = p[-1] + 1 #: small level test = self._index(p[:-1] + [sl]) if test not in self._facesX: # Return early without checking the other faces continue w = self._levelWidth(sl) if self.dim == 2: self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], 0.5], ) self._hangingFx[self._fx2i[self._index([p[0] , p[1] + w, sl])]] = ([self._fx2i[fx], 0.5], ) n0, n1 = fx, self._index([p[0], p[1] + 2*w, p[-1]]) self._hangingN[self._n2i[test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[self._n2i[self._index([p[0] , p[1] + w, sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) self._hangingN[self._n2i[self._index([p[0] , p[1] + 2*w, sl])]] = ([self._n2i[n1], 1.0], ) elif self.dim == 3: ey0 = fx ey1 = self._index([p[0], p[1] , p[2] + 2*w, p[-1]]) ez0 = fx ez1 = self._index([p[0], p[1] + 2*w, p[2] , p[-1]]) n0 = fx n1 = self._index([p[0], p[1] + 2*w, p[2] , p[-1]]) n2 = self._index([p[0], p[1] , p[2] + 2*w, p[-1]]) n3 = self._index([p[0], p[1] + 2*w, p[2] + 2*w, p[-1]]) self._hangingFx[self._fx2i[test ]] = ([self._fx2i[fx], 0.25], ) self._hangingFx[self._fx2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._fx2i[fx], 0.25], ) self._hangingFx[self._fx2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._fx2i[fx], 0.25], ) self._hangingFx[self._fx2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._fx2i[fx], 0.25], ) self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], 0.5], ) self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ey2i[ey0], 0.5], ) self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) self._hangingEy[self._ey2i[self._index([p[0], p[1] , p[2] + 2*w, sl])]] = ([self._ey2i[ey1], 0.5], ) self._hangingEy[self._ey2i[self._index([p[0], p[1] + w, p[2] + 2*w, sl])]] = ([self._ey2i[ey1], 0.5], ) self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], 0.5], ) self._hangingEz[self._ez2i[self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._ez2i[ez0], 0.5], ) self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) self._hangingEz[self._ez2i[self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] , sl])]] = ([self._ez2i[ez1], 0.5], ) self._hangingEz[self._ez2i[self._index([p[0], p[1] + 2*w, p[2] + w, sl])]] = ([self._ez2i[ez1], 0.5], ) self._hangingN[ self._n2i[ test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0], p[1] + w, p[2] , sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) self._hangingN[ self._n2i[ self._index([p[0], p[1] + 2*w, p[2] , sl])]] = ([self._n2i[n1], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0], p[1] , p[2] + w, sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n2], 0.5]) self._hangingN[ self._n2i[ self._index([p[0], p[1] + w, p[2] + w, sl])]] = ([self._n2i[n0], 0.25], [self._n2i[n1], 0.25], [self._n2i[n2], 0.25], [self._n2i[n3], 0.25]) self._hangingN[ self._n2i[ self._index([p[0], p[1] + 2*w, p[2] + w, sl])]] = ([self._n2i[n1], 0.5], [self._n2i[n3], 0.5]) self._hangingN[ self._n2i[ self._index([p[0], p[1] , p[2] + 2*w, sl])]] = ([self._n2i[n2], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0], p[1] + w, p[2] + 2*w, sl])]] = ([self._n2i[n2], 0.5], [self._n2i[n3], 0.5]) self._hangingN[ self._n2i[ self._index([p[0], p[1] + 2*w, p[2] + 2*w, sl])]] = ([self._n2i[n3], 1.0], ) # Compute from y faces for fy in self._facesY: p = self._pointer(fy) if p[-1] + 1 > self.levels: continue sl = p[-1] + 1 #: small level test = self._index(p[:-1] + [sl]) if test not in self._facesY: # Return early without checking the other faces continue w = self._levelWidth(sl) if self.dim == 2: self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], 0.5], ) self._hangingFy[self._fy2i[self._index([p[0] + w, p[1] , sl])]] = ([self._fy2i[fy], 0.5], ) n0, n1 = fy, self._index([p[0] + 2*w, p[1], p[-1]]) self._hangingN[self._n2i[test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[self._n2i[self._index([p[0] + w, p[1] , sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) self._hangingN[self._n2i[self._index([p[0] + 2*w, p[1] , sl])]] = ([self._n2i[n1], 1.0], ) elif self.dim == 3: ex0 = fy ex1 = self._index([p[0] , p[1], p[2] + 2*w, p[-1]]) ez0 = fy ez1 = self._index([p[0] + 2*w, p[1], p[2] , p[-1]]) n0 = fy n1 = self._index([p[0] + 2*w, p[1], p[2] , p[-1]]) n2 = self._index([p[0] , p[1], p[2] + 2*w, p[-1]]) n3 = self._index([p[0] + 2*w, p[1], p[2] + 2*w, p[-1]]) self._hangingFy[self._fy2i[test ]] = ([self._fy2i[fy], 0.25], ) self._hangingFy[self._fy2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._fy2i[fy], 0.25], ) self._hangingFy[self._fy2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._fy2i[fy], 0.25], ) self._hangingFy[self._fy2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._fy2i[fy], 0.25], ) self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], 0.5], ) self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ex2i[ex0], 0.5], ) self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) self._hangingEx[self._ex2i[self._index([p[0] , p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], 0.5], ) self._hangingEx[self._ex2i[self._index([p[0] + w, p[1], p[2] + 2*w, sl])]] = ([self._ex2i[ex1], 0.5], ) self._hangingEz[self._ez2i[test ]] = ([self._ez2i[ez0], 0.5], ) self._hangingEz[self._ez2i[self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], 0.5], ) self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) self._hangingEz[self._ez2i[self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez0], 0.25], [self._ez2i[ez1], 0.25]) self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] , sl])]] = ([self._ez2i[ez1], 0.5], ) self._hangingEz[self._ez2i[self._index([p[0] + 2*w, p[1], p[2] + w, sl])]] = ([self._ez2i[ez1], 0.5], ) self._hangingN[ self._n2i[ test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1], p[2] , sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] + 2*w, p[1], p[2] , sl])]] = ([self._n2i[n1], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] , p[1], p[2] + w, sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n2], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1], p[2] + w, sl])]] = ([self._n2i[n0], 0.25], [self._n2i[n1], 0.25], [self._n2i[n2], 0.25], [self._n2i[n3], 0.25]) self._hangingN[ self._n2i[ self._index([p[0] + 2*w, p[1], p[2] + w, sl])]] = ([self._n2i[n1], 0.5], [self._n2i[n3], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] , p[1], p[2] + 2*w, sl])]] = ([self._n2i[n2], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1], p[2] + 2*w, sl])]] = ([self._n2i[n2], 0.5], [self._n2i[n3], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] + 2*w, p[1], p[2] + 2*w, sl])]] = ([self._n2i[n3], 1.0], ) if self.dim == 2: self.__dirtyHanging__ = False return # Compute from z faces for fz in self._facesZ: p = self._pointer(fz) if p[-1] + 1 > self.levels: continue sl = p[-1] + 1 #: small level test = self._index(p[:-1] + [sl]) if test not in self._facesZ: # Return early without checking the other faces continue w = self._levelWidth(sl) ex0 = fz ex1 = self._index([p[0] , p[1] + 2*w, p[2], p[-1]]) ey0 = fz ey1 = self._index([p[0] + 2*w, p[1] , p[2], p[-1]]) n0 = fz n1 = self._index([p[0] + 2*w, p[1] , p[2], p[-1]]) n2 = self._index([p[0] , p[1] + 2*w, p[2], p[-1]]) n3 = self._index([p[0] + 2*w, p[1] + 2*w, p[2], p[-1]]) self._hangingFz[self._fz2i[test ]] = ([self._fz2i[fz], 0.25], ) self._hangingFz[self._fz2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._fz2i[fz], 0.25], ) self._hangingFz[self._fz2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._fz2i[fz], 0.25], ) self._hangingFz[self._fz2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._fz2i[fz], 0.25], ) self._hangingEx[self._ex2i[test ]] = ([self._ex2i[ex0], 0.5], ) self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ex2i[ex0], 0.5], ) self._hangingEx[self._ex2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ex2i[ex0], 0.25], [self._ex2i[ex1], 0.25]) self._hangingEx[self._ex2i[self._index([p[0] , p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], 0.5], ) self._hangingEx[self._ex2i[self._index([p[0] + w, p[1] + 2*w, p[2], sl])]] = ([self._ex2i[ex1], 0.5], ) self._hangingEy[self._ey2i[test ]] = ([self._ey2i[ey0], 0.5], ) self._hangingEy[self._ey2i[self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], 0.5], ) self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) self._hangingEy[self._ey2i[self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey0], 0.25], [self._ey2i[ey1], 0.25]) self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] , p[2], sl])]] = ([self._ey2i[ey1], 0.5], ) self._hangingEy[self._ey2i[self._index([p[0] + 2*w, p[1] + w, p[2], sl])]] = ([self._ey2i[ey1], 0.5], ) self._hangingN[ self._n2i[ test ]] = ([self._n2i[n0], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1] , p[2], sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n1], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] + 2*w, p[1] , p[2], sl])]] = ([self._n2i[n1], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] , p[1] + w, p[2], sl])]] = ([self._n2i[n0], 0.5], [self._n2i[n2], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1] + w, p[2], sl])]] = ([self._n2i[n0], 0.25], [self._n2i[n1], 0.25], [self._n2i[n2], 0.25], [self._n2i[n3], 0.25]) self._hangingN[ self._n2i[ self._index([p[0] + 2*w, p[1] + w, p[2], sl])]] = ([self._n2i[n1], 0.5], [self._n2i[n3], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] , p[1] + 2*w, p[2], sl])]] = ([self._n2i[n2], 1.0], ) self._hangingN[ self._n2i[ self._index([p[0] + w, p[1] + 2*w, p[2], sl])]] = ([self._n2i[n2], 0.5], [self._n2i[n3], 0.5]) self._hangingN[ self._n2i[ self._index([p[0] + 2*w, p[1] + 2*w, p[2], sl])]] = ([self._n2i[n3], 1.0], ) self.__dirtyHanging__ = False def _numberEdges(self, force=False): if self.dim == 2: return if not self.__dirtyEdges__ and not force: return self._edgesX = set() self._edgesY = set() self._edgesZ = set() for ind in self._treeInds: p = self._asPointer(ind) w = self._levelWidth(p[-1]) self._edgesX.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._edgesX.add(self._index([p[0] , p[1] + w, p[2] , p[3]])) self._edgesX.add(self._index([p[0] , p[1] , p[2] + w, p[3]])) self._edgesX.add(self._index([p[0] , p[1] + w, p[2] + w, p[3]])) self._edgesY.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._edgesY.add(self._index([p[0] + w, p[1] , p[2] , p[3]])) self._edgesY.add(self._index([p[0] , p[1] , p[2] + w, p[3]])) self._edgesY.add(self._index([p[0] + w, p[1] , p[2] + w, p[3]])) self._edgesZ.add(self._index([p[0] , p[1] , p[2] , p[3]])) self._edgesZ.add(self._index([p[0] + w, p[1] , p[2] , p[3]])) self._edgesZ.add(self._index([p[0] , p[1] + w, p[2] , p[3]])) self._edgesZ.add(self._index([p[0] + w, p[1] + w, p[2] , p[3]])) gridEx = [] self._ex2i = dict() for ii, ex in enumerate(sorted(self._edgesX)): self._ex2i[ex] = ii p = self._pointer(ex) n, h = self._cellN(p), self._cellH(p) gridEx.append( [n[0] + h[0]/2.0, n[1], n[2]] ) self._gridEx = np.array(gridEx) gridEy = [] self._ey2i = dict() for ii, ey in enumerate(sorted(self._edgesY)): self._ey2i[ey] = ii p = self._pointer(ey) n, h = self._cellN(p), self._cellH(p) gridEy.append( [n[0], n[1] + h[1]/2.0, n[2]] ) self._gridEy = np.array(gridEy) gridEz = [] self._ez2i = dict() for ii, ez in enumerate(sorted(self._edgesZ)): self._ez2i[ez] = ii p = self._pointer(ez) n, h = self._cellN(p), self._cellH(p) gridEz.append( [n[0], n[1], n[2] + h[2]/2.0] ) self._gridEz = np.array(gridEz) self.__dirtyEdges__ = False def number(self, force=False): if not self.__dirty__ and not force: return self._hanging() return facesX, facesY, facesZ = [], [], [] areaX, areaY, areaZ = [], [], [] hangingFacesX, hangingFacesY, hangingFacesZ = [], [], [] hangingNodes = [] faceXCount, faceYCount, faceZCount = -1, -1, -1 nodeCount = -1 fXm,fXp,fYm,fYp,fZm,fZp = range(6) vol, nodes = [], [] def addXFace(count, p, positive=True): n = self._cellN(p) w = self._cellH(p) areaX.append(w[1] if self.dim == 2 else w[1]*w[2]) if self.dim == 2: facesX.append([n[0] + (w[0] if positive else 0), n[1] + w[1]/2.0]) elif self.dim == 3: facesX.append([n[0] + (w[0] if positive else 0), n[1] + w[1]/2.0, n[2] + w[2]/2.0]) return count + 1 def addYFace(count, p, positive=True): n = self._cellN(p) w = self._cellH(p) areaY.append(w[0] if self.dim == 2 else w[0]*w[2]) if self.dim == 2: facesY.append([n[0] + w[0]/2.0, n[1] + (w[1] if positive else 0)]) elif self.dim == 3: facesY.append([n[0] + w[0]/2.0, n[1] + (w[1] if positive else 0), n[2] + w[2]/2.0]) return count + 1 def addZFace(count, p, positive=True): n = self._cellN(p) w = self._cellH(p) areaZ.append(w[0]*w[1]) facesZ.append([n[0] + w[0]/2.0, n[1] + w[1]/2.0, n[2] + (w[2] if positive else 0)]) return count + 1 def addNode(count, p, loc=[0,0,0]): """loc=[0,0]""" n = self._cellN(p) w = self._cellH(p) if self.dim == 2: nodes.append([n[0] + w[0]*loc[0], n[1] + w[1]*loc[1]]) elif self.dim == 3: nodes.append([n[0] + w[0]*loc[0], n[1] + w[1]*loc[1], n[2] + w[2]*loc[2]]) return count + 1 # c2cn = dict() c2f = dict() def gc2f(ind): if ind in c2f: return c2f[ind] c2f_ind = [list() for _ in xrange(2*self.dim)] c2f[ind] = c2f_ind return c2f_ind c2n = dict() def gc2n(ind): if ind in c2n: return c2n[ind] c2n_ind = [list() for _ in xrange(2**self.dim)] c2n[ind] = c2n_ind return c2n_ind def processCellFace(ind, faceCount, addFace, hangingFaces, DIR=0): fM,fP=(0,1) if DIR == 0 else (2,3) if DIR == 1 else (4,5) p = self._asPointer(ind) if self._getNextCell(p, direction=DIR, positive=False) is None: faceCount = addFace(faceCount, p, positive=False) gc2f(ind)[fM] += [faceCount] nextCell = self._getNextCell(p, direction=DIR) # Add the next Xface if nextCell is None: # on the boundary faceCount = addFace(faceCount, p) gc2f(ind)[fP] += [faceCount] elif type(nextCell) in [int, long] and self._onSameLevel(p,nextCell): # same sized cell faceCount = addFace(faceCount, p) gc2f(ind)[fP] += [faceCount] gc2f(nextCell)[fM] += [faceCount] elif type(nextCell) in [int, long] and not self._onSameLevel(p,nextCell): # the cell is bigger than me faceCount = addFace(faceCount, p) gc2f(ind)[fP] += [faceCount] gc2f(nextCell)[fM] += [faceCount] hangingFaces.append(faceCount) elif type(nextCell) is list: # the cell is smaller than me # TODO: ensure that things are balanced. p0 = self._pointer(nextCell[0]) p1 = self._pointer(nextCell[1]) faceCount = addFace(faceCount, p0, positive=False) gc2f(nextCell[0])[fM] += [faceCount] faceCount = addFace(faceCount, p1, positive=False) gc2f(nextCell[1])[fM] += [faceCount] gc2f(ind)[fP] += [faceCount-1,faceCount] hangingFaces += [faceCount-1, faceCount] return faceCount def processCellNode(ind, nodeCount): MMM, PMM, MPM, PPM, MMP, PMP, MPP, PPP = range(8) p = self._asPointer(ind) xM = self._getNextCell(p, direction=0, positive=False) yM = self._getNextCell(p, direction=1, positive=False) zM = None if self.dim == 2 else self._getNextCell(p, direction=2, positive=False) xP = self._getNextCell(p, direction=0, positive=True) yP = self._getNextCell(p, direction=1, positive=True) zP = None if self.dim == 2 else self._getNextCell(p, direction=2, positive=True) if xM is None and yM is None and zM is None: nodeCount = addNode(nodeCount, p, loc=[0,0,0]) gc2n(ind)[MMM] += [nodeCount] if yM is None: nodeCount = addNode(nodeCount, p, loc=[1,0,0]) gc2n(ind)[PMM] += [nodeCount] if xM is None: nodeCount = addNode(nodeCount, p, loc=[0,1,0]) gc2n(ind)[MPM] += [nodeCount] # Add the next Xface if nextCell is None: # on the boundary pass # nodeCount = addFace(nodeCount, p) # gc2f(ind)[fP] += [nodeCount] elif type(nextCell) in [int, long] and self._onSameLevel(p,nextCell): # same sized cell pass # nodeCount = addFace(nodeCount, p) # gc2f(ind)[fP] += [nodeCount] # gc2f(nextCell)[fM] += [nodeCount] elif type(nextCell) in [int, long] and not self._onSameLevel(p,nextCell): # the cell is bigger than me pass # nodeCount = addFace(nodeCount, p) # gc2f(ind)[fP] += [nodeCount] # gc2f(nextCell)[fM] += [nodeCount] # hangingFaces.append(nodeCount) elif type(nextCell) is list: # the cell is smaller than me pass # TODO: ensure that things are balanced. # p0 = self._pointer(nextCell[0]) # p1 = self._pointer(nextCell[1]) # nodeCount = addFace(nodeCount, p0, positive=False) # gc2f(nextCell[0])[fM] += [nodeCount] # nodeCount = addFace(nodeCount, p1, positive=False) # gc2f(nextCell[1])[fM] += [nodeCount] # gc2f(ind)[fP] += [nodeCount-1,nodeCount] # hangingFaces += [nodeCount-1, nodeCount] return nodeCount for ii, ind in enumerate(self._sortedInds): # c2cn[ind] = ii vol.append(np.prod(self._cellH(ind))) # nodeCount = processCellNode(ind, nodeCount) faceXCount = processCellFace(ind, faceXCount, addXFace, hangingFacesX, DIR=0) faceYCount = processCellFace(ind, faceYCount, addYFace, hangingFacesY, DIR=1) if self.dim == 3: faceZCount = processCellFace(ind, faceZCount, addZFace, hangingFacesZ, DIR=2) self._c2f = c2f self._area = np.array(areaX + areaY + (areaZ if self.dim == 3 else [])) self._vol = np.array(vol) self._gridFx = np.array(facesX) self._gridFy = np.array(facesY) self._gridN = np.array(nodes) self._hangingFx = hangingFacesX self._hangingFy = hangingFacesY if self.dim == 3: self._gridFz = np.array(facesZ) self._nFz = self._gridFz.shape[0] self._hangingFz = hangingFacesZ self._nC = len(self._sortedInds) self._nN = self._gridN.shape[0] self._nFx = self._gridFx.shape[0] self._nFy = self._gridFy.shape[0] self._nF = self._nFx + self._nFy + (self._nFz if self.dim == 3 else 0) # self.__dirty__ = False def _deflationMatrix(self, theSet, theHang, theIndex): reducedInd = dict() # final reduced index ii = 0 I,J,V = [],[],[] for fx in sorted(theSet): if theIndex[fx] not in theHang: reducedInd[theIndex[fx]] = ii I += [theIndex[fx]] J += [ii] V += [1.0] ii += 1 for hfkey in theHang.keys(): hf = theHang[hfkey] I += [hfkey]*len(hf) J += [_[0] for _ in hf] V += [_[1] for _ in hf] return sp.csr_matrix((V,(I,J)), shape=(len(theSet), len(reducedInd))) @property def faceDiv(self): # print self._c2f if getattr(self, '_faceDiv', None) is None: self.number() # TODO: Preallocate! I, J, V = [], [], [] PM = [-1,1]*self.dim # plus / minus offset = [0,0,self.nFx,self.nFx,self.nFx+self.nFy,self.nFx+self.nFy] for ii, ind in enumerate(self._sortedInds): p = self._pointer(ind) w = self._levelWidth(p[-1]) if self.dim == 2: faces = [ self._fx2i[self._index([ p[0] , p[1] , p[2]])], self._fx2i[self._index([ p[0] + w, p[1] , p[2]])], self._fy2i[self._index([ p[0] , p[1] , p[2]])], self._fy2i[self._index([ p[0] , p[1] + w, p[2]])] ] elif self.dim == 3: faces = [ self._fx2i[self._index([ p[0] , p[1] , p[2] , p[3]])], self._fx2i[self._index([ p[0] + w, p[1] , p[2] , p[3]])], self._fy2i[self._index([ p[0] , p[1] , p[2] , p[3]])], self._fy2i[self._index([ p[0] , p[1] + w, p[2] , p[3]])], self._fz2i[self._index([ p[0] , p[1] , p[2] , p[3]])], self._fz2i[self._index([ p[0] , p[1] , p[2] + w, p[3]])] ] for off, pm, face in zip(offset,PM,faces): I += [ii] J += [face + off] V += [pm] # total number of faces tnF = len(self._facesX) + len(self._facesY) + (0 if self.dim == 2 else len(self._facesZ)) D = sp.csr_matrix((V,(I,J)), shape=(self.nC, tnF)) Rlist = [0]*self.dim Rlist[0] = self._deflationMatrix(self._facesX, self._hangingFx, self._fx2i) Rlist[1] = self._deflationMatrix(self._facesY, self._hangingFy, self._fy2i) if self.dim == 3: Rlist[2] = self._deflationMatrix(self._facesZ, self._hangingFz, self._fz2i) R = sp.block_diag(Rlist) # VOL = self.vol # S = self.area self._faceDiv = D * R # self._faceDiv = Utils.sdiag(1.0/VOL)*D*Utils.sdiag(S) return self._faceDiv def plotGrid(self, ax=None, showIt=False, grid=True): self.number() axOpts = {'projection':'3d'} if self.dim == 3 else {} if ax is None: ax = plt.subplot(111, **axOpts) else: assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" fig = ax.figure if grid: for ind in self._sortedInds: p = self._asPointer(ind) n = self._cellN(p) h = self._cellH(p) x = [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0]] y = [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1]] if self.dim == 2: ax.plot(x,y, 'b-') elif self.dim == 3: ax.plot(x,y, 'b-', zs=[n[2]]*5) z = [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2]] ax.plot(x,y, 'b-', zs=z) sides = [0,0], [h[0],0], [0,h[1]], [h[0],h[1]] for s in sides: x = [n[0] + s[0], n[0] + s[0]] y = [n[1] + s[1], n[1] + s[1]] z = [n[2] , n[2] + h[2]] ax.plot(x,y, 'b-', zs=z) if self.dim == 2: ax.plot(self.gridCC[[0,-1],0], self.gridCC[[0,-1],1], 'ro') ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.') ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r:') 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.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.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: ax.plot(self.gridCC[[0,-1],0], self.gridCC[[0,-1],1], 'ro', zs=self.gridCC[[0,-1],2]) ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.', zs=self.gridCC[:,2]) ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r:', zs=self.gridCC[:,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.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.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.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.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.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.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.axis('equal') if showIt:plt.show() if __name__ == '__main__': def function(xc): r = xc - np.r_[0.5,0.5] dist = np.sqrt(r.dot(r)) # if dist < 0.05: # return 5 if dist < 0.1: return 4 if dist < 0.3: return 3 if dist < 1.0: return 2 else: return 0 # T = Tree([[(1,8)],[(1,8)],[(1,8)]],levels=3) T = Tree([[(1,16)],[(1,16)]],levels=4) T.refine(lambda xc:1) T._refineCell([0,0,1]) T._refineCell([8,8,1]) # T._refineCell([4,4,4,1]) ax = plt.subplot(211) ax.spy(T.faceDiv) # R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i) # print R ax = plt.subplot(212) # ax.spy(R) # ax = plt.subplot(313) # ax.spy(T.faceDiv[:,:T.nFx] * R) T.plotGrid(ax=ax) # print T.nN plt.show()