wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-29 13:50:06 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #96 from simpeg/Curvilinear

Browse Source 
Changed LogicallyRectMesh to CurvilinearMesh
This commit is contained in:
Lindsey
2015-05-15 08:44:28 -07:00
parent 8a36cbab3b e9fea3bad2
commit 1a2edfa8f1
16 changed files with 198 additions and 198 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/LogicallyRectMesh.py → SimPEG/Mesh/CurvilinearMesh.py
+8 -8
View File
@@ -10,9 +10,9 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
"""
LogicallyRectMesh is a mesh class that deals with logically rectangular meshes.
CurvilinearMesh is a mesh class that deals with logically rectangular meshes.
Example of a logically rectangular mesh:
@@ -21,13 +21,13 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
from SimPEG import Mesh, Utils
X, Y = Utils.exampleLrmGrid([3,3],'rotate')
M = Mesh.LogicallyRectMesh([X, Y])
M = Mesh.CurvilinearMesh([X, Y])
M.plotGrid(showIt=True)
"""
__metaclass__ = Utils.SimPEGMetaClass
_meshType = 'LRM'
_meshType = 'Curv'
def __init__(self, nodes):
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
@@ -38,7 +38,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
assert len(nodes[0].shape) > 1, "Not worth using LRM for a 1D mesh."
assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh."
BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None)
@@ -343,7 +343,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
from SimPEG import Mesh, Utils
X, Y = Utils.exampleLrmGrid([3,3],'rotate')
M = Mesh.LogicallyRectMesh([X, Y])
M = Mesh.CurvilinearMesh([X, Y])
M.plotGrid(showIt=True)
"""
@@ -435,9 +435,9 @@ if __name__ == '__main__':
dee3 = True
if dee3:
X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False)
M = LogicallyRectMesh([X, Y, Z])
M = CurvilinearMesh([X, Y, Z])
else:
X, Y = Utils.ndgrid(h1, h2, vector=False)
M = LogicallyRectMesh([X, Y])
M = CurvilinearMesh([X, Y])
print M.r(M.normals, 'F', 'Fx', 'V')
SimPEG/Mesh/InnerProducts.py
+8 -8
View File
@@ -328,7 +328,7 @@ class InnerProducts(object):
iijj = ndgrid(i, j)
ii, jj = iijj[:, 0], iijj[:, 1]
if M._meshType == 'LRM':
if M._meshType == 'Curv':
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
@@ -353,7 +353,7 @@ class InnerProducts(object):
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
if M._meshType == 'LRM':
if M._meshType == 'Curv':
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
PXX = I2x2 * PXX
@@ -376,7 +376,7 @@ class InnerProducts(object):
iijjkk = ndgrid(i, j, k)
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
if M._meshType == 'LRM':
if M._meshType == 'Curv':
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
fN3 = M.r(M.normals, 'F', 'Fz', 'M')
@@ -410,7 +410,7 @@ class InnerProducts(object):
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
if M._meshType == 'LRM':
if M._meshType == 'Curv':
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
@@ -432,7 +432,7 @@ class InnerProducts(object):
iijj = ndgrid(i, j)
ii, jj = iijj[:, 0], iijj[:, 1]
if M._meshType == 'LRM':
if M._meshType == 'Curv':
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
@@ -452,7 +452,7 @@ class InnerProducts(object):
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
if M._meshType == 'LRM':
if M._meshType == 'Curv':
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
PXX = I2x2 * PXX
@@ -466,7 +466,7 @@ class InnerProducts(object):
iijjkk = ndgrid(i, j, k)
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
if M._meshType == 'LRM':
if M._meshType == 'Curv':
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
@@ -495,7 +495,7 @@ class InnerProducts(object):
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
if M._meshType == 'LRM':
if M._meshType == 'Curv':
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
SimPEG/Mesh/__init__.py
+1 -1
View File
@@ -1,5 +1,5 @@
from TensorMesh import TensorMesh
from CylMesh import CylMesh
from LogicallyRectMesh import LogicallyRectMesh
from CurvilinearMesh import CurvilinearMesh
from TreeMesh import TreeMesh
from BaseMesh import BaseMesh
SimPEG/Tests/TestUtils.py
+4 -4
View File
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
from numpy.linalg import norm
from SimPEG.Utils import mkvc, sdiag, diagEst
from SimPEG import Utils
from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh
from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
import numpy as np
import scipy.sparse as sp
import unittest
@@ -115,7 +115,7 @@ class OrderTest(unittest.TestCase):
max_h = max([np.max(hi) for hi in self.M.h])
return max_h
elif 'LRM' in self._meshType:
elif 'Curv' in self._meshType:
if 'uniform' in self._meshType:
kwrd = 'rect'
elif 'rotate' in self._meshType:
@@ -126,10 +126,10 @@ class OrderTest(unittest.TestCase):
raise Exception('Lom not supported for 1D')
elif self.meshDimension == 2:
X, Y = Utils.exampleLrmGrid([nc, nc], kwrd)
self.M = LogicallyRectMesh([X, Y])
self.M = CurvilinearMesh([X, Y])
elif self.meshDimension == 3:
X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd)
self.M = LogicallyRectMesh([X, Y, Z])
self.M = CurvilinearMesh([X, Y, Z])
return 1./nc
def getError(self):
SimPEG/Tests/test_CurvilinearMesh.py
+104
View File
@@ -0,0 +1,104 @@
import numpy as np
import unittest
from SimPEG.Mesh import TensorMesh, CurvilinearMesh
from SimPEG.Utils import ndgrid
class BasicCurvTests(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.Curv2 = CurvilinearMesh([X, Y])
X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
self.TM3 = TensorMesh([a, b, c])
self.Curv3 = CurvilinearMesh([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.Curv3.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.Curv3.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.Curv2.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.Curv3.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.Curv2.edge == test_edge)
self.assertTrue(t1)
def test_tangents(self):
T = self.Curv2.tangents
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx)))
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx)))
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy)))
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy)))
T = self.Curv3.tangents
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz)))
def test_normals(self):
N = self.Curv2.normals
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx)))
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx)))
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy)))
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy)))
N = self.Curv3.normals
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz)))
def test_grid(self):
self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))
self.assertTrue(np.all(self.Curv2.gridN == self.TM2.gridN))
self.assertTrue(np.all(self.Curv2.gridFx == self.TM2.gridFx))
self.assertTrue(np.all(self.Curv2.gridFy == self.TM2.gridFy))
self.assertTrue(np.all(self.Curv2.gridEx == self.TM2.gridEx))
self.assertTrue(np.all(self.Curv2.gridEy == self.TM2.gridEy))
self.assertTrue(np.all(self.Curv3.gridCC == self.TM3.gridCC))
self.assertTrue(np.all(self.Curv3.gridN == self.TM3.gridN))
self.assertTrue(np.all(self.Curv3.gridFx == self.TM3.gridFx))
self.assertTrue(np.all(self.Curv3.gridFy == self.TM3.gridFy))
self.assertTrue(np.all(self.Curv3.gridFz == self.TM3.gridFz))
self.assertTrue(np.all(self.Curv3.gridEx == self.TM3.gridEx))
self.assertTrue(np.all(self.Curv3.gridEy == self.TM3.gridEy))
self.assertTrue(np.all(self.Curv3.gridEz == self.TM3.gridEz))
if __name__ == '__main__':
unittest.main()
SimPEG/Tests/test_LogicallyRectMesh.py
-104
View File
@@ -1,104 +0,0 @@
import numpy as np
import unittest
from SimPEG.Mesh import TensorMesh, LogicallyRectMesh
from SimPEG.Utils import ndgrid
class BasicLRMTests(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.LRM2 = LogicallyRectMesh([X, Y])
X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
self.TM3 = TensorMesh([a, b, c])
self.LRM3 = LogicallyRectMesh([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.LRM3.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.LRM3.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.LRM2.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.LRM3.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.LRM2.edge == test_edge)
self.assertTrue(t1)
def test_tangents(self):
T = self.LRM2.tangents
self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM2.nEx)))
self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM2.nEx)))
self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM2.nEy)))
self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM2.nEy)))
T = self.LRM3.tangents
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM3.nEx)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM3.nEx)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LRM3.nEx)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM3.nEy)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM3.nEy)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LRM3.nEy)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LRM3.nEz)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LRM3.nEz)))
self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LRM3.nEz)))
def test_normals(self):
N = self.LRM2.normals
self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM2.nFx)))
self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM2.nFx)))
self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM2.nFy)))
self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM2.nFy)))
N = self.LRM3.normals
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM3.nFx)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM3.nFx)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LRM3.nFx)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM3.nFy)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM3.nFy)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LRM3.nFy)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LRM3.nFz)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LRM3.nFz)))
self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LRM3.nFz)))
def test_grid(self):
self.assertTrue(np.all(self.LRM2.gridCC == self.TM2.gridCC))
self.assertTrue(np.all(self.LRM2.gridN == self.TM2.gridN))
self.assertTrue(np.all(self.LRM2.gridFx == self.TM2.gridFx))
self.assertTrue(np.all(self.LRM2.gridFy == self.TM2.gridFy))
self.assertTrue(np.all(self.LRM2.gridEx == self.TM2.gridEx))
self.assertTrue(np.all(self.LRM2.gridEy == self.TM2.gridEy))
self.assertTrue(np.all(self.LRM3.gridCC == self.TM3.gridCC))
self.assertTrue(np.all(self.LRM3.gridN == self.TM3.gridN))
self.assertTrue(np.all(self.LRM3.gridFx == self.TM3.gridFx))
self.assertTrue(np.all(self.LRM3.gridFy == self.TM3.gridFy))
self.assertTrue(np.all(self.LRM3.gridFz == self.TM3.gridFz))
self.assertTrue(np.all(self.LRM3.gridEx == self.TM3.gridEx))
self.assertTrue(np.all(self.LRM3.gridEy == self.TM3.gridEy))
self.assertTrue(np.all(self.LRM3.gridEz == self.TM3.gridEz))
if __name__ == '__main__':
unittest.main()
SimPEG/Tests/test_innerProduct.py
+2 -2
View File
@@ -7,7 +7,7 @@ from SimPEG import Utils
class TestInnerProducts(OrderTest):
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
meshDimension = 3
meshSizes = [16, 32]
@@ -154,7 +154,7 @@ class TestInnerProducts(OrderTest):
class TestInnerProducts2D(OrderTest):
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
meshDimension = 2
meshSizes = [4, 8, 16, 32, 64, 128]
SimPEG/Tests/test_innerProductDerivs.py
+60 -60
View File
@@ -7,9 +7,9 @@ from TestUtils import checkDerivative
class TestInnerProductsDerivs(unittest.TestCase):
def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False):
if meshType == 'LRM':
if meshType == 'Curv':
hRect = Utils.exampleLrmGrid(h,'rotate')
mesh = Mesh.LogicallyRectMesh(hRect)
mesh = Mesh.CurvilinearMesh(hRect)
elif meshType == 'Tree':
mesh = Mesh.TreeMesh(h)
elif meshType == 'Tensor':
@@ -24,9 +24,9 @@ class TestInnerProductsDerivs(unittest.TestCase):
return checkDerivative(fun, sig, num=5, plotIt=False)
def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
if meshType == 'LRM':
if meshType == 'Curv':
hRect = Utils.exampleLrmGrid(h,'rotate')
mesh = Mesh.LogicallyRectMesh(hRect)
mesh = Mesh.CurvilinearMesh(hRect)
elif meshType == 'Tree':
mesh = Mesh.TreeMesh(h)
elif meshType == 'Tensor':
@@ -137,65 +137,65 @@ class TestInnerProductsDerivs(unittest.TestCase):
def test_FaceIP_2D_float_LRM(self):
self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM'))
def test_FaceIP_3D_float_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM'))
def test_FaceIP_2D_isotropic_LRM(self):
self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM'))
def test_FaceIP_3D_isotropic_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM'))
def test_FaceIP_2D_anisotropic_LRM(self):
self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM'))
def test_FaceIP_3D_anisotropic_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM'))
def test_FaceIP_2D_tensor_LRM(self):
self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM'))
def test_FaceIP_3D_tensor_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM'))
def test_FaceIP_2D_float_Curv(self):
self.assertTrue(self.doTestFace([10, 4],0, False, 'Curv'))
def test_FaceIP_3D_float_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Curv'))
def test_FaceIP_2D_isotropic_Curv(self):
self.assertTrue(self.doTestFace([10, 4],1, False, 'Curv'))
def test_FaceIP_3D_isotropic_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Curv'))
def test_FaceIP_2D_anisotropic_Curv(self):
self.assertTrue(self.doTestFace([10, 4],2, False, 'Curv'))
def test_FaceIP_3D_anisotropic_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Curv'))
def test_FaceIP_2D_tensor_Curv(self):
self.assertTrue(self.doTestFace([10, 4],3, False, 'Curv'))
def test_FaceIP_3D_tensor_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Curv'))
def test_FaceIP_2D_float_fast_LRM(self):
self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM'))
def test_FaceIP_3D_float_fast_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM'))
def test_FaceIP_2D_isotropic_fast_LRM(self):
self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM'))
def test_FaceIP_3D_isotropic_fast_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM'))
def test_FaceIP_2D_anisotropic_fast_LRM(self):
self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM'))
def test_FaceIP_3D_anisotropic_fast_LRM(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM'))
def test_FaceIP_2D_float_fast_Curv(self):
self.assertTrue(self.doTestFace([10, 4],0, True, 'Curv'))
def test_FaceIP_3D_float_fast_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Curv'))
def test_FaceIP_2D_isotropic_fast_Curv(self):
self.assertTrue(self.doTestFace([10, 4],1, True, 'Curv'))
def test_FaceIP_3D_isotropic_fast_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Curv'))
def test_FaceIP_2D_anisotropic_fast_Curv(self):
self.assertTrue(self.doTestFace([10, 4],2, True, 'Curv'))
def test_FaceIP_3D_anisotropic_fast_Curv(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Curv'))
def test_EdgeIP_2D_float_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM'))
def test_EdgeIP_3D_float_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM'))
def test_EdgeIP_2D_isotropic_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM'))
def test_EdgeIP_3D_isotropic_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM'))
def test_EdgeIP_2D_anisotropic_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM'))
def test_EdgeIP_3D_anisotropic_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM'))
def test_EdgeIP_2D_tensor_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM'))
def test_EdgeIP_3D_tensor_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM'))
def test_EdgeIP_2D_float_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],0, False, 'Curv'))
def test_EdgeIP_3D_float_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Curv'))
def test_EdgeIP_2D_isotropic_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],1, False, 'Curv'))
def test_EdgeIP_3D_isotropic_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Curv'))
def test_EdgeIP_2D_anisotropic_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],2, False, 'Curv'))
def test_EdgeIP_3D_anisotropic_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Curv'))
def test_EdgeIP_2D_tensor_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],3, False, 'Curv'))
def test_EdgeIP_3D_tensor_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Curv'))
def test_EdgeIP_2D_float_fast_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM'))
def test_EdgeIP_3D_float_fast_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM'))
def test_EdgeIP_2D_isotropic_fast_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM'))
def test_EdgeIP_3D_isotropic_fast_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM'))
def test_EdgeIP_2D_anisotropic_fast_LRM(self):
self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM'))
def test_EdgeIP_3D_anisotropic_fast_LRM(self):
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM'))
def test_EdgeIP_2D_float_fast_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],0, True, 'Curv'))
def test_EdgeIP_3D_float_fast_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Curv'))
def test_EdgeIP_2D_isotropic_fast_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],1, True, 'Curv'))
def test_EdgeIP_3D_isotropic_fast_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Curv'))
def test_EdgeIP_2D_anisotropic_fast_Curv(self):
self.assertTrue(self.doTestEdge([10, 4],2, True, 'Curv'))
def test_EdgeIP_3D_anisotropic_fast_Curv(self):
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Curv'))
SimPEG/Tests/test_operators.py
+3 -3
View File
@@ -4,7 +4,7 @@ from TestUtils import OrderTest
import matplotlib.pyplot as plt
#TODO: 'randomTensorMesh'
MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
MESHTYPES = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
@@ -38,7 +38,7 @@ class TestCurl(OrderTest):
curlE_ana = self.M.projectFaceVector(Fc)
curlE = self.M.edgeCurl.dot(E)
if self._meshType == 'rotateLRM':
if self._meshType == 'rotateCurv':
# Really it is the integration we should be caring about:
# So, let us look at the l2 norm.
err = np.linalg.norm(self.M.area*(curlE - curlE_ana), 2)
@@ -229,7 +229,7 @@ class TestFaceDiv3D(OrderTest):
divF = self.M.faceDiv.dot(F)
divF_ana = call3(sol, self.M.gridCC)
if self._meshType == 'rotateLRM':
if self._meshType == 'rotateCurv':
# Really it is the integration we should be caring about:
# So, let us look at the l2 norm.
err = np.linalg.norm(self.M.vol*(divF-divF_ana), 2)
SimPEG/Utils/__init__.py
+1 -1
View File
@@ -1,7 +1,7 @@
from matutils import *
from codeutils import *
from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
from lrmutils import volTetra, faceInfo, indexCube
from curvutils import volTetra, faceInfo, indexCube
from interputils import interpmat
from ipythonutils import easyAnimate as animate
from CounterUtils import *
SimPEG/Utils/lrmutils.py → SimPEG/Utils/curvutils.py
View File
docs/SimPEGFrameworkRevised.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 30 KiB

docs/api_InnerProducts.rst
+2 -2
View File
@@ -75,7 +75,7 @@ We multiply by square-root of volume on each side of the tensor conductivity to
.. math::
\mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
\mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM}
\mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv}
Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above.
We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
@@ -114,7 +114,7 @@ Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of
.. math::
\mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)}
\mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{Curv only}} \mathbf{Q}_{(i)}
.. note::
docs/api_Mesh.rst
+2 -2
View File
@@ -29,7 +29,7 @@ the implementations.
tM = Mesh.TensorMesh(sz)
qM = Mesh.TreeMesh(sz)
qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0)
rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
fig, axes = plt.subplots(1,3,figsize=(14,5))
opts = {}
@@ -38,7 +38,7 @@ the implementations.
qM.plotGrid(ax=axes[1], **opts)
axes[1].set_title('TreeMesh')
rM.plotGrid(ax=axes[2], **opts)
axes[2].set_title('LogicallyRectMesh')
axes[2].set_title('CurvilinearMesh')
plt.show()
docs/api_MeshCode.rst
+1 -1
View File
@@ -21,7 +21,7 @@ Tree Mesh
Logically Rectangular Mesh
==========================
.. automodule:: SimPEG.Mesh.LogicallyRectMesh
.. automodule:: SimPEG.Mesh.Curvilinear
:show-inheritance:
:members:
:undoc-members:
docs/api_Utils.rst
+2 -2
View File
@@ -23,10 +23,10 @@ Solver Utilities
:members:
:undoc-members:
LRM Utilities
Curv Utilities
=============
.. automodule:: SimPEG.Utils.lrmutils
.. automodule:: SimPEG.Utils.curvutils
:members:
:undoc-members: