wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 04:02:43 +08:00
Code Issues Packages Projects Releases Wiki Activity

FaceInnerProducts for TreeMesh

Browse Source
This commit is contained in:
rowanc1
2014-02-10 22:35:28 -08:00
parent 61e86c2dd9
commit a808aa5774
3 changed files with 58 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/InnerProducts.py
+2 -2
View File
@@ -101,7 +101,7 @@ class InnerProducts(object):
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P000, P001, P010, P011, P100, P101, P110, P111]
P = [P000, P100, P010, P110, P001, P101, P011, P111]
Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
@@ -205,7 +205,7 @@ class InnerProducts(object):
If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
P = [P000, P001, P010, P011, P100, P101, P110, P111]
P = [P000, P100, P010, P110, P001, P101, P011, P111]
Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
SimPEG/Mesh/TreeMesh.py
+51 -7
View File
@@ -1,4 +1,6 @@
from SimPEG import np, sp, Utils, Solver
from BaseMesh import BaseMesh
from InnerProducts import InnerProducts
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
@@ -582,15 +584,16 @@ class TreeCell(TreeObject):
self.mesh.cells.remove(self)
@property
def faceIndex(self):
def faceIndex(self, theseFaces='all', addDirection=True):
#TODO: preallocate
I, J, V = np.empty(0,dtype=float), np.empty(0,dtype=float), np.empty(0,dtype=float)
for face in self.faces:
thisFace = 'all' == theseFaces or face in theseFaces
if not thisFace: continue
j = self.faces[face].index
i = j*0+self.num
v = j*0+1
if 'm' in face:
if addDirection and 'm' in face:
v *= -1
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
return I, J, V
@@ -598,7 +601,6 @@ class TreeCell(TreeObject):
@property
def vol(self): return self.sz.prod()
def viz(self, ax, color='none', text=False):
if not self.isleaf: return
x0, sz = self.x0, self.sz
@@ -615,8 +617,11 @@ class TreeCell(TreeObject):
if text: ax.text(self.center[0], self.center[1], self.center[2], self.num)
class TreeMesh(object):
class TreeMesh(InnerProducts, BaseMesh):
"""TreeMesh"""
_meshType = 'TREE'
def __init__(self, h_in, x0=None):
assert type(h_in) is list, 'h_in must be a list'
h = range(len(h_in))
@@ -635,7 +640,9 @@ class TreeMesh(object):
assert type(x0) in [list, np.ndarray], 'x0 must be a numpy array or a list'
x0 = np.array(x0, dtype=float)
assert len(x0) == self.dim, 'x0 must have the same dimensions as the mesh'
self.x0 = x0
# TODO: this has a lot of stuff which doesn't work for this style of mesh...
BaseMesh.__init__(self, np.array([x.size for x in h]), x0)
# set the sets for holding the cells, nodes, faces, and edges
self.cells = set()
@@ -846,7 +853,7 @@ class TreeMesh(object):
# TODO: Preallocate!
I, J, V = np.empty(0), np.empty(0), np.empty(0)
for cell in self.sortedCells:
i, j, v = cell.faceIndex
i, j, v = cell.faceIndex('all')
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
VOL = self.vol
@@ -855,6 +862,43 @@ class TreeMesh(object):
self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S)
return self._faceDiv
def _getFacePxx(self, xFace, yFace):
self.number()
I, J, V = np.empty(0), np.empty(0), np.empty(0)
for cell in self.sortedCells:
i, j, v = cell.faceIndex('fX'+xFace)
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
xP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF))
I, J, V = np.empty(0), np.empty(0), np.empty(0)
for cell in self.sortedCells:
i, j, v = cell.faceIndex('fY'+yFace)
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
yP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF))
return sp.vstack((xP, yP))
def _getFacePxxx(self, xFace, yFace, zFace):
self.number()
I, J, V = np.empty(0), np.empty(0), np.empty(0)
for cell in self.sortedCells:
i, j, v = cell.faceIndex('fX'+xFace)
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
xP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF))
I, J, V = np.empty(0), np.empty(0), np.empty(0)
for cell in self.sortedCells:
i, j, v = cell.faceIndex('fY'+yFace)
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
yP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF))
I, J, V = np.empty(0), np.empty(0), np.empty(0)
for cell in self.sortedCells:
i, j, v = cell.faceIndex('fZ'+zFace)
I, J, V = np.r_[I,i], np.r_[J,j], np.r_[V,v]
zP = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF))
return sp.vstack((xP, yP, zP))
def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, plotE=False, plotEx=False, plotEy=False, plotEz=False, showIt=False):
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
SimPEG/Tests/test_TreeMesh.py
+5 -3
View File
@@ -489,13 +489,15 @@ class SimpleOctreeOperatorTests(unittest.TestCase):
h1 = np.random.rand(5)
h2 = np.random.rand(7)
h3 = np.random.rand(3)
self.tM = TensorMesh([h1,h2,1])
self.oM = TreeMesh([h1,h2,1])
self.tM = TensorMesh([h1,h2,h3])
self.oM = TreeMesh([h1,h2,h3])
def test_faceDiv(self):
print (self.tM.faceDiv - self.oM.faceDiv)
self.assertTrue((self.tM.faceDiv - self.oM.faceDiv).toarray().sum() == 0)
def test_InnerProducts(self):
self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() == 0)
if __name__ == '__main__':
unittest.main()