wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-12 12:30:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

updates to edges and areas

Browse Source
This commit is contained in:
Rowan Cockett
2015-02-09 21:30:54 -08:00
parent a7a6671ad4
commit 310f327b5a
1 changed files with 77 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/NewTreeMesh.py
+77 -43
View File
@@ -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()