mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-01 08:37:00 +08:00
LRM --> Curv
This commit is contained in:
@@ -27,7 +27,7 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
|
||||
|
||||
__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 CurvilinearMesh(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)
|
||||
|
||||
|
||||
@@ -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]))
|
||||
|
||||
@@ -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:
|
||||
|
||||
@@ -4,7 +4,7 @@ from SimPEG.Mesh import TensorMesh, CurvilinearMesh
|
||||
from SimPEG.Utils import ndgrid
|
||||
|
||||
|
||||
class BasicLRMTests(unittest.TestCase):
|
||||
class BasicCurvTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
a = np.array([1, 1, 1])
|
||||
@@ -13,91 +13,91 @@ class BasicLRMTests(unittest.TestCase):
|
||||
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 = CurvilinearMesh([X, Y])
|
||||
self.Curv2 = CurvilinearMesh([X, Y])
|
||||
X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
|
||||
self.TM3 = TensorMesh([a, b, c])
|
||||
self.LRM3 = CurvilinearMesh([X, Y, Z])
|
||||
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.LRM3.area == test_area))
|
||||
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.LRM3.vol, test_vol)
|
||||
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.LRM2.vol == test_vol)
|
||||
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.LRM3.edge == test_edge)
|
||||
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.LRM2.edge == test_edge)
|
||||
t1 = np.all(self.Curv2.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.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.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)))
|
||||
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.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.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.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)))
|
||||
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.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.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.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)))
|
||||
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.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.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.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)))
|
||||
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.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.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.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))
|
||||
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__':
|
||||
|
||||
@@ -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]
|
||||
|
||||
|
||||
@@ -7,7 +7,7 @@ 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.CurvilinearMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
@@ -24,7 +24,7 @@ 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.CurvilinearMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
@@ -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'))
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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)
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 30 KiB |
@@ -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::
|
||||
|
||||
|
||||
+1
-1
@@ -23,7 +23,7 @@ Solver Utilities
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
LRM Utilities
|
||||
Curv Utilities
|
||||
=============
|
||||
|
||||
.. automodule:: SimPEG.Utils.lrmutils
|
||||
|
||||
Reference in New Issue
Block a user