From 310f327b5ad1313c63242d7acb639e32067bf437 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 9 Feb 2015 21:30:54 -0800 Subject: [PATCH] updates to edges and areas --- SimPEG/Mesh/NewTreeMesh.py | 120 ++++++++++++++++++++++++------------- 1 file changed, 77 insertions(+), 43 deletions(-) diff --git a/SimPEG/Mesh/NewTreeMesh.py b/SimPEG/Mesh/NewTreeMesh.py index f7282c10..a8741c0f 100644 --- a/SimPEG/Mesh/NewTreeMesh.py +++ b/SimPEG/Mesh/NewTreeMesh.py @@ -2,9 +2,10 @@ import numpy as np, scipy.sparse as sp from SimPEG.Utils import ndgrid, mkvc -NUM, PARENT, ACTIVE, EDIR, ENODE0, ENODE1 = range(6) -NUM, PARENT, ACTIVE, FDIR, FEDGE0, FEDGE1, FEDGE2, FEDGE3 = range(8) -NUM, PARENT, ACTIVE, CFACE0, CFACE1, CFACE2, CFACE3, CFACE4, CFACE5 = range(9) +NUM, ACTIVE, NX, NY, NZ = range(5) +NUM, ACTIVE, PARENT, EDIR, ENODE0, ENODE1 = range(6) +NUM, ACTIVE, PARENT, FDIR, FEDGE0, FEDGE1, FEDGE2, FEDGE3 = range(8) +NUM, ACTIVE, PARENT, CFACE0, CFACE1, CFACE2, CFACE3, CFACE4, CFACE5 = range(9) def SortByX0(grid): dtype=[('x',float),('y',float)] @@ -24,9 +25,7 @@ class TreeMesh(object): vnN = [nx.size, ny.size] XY = ndgrid(nx, ny) - N = np.c_[np.arange(XY.shape[0]), XY] - - N.astype([('num',int),('x',float),('y',float),('z',float)]) + N = np.c_[np.arange(XY.shape[0]), np.ones(XY.shape[0]), XY] I = np.arange(nx.size * ny.size, dtype=int).reshape(vnN, order='F') @@ -39,13 +38,13 @@ class TreeMesh(object): zEx = np.zeros(nEx.size, dtype=int) zEy = np.zeros(nEy.size, dtype=int) - # # parent active dir, n1,n2 - Ex = np.c_[mkvc(nEx), zEx-1, zEx+1, zEx+0, vEx] - Ey = np.c_[mkvc(nEy), zEy-1, zEy+1, zEy+1, vEy] + # # active parent dir, n1,n2 + Ex = np.c_[mkvc(nEx), zEx+1, zEx-1, zEx+0, vEx] + Ey = np.c_[mkvc(nEy), zEy+1, zEy-1, zEy+1, vEy] nC = np.arange(np.prod(vnC), dtype=int) - C = np.c_[nC, nC*0-1, nC*0+1, nC*0+2, mkvc(nEx[:,:-1]), mkvc(nEx[:,1:]), mkvc(nEy[:-1,:]), mkvc(nEy[1:,:])] + C = np.c_[nC, nC*0+1, nC*0-1, nC*0+2, mkvc(nEx[:,:-1]), mkvc(nEx[:,1:]), mkvc(nEy[:-1,:]), mkvc(nEy[1:,:])] self._nodes = N self._edges = np.r_[Ex, Ey] @@ -55,13 +54,16 @@ class TreeMesh(object): @property def isNumbered(self): - return self._numberedCC and self._numberedFx and self._numberedFy + return self._numberedCC and self._numberedEx and self._numberedEy @isNumbered.setter def isNumbered(self, value): assert value is False, 'Can only set to False.' self._numberedCC = False self._numberedEx = False self._numberedEy = False + for name in ['vol', 'area', 'edge', 'gridCC', 'gridN', 'gridEx', 'gridEy', 'gridEz', 'gridFx', 'gridFy', 'gridFz']: + if hasattr(self, '_'+name): + delattr(self, '_'+name) @property def dim(self): @@ -84,6 +86,7 @@ class TreeMesh(object): between = np.array(between).flatten() nodes = self._nodes[between.astype(int), :] newNode = np.mean(nodes, axis=0) + newNode[ACTIVE] = 1 return self._push('_nodes', newNode) def refineEdge(self, index): @@ -162,6 +165,10 @@ class TreeMesh(object): return np.sum(self._faces[:,ACTIVE] == 1) return np.sum(self._cells[:,ACTIVE] == 1) + @property + def nN(self): + return np.sum(self._cells[:,ACTIVE] == 1) + @property def nE(self): return np.sum(self._edges[:,ACTIVE] == 1) @@ -205,16 +212,29 @@ class TreeMesh(object): return np.sum((self._faces[:,ACTIVE] == 1) & (self._faces[:,FDIR] == 1)) @property - def area(self): - if getattr(self, '_area', None) is None: + def edge(self): + if getattr(self, '_edge', None) is None: + self.number() N = self._nodes E = self._edges activeEdges = E[:,ACTIVE] == 1 - e0xy = N[E[activeEdges,ENODE0],:][:,[1,2]] - e1xy = N[E[activeEdges,ENODE1],:][:,[1,2]] + e0xy = N[E[activeEdges,ENODE0],:][:,[NX,NY]] + e1xy = N[E[activeEdges,ENODE1],:][:,[NX,NY]] - self._area = np.sum((e1xy - e0xy)**2,axis=1)**0.5 + A = np.sum((e1xy - e0xy)**2,axis=1)**0.5 + + P = np.argsort(E[activeEdges,NUM]) + self._edge = A[P] + + return self._edge + + @property + def area(self): + if getattr(self, '_area', None) is None: + self.number() + if self.dim == 2: + self._area = np.r_[self.edge[self.nEx:], self.edge[:self.nEx]] return self._area @@ -222,6 +242,7 @@ class TreeMesh(object): @property def vol(self): if getattr(self, '_vol', None) is None: + self.number() N = self._nodes E = self._edges @@ -229,10 +250,10 @@ class TreeMesh(object): activeCells = C[:,ACTIVE] == 1 nInds1 = E[C[activeCells,FEDGE0],:][:,[ENODE0,ENODE1]] nInds2 = E[C[activeCells,FEDGE1],:][:,[ENODE0,ENODE1]] - n0 = N[nInds1[:,0],:][:,[1,2]] # 2------3 3------2 - n1 = N[nInds1[:,1],:][:,[1,2]] # | | --> | | - n3 = N[nInds2[:,0],:][:,[1,2]] # | | | | - n2 = N[nInds2[:,1],:][:,[1,2]] # 0------1 0------1 + n0 = N[nInds1[:,0],:][:,[NX,NY]] # 2------3 3------2 + n1 = N[nInds1[:,1],:][:,[NX,NY]] # | | --> | | + n3 = N[nInds2[:,0],:][:,[NX,NY]] # | | | | + n2 = N[nInds2[:,1],:][:,[NX,NY]] # 0------1 0------1 a = np.sum((n1 - n0)**2,axis=1)**0.5 b = np.sum((n2 - n1)**2,axis=1)**0.5 @@ -256,10 +277,10 @@ class TreeMesh(object): activeCells = C[:,ACTIVE] == 1 nInds1 = E[C[activeCells,FEDGE0],:][:,[ENODE0,ENODE1]] nInds2 = E[C[activeCells,FEDGE1],:][:,[ENODE0,ENODE1]] - Cx = (N[nInds1[:,0],1] + N[nInds1[:,1],1] + N[nInds2[:,0],1] + N[nInds2[:,1],1])/4.0 - Cy = (N[nInds1[:,0],2] + N[nInds1[:,1],2] + N[nInds2[:,0],2] + N[nInds2[:,1],2])/4.0 + Cx = (N[nInds1[:,0],NX] + N[nInds1[:,1],NX] + N[nInds2[:,0],NX] + N[nInds2[:,1],NX])/4.0 + Cy = (N[nInds1[:,0],NY] + N[nInds1[:,1],NY] + N[nInds2[:,0],NY] + N[nInds2[:,1],NY])/4.0 - P = SortByX0(np.c_[N[nInds1[:,0],1], N[nInds1[:,0],2]]) + P = SortByX0(np.c_[N[nInds1[:,0],NX], N[nInds1[:,0],NY]]) if not self._numberedCC: cnt = np.zeros(P.size, dtype=int) cnt[P] = np.arange(P.size) @@ -275,10 +296,10 @@ class TreeMesh(object): C = self._faces activeEdges = (E[:,ACTIVE] == 1) & (E[:,EDIR] == 0) nInds = E[activeEdges,:][:,[ENODE0,ENODE1]] - Ex = (N[nInds[:,0],1] + N[nInds[:,1],1])/2.0 - Ey = (N[nInds[:,0],2] + N[nInds[:,1],2])/2.0 + Ex = (N[nInds[:,0],NX] + N[nInds[:,1],NX])/2.0 + Ey = (N[nInds[:,0],NY] + N[nInds[:,1],NY])/2.0 - P = SortByX0(np.c_[N[nInds[:,0],1], N[nInds[:,0],2]]) + P = SortByX0(np.c_[N[nInds[:,0],NX], N[nInds[:,0],NY]]) if not self._numberedEx: cnt = np.zeros(P.size, dtype=int) cnt[P] = np.arange(P.size) @@ -294,10 +315,10 @@ class TreeMesh(object): C = self._faces activeEdges = (E[:,ACTIVE] == 1) & (E[:,EDIR] == 1) nInds = E[activeEdges,:][:,[ENODE0,ENODE1]] - Ex = (N[nInds[:,0],1] + N[nInds[:,1],1])/2.0 - Ey = (N[nInds[:,0],2] + N[nInds[:,1],2])/2.0 + Ex = (N[nInds[:,0],NX] + N[nInds[:,1],NX])/2.0 + Ey = (N[nInds[:,0],NY] + N[nInds[:,1],NY])/2.0 - P = SortByX0(np.c_[N[nInds[:,0],1], N[nInds[:,0],2]]) + P = SortByX0(np.c_[N[nInds[:,0],NX], N[nInds[:,0],NY]]) if not self._numberedEy: cnt = np.zeros(P.size, dtype=int) cnt[P] = np.arange(P.size) @@ -306,6 +327,17 @@ class TreeMesh(object): return np.c_[Ex,Ey][P, :] + + @property + def gridFx(self): + if self.dim == 2: + return self.gridEy + + @property + def gridFy(self): + if self.dim == 2: + return self.gridEx + def _index(self, attr, index): index = [index] if np.isscalar(index) else list(index) C = getattr(self, attr) @@ -328,20 +360,15 @@ class TreeMesh(object): self.number() # TODO: Preallocate! I, J, V = [], [], [] - nEx, nFx = self.nEx, self.nFx - offset = np.r_[nFx, -nEx] - N = self._nodes - E = self._edges + offset = np.r_[self.nFx, -self.nEx] # this switches from edge to face numbering C = self._faces activeCells = C[:,ACTIVE] == 1 for cell in C[activeCells]: for sign, face in zip([-1,1,-1,1],[FEDGE0, FEDGE1, FEDGE2, FEDGE3]): ij, jrow = self._index('_edges', cell[face]) I += [cell[NUM]]*len(ij) - print jrow - J += list(jrow[:,0] + offset[jrow[:,EDIR]])# + nFx) - # J += list(jrow[:,0] - nEx) + J += list(jrow[:,0] + offset[jrow[:,EDIR]]) V += [sign]*len(ij) VOL = self.vol D = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) @@ -350,6 +377,8 @@ class TreeMesh(object): return self._faceDiv def number(self): + if self.isNumbered: + return self._nodes[:,NUM] = -1 self._edges[:,NUM] = -1 self._faces[:,NUM] = -1 @@ -373,16 +402,16 @@ class TreeMesh(object): activeCells = C[:,ACTIVE] == 1 for FEDGE in [FEDGE0, FEDGE1, FEDGE2, FEDGE3]: nInds = E[C[activeCells,FEDGE],:][:,[ENODE0,ENODE1]] - eX = np.c_[N[nInds[:,0],1], N[nInds[:,1],1], [np.nan]*nInds.shape[0]] - eY = np.c_[N[nInds[:,0],2], N[nInds[:,1],2], [np.nan]*nInds.shape[0]] + eX = np.c_[N[nInds[:,0],NX], N[nInds[:,1],NX], [np.nan]*nInds.shape[0]] + eY = np.c_[N[nInds[:,0],NY], N[nInds[:,1],NY], [np.nan]*nInds.shape[0]] plt.plot(eX.flatten(), eY.flatten(), 'b-') gridCC = self.gridCC if text: [ax.text(cc[0], cc[1],i) for i, cc in enumerate(gridCC)] plt.plot(gridCC[:,0], gridCC[:,1], 'r.') - gridFx = self.gridEy - gridFy = self.gridEx + gridFx = self.gridFx + gridFy = self.gridFy if text: [ax.text(cc[0], cc[1],i) for i, cc in enumerate(np.vstack((gridFx,gridFy)))] gridEx = self.gridEx @@ -392,8 +421,8 @@ class TreeMesh(object): # for E in self._edges: # if E[ACTIVE] == 0: continue - # ex = N[E[[ENODE0,ENODE1]],1] - # ey = N[E[[ENODE0,ENODE1]],2] + # ex = N[E[[ENODE0,ENODE1]],NX] + # ey = N[E[[ENODE0,ENODE1]],NY] # ax.plot(ex, ey, 'b-') # ax.text(ex.mean(), ey.mean(), E[NUM]) @@ -417,13 +446,18 @@ if __name__ == '__main__': # tM.number() # print tM._index('_edges',3)[1] - print tM.vol # print tM._edges[:,[0,1,3, 4,5 ]] plt.subplot(211) plt.spy(tM.faceDiv) tM.plotGrid(ax=plt.subplot(212)) + + + print tM.vol + print tM.area + print tM.edge + # plt.figure(2) # plt.plot(SortByX0(tM.gridCC),'b.') plt.show()