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

edge and tangent calculating, along with testing functions.

Browse Source
This commit is contained in:
Rowan Cockett
2013-07-31 15:07:14 -07:00
parent b07256bbc5
commit 34e380ccb9
3 changed files with 104 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/LogicallyOrthogonalMesh.py
+51 -6
View File
@@ -3,12 +3,20 @@ from BaseMesh import BaseMesh
from DiffOperators import DiffOperators
from utils import mkvc, ndgrid, volTetra, indexCube, faceInfo
# Some helper functions.
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
normalize2D = lambda x: x/np.kron(np.ones((1, 2)), mkvc(length2D(x), 2))
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), mkvc(length3D(x), 2))
class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid
"""
LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes.
"""
_meshType = 'LOM'
def __init__(self, nodes):
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
@@ -117,12 +125,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid
# Compute areas of cell faces
if(self.dim == 2):
xy = self.gridN
length = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
area1 = length(xy[B, :] - xy[A, :])
area1 = length2D(xy[B, :] - xy[A, :])
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
area2 = length(xy[D, :] - xy[A, :])
area2 = length2D(xy[D, :] - xy[A, :])
self._area = np.r_[mkvc(area1), mkvc(area2)]
elif(self.dim == 3):
@@ -141,6 +147,45 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid
_area = None
area = property(**area())
def edge():
doc = "Edge legnths."
def fget(self):
if(self._edge is None or self._tangents is None):
if(self.dim == 2):
xy = self.gridN
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
edge1 = xy[D, :] - xy[A, :]
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
edge2 = xy[B, :] - xy[A, :]
self._edge = np.r_[mkvc(length2D(edge1)), mkvc(length2D(edge2))]
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
elif(self.dim == 3):
xyz = self.gridN
A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
edge1 = xyz[D, :] - xyz[A, :]
A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
edge2 = xyz[B, :] - xyz[A, :]
A, E = indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
edge3 = xyz[E, :] - xyz[A, :]
self._edge = np.r_[mkvc(length3D(edge1)), mkvc(length3D(edge2)), mkvc(length3D(edge3))]
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
return self._edge
return locals()
_edge = None
edge = property(**edge())
def tangents():
doc = "Edge tangents."
def fget(self):
if(self._tangents is None):
self.edge # calling .edge will create the tangents
return self._tangents
return locals()
_tangents = None
tangents = property(**tangents())
if __name__ == '__main__':
nc = 5
@@ -148,7 +193,7 @@ if __name__ == '__main__':
nc = 7
h2 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
dee3 = True
dee3 = False
if dee3:
X, Y, Z = ndgrid(h1, h2, h3, vector=False)
M = LogicallyOrthogonalMesh([X, Y, Z])
@@ -159,4 +204,4 @@ if __name__ == '__main__':
# print M.r(M.gridCC, format='M')
# print M.gridN[:, 0]
print M.nE
print M.area
print M.r(M.tangents, 'E', 'Ex', 'M')
SimPEG/tests/test_LogicallyOrthogonalMesh.py
+52
View File
@@ -0,0 +1,52 @@
import numpy as np
import unittest
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
from OrderTest import OrderTest
from scipy.sparse.linalg import dsolve
from utils import ndgrid
class BasicLOMTests(unittest.TestCase):
def setUp(self):
a = np.array([1, 1, 1])
b = np.array([1, 2])
c = np.array([1, 4])
gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
X, Y = ndgrid(gridIt([a, b]), vector=False)
self.TM2 = TensorMesh([a, b])
self.LOM2 = LogicallyOrthogonalMesh([X, Y])
X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
self.TM3 = TensorMesh([a, b, c])
self.LOM3 = LogicallyOrthogonalMesh([X, Y, Z])
def test_area_3D(self):
test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
self.assertTrue(np.all(self.LOM3.area == test_area))
def test_vol_3D(self):
test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
np.testing.assert_almost_equal(self.LOM3.vol, test_vol)
self.assertTrue(True) # Pass if you get past the assertion.
def test_vol_2D(self):
test_vol = np.array([1, 1, 1, 2, 2, 2])
t1 = np.all(self.LOM2.vol == test_vol)
self.assertTrue(t1)
def test_edge_3D(self):
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
t1 = np.all(self.LOM3.edge == test_edge)
self.assertTrue(t1)
def test_edge_2D(self):
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
t1 = np.all(self.LOM2.edge == test_edge)
self.assertTrue(t1)
if __name__ == '__main__':
unittest.main()
SimPEG/utils.py
+1 -1
View File
@@ -259,7 +259,7 @@ def faceInfo(xyz, A, B, C, D, average=True):
nD = cross(DA, CD)
length = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(N), 2))
normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(x), 2))
if average:
# average the normals at each vertex.
N = (nA + nB + nC + nD)/4 # this is intrinsically weighted by area