mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-11 11:26:01 +08:00
Separate tests into folders.
Build in a matrix?
This commit is contained in:
Rowan Cockett
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -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()
|
||||
@@ -0,0 +1,505 @@
|
||||
from SimPEG.Mesh import TensorMesh
|
||||
from SimPEG.Mesh.TreeMesh import TreeMesh, TreeFace, TreeCell
|
||||
import numpy as np
|
||||
import unittest
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
TOL = 1e-10
|
||||
|
||||
class TestOcTreeObjects(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.M = TreeMesh([2,1,1])
|
||||
self.M.number()
|
||||
|
||||
self.Mr = TreeMesh([2,1,1])
|
||||
self.Mr.children[0,0,0].refine()
|
||||
self.Mr.number()
|
||||
|
||||
def q(s):
|
||||
if s[0] == 'M':
|
||||
m = self.M
|
||||
s = s[1:]
|
||||
else:
|
||||
m = self.Mr
|
||||
c = m.sortedCells[int(s[1])]
|
||||
if len(s) == 2: return c
|
||||
if s[2] == 'f' and len(s) == 5: return c.faceDict[s[2:]]
|
||||
if s[2] == 'f': return getattr(c.faceDict[s[2:5]], 'edg' +s[5:])
|
||||
if s[2] == 'e': return getattr(c,s[2:])
|
||||
if s[2] == 'n': return getattr(c,'node'+s[3:])
|
||||
|
||||
self.q = q
|
||||
|
||||
def test_counts(self):
|
||||
self.assertTrue(self.M.nC == 2)
|
||||
self.assertTrue(self.M.nFx == 3)
|
||||
self.assertTrue(self.M.nFy == 4)
|
||||
self.assertTrue(self.M.nFz == 4)
|
||||
self.assertTrue(self.M.nF == 11)
|
||||
self.assertTrue(self.M.nEx == 8)
|
||||
self.assertTrue(self.M.nEy == 6)
|
||||
self.assertTrue(self.M.nEz == 6)
|
||||
self.assertTrue(self.M.nE == 20)
|
||||
self.assertTrue(self.M.nN == 12)
|
||||
|
||||
self.assertTrue(self.Mr.nC == 9)
|
||||
self.assertTrue(self.Mr.nFx == 13)
|
||||
self.assertTrue(self.Mr.nFy == 14)
|
||||
self.assertTrue(self.Mr.nFz == 14)
|
||||
self.assertTrue(self.Mr.nF == 41)
|
||||
|
||||
|
||||
for cell in self.Mr.sortedCells:
|
||||
for e in cell.edgeDict:
|
||||
self.assertTrue(cell.edgeDict[e].edgeType==e[1].lower())
|
||||
|
||||
self.assertTrue(self.Mr.nN == 31)
|
||||
self.assertTrue(self.Mr.nEx == 22)
|
||||
self.assertTrue(self.Mr.nEy == 20)
|
||||
self.assertTrue(self.Mr.nEz == 20)
|
||||
|
||||
def test_sizes(self):
|
||||
q = self.q
|
||||
|
||||
for key in ['Mc0','Mc1']:
|
||||
self.assertTrue(q(key).vol == 0.5)
|
||||
self.assertTrue(q(key+'fXm').area == 1.)
|
||||
self.assertTrue(q(key+'fXp').area == 1.)
|
||||
self.assertTrue(q(key+'fYm').area == 0.5)
|
||||
self.assertTrue(q(key+'fYp').area == 0.5)
|
||||
self.assertTrue(q(key+'fZm').area == 0.5)
|
||||
self.assertTrue(q(key+'fZp').area == 0.5)
|
||||
|
||||
def test_pointersM(self):
|
||||
q = self.q
|
||||
|
||||
self.assertTrue(q('Mc0fXp') is q('Mc1fXm'))
|
||||
self.assertTrue(q('Mc0fXpe0') is q('Mc1fXme0'))
|
||||
self.assertTrue(q('Mc0fXpe1') is q('Mc1fXme1'))
|
||||
self.assertTrue(q('Mc0fXpe2') is q('Mc1fXme2'))
|
||||
self.assertTrue(q('Mc0fXpe3') is q('Mc1fXme3'))
|
||||
self.assertTrue(q('Mc0fYp') is not q('c1fYm'))
|
||||
self.assertTrue(q('Mc0fXm') is not q('c1fXm'))
|
||||
|
||||
# Test connectivity of shared edges
|
||||
self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe0'))
|
||||
self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe1'))
|
||||
self.assertTrue(q('Mc0fZpe3') is q('Mc1fZpe2'))
|
||||
self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe3'))
|
||||
|
||||
self.assertTrue(q('Mc0fZme3') is not q('c1fZme0'))
|
||||
self.assertTrue(q('Mc0fZme3') is not q('c1fZme1'))
|
||||
self.assertTrue(q('Mc0fZme3') is q('Mc1fZme2'))
|
||||
self.assertTrue(q('Mc0fZme3') is not q('c1fZme3'))
|
||||
|
||||
self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe0'))
|
||||
self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe1'))
|
||||
self.assertTrue(q('Mc0fYpe3') is q('Mc1fYpe2'))
|
||||
self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe3'))
|
||||
|
||||
self.assertTrue(q('Mc0fYme3') is not q('c1fYme0'))
|
||||
self.assertTrue(q('Mc0fYme3') is not q('c1fYme1'))
|
||||
self.assertTrue(q('Mc0fYme3') is q('Mc1fYme2'))
|
||||
self.assertTrue(q('Mc0fYme3') is not q('c1fYme3'))
|
||||
|
||||
self.assertTrue(q('Mc0fZme3') is q('Mc1fXme0'))
|
||||
self.assertTrue(q('Mc0fZpe3') is q('Mc1fXme1'))
|
||||
self.assertTrue(q('Mc0fYme3') is q('Mc1fXme2'))
|
||||
self.assertTrue(q('Mc0fYpe3') is q('Mc1fXme3'))
|
||||
|
||||
self.assertTrue(q('Mc0fZme3') is q('Mc0fXpe0'))
|
||||
self.assertTrue(q('Mc0fZpe3') is q('Mc0fXpe1'))
|
||||
self.assertTrue(q('Mc0fYme3') is q('Mc0fXpe2'))
|
||||
self.assertTrue(q('Mc0fYpe3') is q('Mc0fXpe3'))
|
||||
|
||||
self.assertTrue(q('Mc1fZme2') is q('Mc1fXme0'))
|
||||
self.assertTrue(q('Mc1fZpe2') is q('Mc1fXme1'))
|
||||
self.assertTrue(q('Mc1fYme2') is q('Mc1fXme2'))
|
||||
self.assertTrue(q('Mc1fYpe2') is q('Mc1fXme3'))
|
||||
|
||||
self.assertTrue(q('Mc1fZme2') is q('Mc0fXpe0'))
|
||||
self.assertTrue(q('Mc1fZpe2') is q('Mc0fXpe1'))
|
||||
self.assertTrue(q('Mc1fYme2') is q('Mc0fXpe2'))
|
||||
self.assertTrue(q('Mc1fYpe2') is q('Mc0fXpe3'))
|
||||
|
||||
|
||||
def test_nodePointers(self):
|
||||
q = self.q
|
||||
c0 = self.Mr.sortedCells[0]
|
||||
c0n0 = c0.node0
|
||||
self.assertTrue(c0n0 is q('c0n0'))
|
||||
self.assertTrue(np.all(q('c0n0').center == np.r_[0,0,0.]))
|
||||
self.assertTrue(q('c0n0').num == 0)
|
||||
self.assertTrue(q('c0n1').num == 1)
|
||||
self.assertTrue(q('c0n2').num == 4)
|
||||
self.assertTrue(q('c0n3').num == 5)
|
||||
self.assertTrue(q('c0n4').num == 11)
|
||||
self.assertTrue(q('c0n5').num == 12)
|
||||
self.assertTrue(q('c0n6').num == 14)
|
||||
self.assertTrue(q('c0n7').num == 15)
|
||||
|
||||
def test_pointersMr(self):
|
||||
q = self.q
|
||||
|
||||
c0 = self.Mr.sortedCells[0]
|
||||
c0fXm = c0.fXm
|
||||
c0eX0 = c0.eX0
|
||||
c0fYme0 = c0.fYm.edge0
|
||||
self.assertTrue(c0 is q('c0'))
|
||||
self.assertTrue(c0fXm is q('c0fXm'))
|
||||
self.assertTrue(c0eX0 is q('c0eX0'))
|
||||
self.assertTrue(c0fYme0 is q('c0fYme0'))
|
||||
|
||||
self.assertTrue(q('c0').depth == 1)
|
||||
self.assertTrue(q('c1').depth == 1)
|
||||
self.assertTrue(q('c2').depth == 0)
|
||||
|
||||
# Make sure we know where the center of the cells are.
|
||||
self.assertTrue(np.all(q('c0').center == np.r_[0.125,0.25,0.25]))
|
||||
self.assertTrue(np.all(q('c1').center == np.r_[0.375,0.25,0.25]))
|
||||
self.assertTrue(np.all(q('c2').center == np.r_[0.75,0.5,0.5]))
|
||||
self.assertTrue(np.all(q('c3').center == np.r_[0.125,0.75,0.25]))
|
||||
self.assertTrue(np.all(q('c4').center == np.r_[0.375,0.75,0.25]))
|
||||
self.assertTrue(np.all(q('c5').center == np.r_[0.125,0.25,0.75]))
|
||||
self.assertTrue(np.all(q('c6').center == np.r_[0.375,0.25,0.75]))
|
||||
self.assertTrue(np.all(q('c7').center == np.r_[0.125,0.75,0.75]))
|
||||
self.assertTrue(np.all(q('c8').center == np.r_[0.375,0.75,0.75]))
|
||||
|
||||
# Test X face connectivity and locations and stuff...
|
||||
self.assertTrue(np.all(q('c0fXm').center == np.r_[0,0.25,0.25]))
|
||||
self.assertTrue(np.all(q('c0fXp').center == np.r_[0.25,0.25,0.25]))
|
||||
self.assertTrue(q('c0fXp') is q('c1fXm'))
|
||||
self.assertTrue(np.all(q('c1fXp').center == np.r_[0.5,0.25,0.25]))
|
||||
self.assertTrue(np.all(q('c2fXm').center == np.r_[0.5,0.5,0.5]))
|
||||
self.assertTrue(q('c2fXm').branchdepth == 1)
|
||||
self.assertTrue(q('c2fXm').children[0,0] is q('c1fXp'))
|
||||
self.assertTrue(np.all(q('c3fXm').center == np.r_[0,0.75,0.25]))
|
||||
self.assertTrue(np.all(q('c3fXp').center == np.r_[0.25,0.75,0.25]))
|
||||
self.assertTrue(q('c4fXm') is q('c3fXp'))
|
||||
self.assertTrue(q('c2fXm').children[1,0] is q('c4fXp'))
|
||||
|
||||
#Test some internal stuff (edges held by cell should be same as inside)
|
||||
for key in ['Mc0', 'Mc1'] + ['c%d'%i for i in range(9)]:
|
||||
self.assertTrue(q(key+'eX0') is q(key+'fZme0'))
|
||||
self.assertTrue(q(key+'eX1') is q(key+'fZme1'))
|
||||
self.assertTrue(q(key+'eX2') is q(key+'fZpe0'))
|
||||
self.assertTrue(q(key+'eX3') is q(key+'fZpe1'))
|
||||
|
||||
self.assertTrue(q(key+'eX0') is q(key+'fYme0'))
|
||||
self.assertTrue(q(key+'eX1') is q(key+'fYpe0'))
|
||||
self.assertTrue(q(key+'eX2') is q(key+'fYme1'))
|
||||
self.assertTrue(q(key+'eX3') is q(key+'fYpe1'))
|
||||
|
||||
self.assertTrue(q(key+'eY0') is q(key+'fXme0'))
|
||||
self.assertTrue(q(key+'eY1') is q(key+'fXpe0'))
|
||||
self.assertTrue(q(key+'eY2') is q(key+'fXme1'))
|
||||
self.assertTrue(q(key+'eY3') is q(key+'fXpe1'))
|
||||
|
||||
self.assertTrue(q(key+'eY0') is q(key+'fZme2'))
|
||||
self.assertTrue(q(key+'eY1') is q(key+'fZme3'))
|
||||
self.assertTrue(q(key+'eY2') is q(key+'fZpe2'))
|
||||
self.assertTrue(q(key+'eY3') is q(key+'fZpe3'))
|
||||
|
||||
self.assertTrue(q(key+'eZ0') is q(key+'fXme2'))
|
||||
self.assertTrue(q(key+'eZ1') is q(key+'fXpe2'))
|
||||
self.assertTrue(q(key+'eZ2') is q(key+'fXme3'))
|
||||
self.assertTrue(q(key+'eZ3') is q(key+'fXpe3'))
|
||||
|
||||
self.assertTrue(q(key+'eZ0') is q(key+'fYme2'))
|
||||
self.assertTrue(q(key+'eZ1') is q(key+'fYme3'))
|
||||
self.assertTrue(q(key+'eZ2') is q(key+'fYpe2'))
|
||||
self.assertTrue(q(key+'eZ3') is q(key+'fYpe3'))
|
||||
|
||||
#Test some edge stuff
|
||||
self.assertTrue(np.all(q('c0eX0').center == np.r_[0.125,0,0]))
|
||||
self.assertTrue(np.all(q('c0eX1').center == np.r_[0.125,0.5,0]))
|
||||
self.assertTrue(np.all(q('c0eX2').center == np.r_[0.125,0,0.5]))
|
||||
self.assertTrue(np.all(q('c0eX3').center == np.r_[0.125,0.5,0.5]))
|
||||
|
||||
self.assertTrue(np.all(q('c5eX0').center == np.r_[0.125,0,0.5]))
|
||||
self.assertTrue(np.all(q('c5eX1').center == np.r_[0.125,0.5,0.5]))
|
||||
self.assertTrue(q('c5eX0') is q('c0eX2'))
|
||||
self.assertTrue(q('c5eX1') is q('c0eX3'))
|
||||
|
||||
self.assertTrue(np.all(q('c0eY0').center == np.r_[0,0.25,0]))
|
||||
self.assertTrue(np.all(q('c0eY1').center == np.r_[0.25,0.25,0]))
|
||||
self.assertTrue(np.all(q('c0eY2').center == np.r_[0,0.25,0.5]))
|
||||
self.assertTrue(np.all(q('c0eY3').center == np.r_[0.25,0.25,0.5]))
|
||||
|
||||
self.assertTrue(np.all(q('c1eY0').center == np.r_[0.25,0.25,0]))
|
||||
self.assertTrue(np.all(q('c1eY2').center == np.r_[0.25,0.25,0.5]))
|
||||
self.assertTrue(q('c1eY0') is q('c0eY1'))
|
||||
self.assertTrue(q('c1eY2') is q('c0eY3'))
|
||||
|
||||
|
||||
self.assertTrue(np.all(q('c0eZ0').center == np.r_[0,0,0.25]))
|
||||
self.assertTrue(np.all(q('c0eZ1').center == np.r_[0.25,0,0.25]))
|
||||
self.assertTrue(np.all(q('c0eZ2').center == np.r_[0,0.5,0.25]))
|
||||
self.assertTrue(np.all(q('c0eZ3').center == np.r_[0.25,0.5,0.25]))
|
||||
|
||||
self.assertTrue(np.all(q('c3eZ0').center == np.r_[0,0.5,0.25]))
|
||||
self.assertTrue(np.all(q('c3eZ1').center == np.r_[0.25,0.5,0.25]))
|
||||
self.assertTrue(q('c3eZ0') is q('c0eZ2'))
|
||||
self.assertTrue(q('c3eZ1') is q('c0eZ3'))
|
||||
|
||||
|
||||
self.assertTrue(q('c0fXp') is q('c1fXm'))
|
||||
self.assertTrue(q('c0fYp') is not q('c1fYm'))
|
||||
self.assertTrue(q('c0fXm') is not q('c1fXm'))
|
||||
|
||||
self.assertTrue(q('c1fXp') is q('c2fXm').children[0,0])
|
||||
|
||||
self.assertTrue(q('c1fYp') is q('c4fYm'))
|
||||
self.assertTrue(q('c1fZp') is q('c6fZm'))
|
||||
|
||||
self.assertTrue(q('c6fXp') is q('c2fXm').children[0,1])
|
||||
|
||||
self.assertTrue(q('c4fXp') is q('c2fXm').children[1,0])
|
||||
|
||||
|
||||
def test_gridCC(self):
|
||||
x = np.r_[0.25,0.75]
|
||||
y = np.r_[0.5,0.5]
|
||||
z = np.r_[0.5,0.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridCC).flatten()) == 0)
|
||||
|
||||
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375]
|
||||
y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75]
|
||||
z = np.r_[0.25,0.25,0.5,0.25,0.25,0.75,0.75,0.75,0.75]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridCC).flatten()) == 0)
|
||||
|
||||
def test_gridN(self):
|
||||
x = np.r_[0,0.5,1,0,0.5,1,0,0.5,1,0,0.5,1]
|
||||
y = np.r_[0,0,0,1,1,1,0,0,0,1,1,1.]
|
||||
z = np.r_[0,0,0,0,0,0,1,1,1,1,1,1.]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridN).flatten()) == 0)
|
||||
|
||||
x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1]
|
||||
y = np.r_[0,0,0,0,0.5,0.5,0.5,1,1,1,1,0,0,0,0.5,0.5,0.5,1,1,1,0,0,0,0,0.5,0.5,0.5,1,1,1,1]
|
||||
z = np.r_[0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1,1,1,1]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridN).flatten()) == 0)
|
||||
|
||||
def test_gridFx(self):
|
||||
x = np.r_[0.0,0.5,1.0]
|
||||
y = np.r_[0.5,0.5,0.5]
|
||||
z = np.r_[0.5,0.5,0.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFx).flatten()) == 0)
|
||||
|
||||
x = np.r_[0.0,0.25,0.5,1.0,0.0,0.25,0.5,0.0,0.25,0.5,0.0,0.25,0.5]
|
||||
y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75]
|
||||
z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.75,0.75,0.75,0.75,0.75,0.75]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFx).flatten()) == 0)
|
||||
|
||||
def test_gridFy(self):
|
||||
x = np.r_[0.25,0.75,0.25,0.75]
|
||||
y = np.r_[0,0,1.,1.]
|
||||
z = np.r_[0.5,0.5,0.5,0.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFy).flatten()) == 0)
|
||||
|
||||
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375]
|
||||
y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1]
|
||||
z = np.r_[0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFy).flatten()) == 0)
|
||||
|
||||
def test_gridFz(self):
|
||||
x = np.r_[0.25,0.75,0.25,0.75]
|
||||
y = np.r_[0.5,0.5,0.5,0.5]
|
||||
z = np.r_[0,0,1.,1.]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFz).flatten()) == 0)
|
||||
|
||||
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375]
|
||||
y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75,0.25,0.25,0.5,0.75,0.75]
|
||||
z = np.r_[0,0,0,0,0,0.5,0.5,0.5,0.5,1,1,1,1,1]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFz).flatten()) == 0)
|
||||
|
||||
|
||||
def test_gridEx(self):
|
||||
x = np.r_[0.25,0.75,0.25,0.75,0.25,0.75,0.25,0.75]
|
||||
y = np.r_[0,0,1.,1.,0,0,1.,1.]
|
||||
z = np.r_[0,0,0,0,1.,1.,1.,1.]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEx).flatten()) == 0)
|
||||
|
||||
x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75]
|
||||
y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1,0,0,0,0.5,0.5,1,1,1]
|
||||
z = np.r_[0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEx).flatten()) == 0)
|
||||
|
||||
def test_gridEy(self):
|
||||
x = np.r_[0,0.5,1,0,0.5,1]
|
||||
y = np.r_[0.5,0.5,0.5,0.5,0.5,0.5]
|
||||
z = np.r_[0,0,0,1.,1.,1.]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEy).flatten()) == 0)
|
||||
|
||||
x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5]
|
||||
y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75,0.25,0.25,0.25,0.5,0.75,0.75,0.75]
|
||||
z = np.r_[0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEy).flatten()) == 0)
|
||||
|
||||
def test_gridEz(self):
|
||||
x = np.r_[0,0.5,1,0,0.5,1]
|
||||
y = np.r_[0,0,0,1.,1.,1.]
|
||||
z = np.r_[0.5,0.5,0.5,0.5,0.5,0.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEz).flatten()) == 0)
|
||||
|
||||
x = np.r_[0,0.25,0.5,1,0 ,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0 ,0.25,0.5,0 ,0.25,0.5]
|
||||
y = np.r_[0,0 ,0 ,0,0.5,0.5 ,0.5,1,1 ,1 ,1,0,0 ,0 ,0.5,0.5 ,0.5,1 ,1 ,1 ]
|
||||
z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEz).flatten()) == 0)
|
||||
|
||||
|
||||
class TestQuadTreeObjects(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.M = TreeMesh([2,1])
|
||||
self.Mr = TreeMesh([2,1])
|
||||
self.Mr.children[0,0].refine()
|
||||
self.Mr.number()
|
||||
# self.Mr.plotGrid(showIt=True)
|
||||
|
||||
def test_pointersM(self):
|
||||
c0 = self.M.children[0,0]
|
||||
c0fXm = c0.fXm
|
||||
c0fXp = c0.fXp
|
||||
c0fYm = c0.fYm
|
||||
c0fYp = c0.fYp
|
||||
|
||||
c1 = self.M.children[1,0]
|
||||
c1fXm = c1.fXm
|
||||
c1fXp = c1.fXp
|
||||
c1fYm = c1.fYm
|
||||
c1fYp = c1.fYp
|
||||
|
||||
self.assertTrue(c0fXp is c1fXm)
|
||||
self.assertTrue(c0fYp is not c1fYm)
|
||||
self.assertTrue(c0fXm is not c1fXm)
|
||||
|
||||
self.assertTrue(c0fXm.area == 1)
|
||||
self.assertTrue(c0fYm.area == 0.5)
|
||||
|
||||
self.assertTrue(c0.node1 is c1.node0)
|
||||
self.assertTrue(c0.node3 is c1.node2)
|
||||
self.assertTrue(self.M.nN == 6)
|
||||
|
||||
|
||||
def test_pointersMr(self):
|
||||
c0 = self.Mr.sortedCells[0]
|
||||
c0fXm = c0.fXm
|
||||
c0fXp = c0.fXp
|
||||
c0fYm = c0.fYm
|
||||
c0fYp = c0.fYp
|
||||
|
||||
c1 = self.Mr.sortedCells[1]
|
||||
c1fXm = c1.fXm
|
||||
c1fXp = c1.fXp
|
||||
c1fYm = c1.fYm
|
||||
c1fYp = c1.fYp
|
||||
|
||||
c2 = self.Mr.sortedCells[2]
|
||||
c2fXm = c2.fXm
|
||||
c2fXp = c2.fXp
|
||||
c2fYm = c2.fYm
|
||||
c2fYp = c2.fYp
|
||||
|
||||
c4 = self.Mr.sortedCells[4]
|
||||
c4fXm = c4.fXm
|
||||
c4fXp = c4.fXp
|
||||
c4fYm = c4.fYm
|
||||
c4fYp = c4.fYp
|
||||
|
||||
self.assertTrue(c0fXp is c1fXm)
|
||||
self.assertTrue(c1fXp.node0 is c2fXm.node0)
|
||||
self.assertTrue(c1fXp.node0 is c2fXm.node0)
|
||||
self.assertTrue(c4fYm is c1fYp)
|
||||
self.assertTrue(c4fXp.node1 is c2fXm.node1)
|
||||
self.assertTrue(c4fXp.node0 is c1fYp.node1)
|
||||
self.assertTrue(c0fXp.node1 is c4fYm.node0)
|
||||
|
||||
self.assertTrue(self.Mr.nN == 11)
|
||||
|
||||
self.assertTrue(np.all(c1fXp.node0.x0 == np.r_[0.5,0]))
|
||||
self.assertTrue(np.all(c1fYp.node0.x0 == np.r_[0.25,0.5]))
|
||||
|
||||
|
||||
class TestQuadTreeMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
M = TreeMesh([np.ones(x) for x in [3,2]])
|
||||
for ii in range(1):
|
||||
M.children[ii,ii].refine()
|
||||
self.M = M
|
||||
M.number()
|
||||
# M.plotGrid(showIt=True)
|
||||
|
||||
def test_MeshSizes(self):
|
||||
self.assertTrue(self.M.nC==9)
|
||||
self.assertTrue(self.M.nF==25)
|
||||
self.assertTrue(self.M.nFx==12)
|
||||
self.assertTrue(self.M.nFy==13)
|
||||
self.assertTrue(self.M.nE==25)
|
||||
self.assertTrue(self.M.nEx==13)
|
||||
self.assertTrue(self.M.nEy==12)
|
||||
|
||||
def test_gridCC(self):
|
||||
x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.5,1.5,2.5]
|
||||
y = np.r_[0.25,0.25,0.5,0.5,0.75,0.75,1.5,1.5,1.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridCC).flatten()) == 0)
|
||||
|
||||
def test_gridN(self):
|
||||
x = np.r_[0,0.5,1,2,3,0,0.5,1,0,0.5,1,2,3,0,1,2,3]
|
||||
y = np.r_[0,0,0,0,0,.5,.5,.5,1,1,1,1,1,2,2,2,2]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridN).flatten()) == 0)
|
||||
|
||||
def test_gridFx(self):
|
||||
x = np.r_[0.0,0.5,1.0,2.0,3.0,0.0,0.5,1.0,0.0,1.0,2.0,3.0]
|
||||
y = np.r_[0.25,0.25,0.25,0.5,0.5,0.75,0.75,0.75,1.5,1.5,1.5,1.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridFx).flatten()) == 0)
|
||||
|
||||
def test_gridFy(self):
|
||||
x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.25,0.75,1.5,2.5,0.5,1.5,2.5]
|
||||
y = np.r_[0,0,0,0,0.5,0.5,1,1,1,1,2,2,2]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridFy).flatten()) == 0)
|
||||
|
||||
def test_gridEx(self):
|
||||
x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.25,0.75,1.5,2.5,0.5,1.5,2.5]
|
||||
y = np.r_[0,0,0,0,0.5,0.5,1,1,1,1,2,2,2]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridEx).flatten()) == 0)
|
||||
|
||||
def test_gridEy(self):
|
||||
x = np.r_[0.0,0.5,1.0,2.0,3.0,0.0,0.5,1.0,0.0,1.0,2.0,3.0]
|
||||
y = np.r_[0.25,0.25,0.25,0.5,0.5,0.75,0.75,0.75,1.5,1.5,1.5,1.5]
|
||||
self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridEy).flatten()) == 0)
|
||||
|
||||
|
||||
class SimpleOctreeOperatorTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
h1 = np.random.rand(5)
|
||||
h2 = np.random.rand(7)
|
||||
h3 = np.random.rand(3)
|
||||
self.tM = TensorMesh([h1,h2,h3])
|
||||
self.oM = TreeMesh([h1,h2,h3])
|
||||
self.tM2 = TensorMesh([h1,h2])
|
||||
self.oM2 = TreeMesh([h1,h2])
|
||||
|
||||
def test_faceDiv(self):
|
||||
self.assertAlmostEqual((self.tM.faceDiv - self.oM.faceDiv).toarray().sum(), 0)
|
||||
self.assertAlmostEqual((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum(), 0)
|
||||
|
||||
def test_nodalGrad(self):
|
||||
self.assertAlmostEqual((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum(), 0)
|
||||
self.assertAlmostEqual((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum(), 0)
|
||||
|
||||
def test_edgeCurl(self):
|
||||
self.assertAlmostEqual((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum(), 0)
|
||||
# self.assertAlmostEqual((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum(), 0)
|
||||
|
||||
def test_InnerProducts(self):
|
||||
self.assertAlmostEqual((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum(), 0)
|
||||
self.assertAlmostEqual((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum(), 0)
|
||||
self.assertAlmostEqual((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum(), 0)
|
||||
self.assertAlmostEqual((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum(), 0)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,202 @@
|
||||
import unittest
|
||||
import sys
|
||||
from SimPEG.Mesh.BaseMesh import BaseRectangularMesh
|
||||
import numpy as np
|
||||
|
||||
|
||||
class TestBaseMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.mesh = BaseRectangularMesh([6, 2, 3])
|
||||
|
||||
def test_meshDimensions(self):
|
||||
self.assertTrue(self.mesh.dim, 3)
|
||||
|
||||
def test_mesh_nc(self):
|
||||
self.assertTrue(self.mesh.nC == 36)
|
||||
self.assertTrue(np.all(self.mesh.vnC == [6, 2, 3]))
|
||||
|
||||
def test_mesh_nc_xyz(self):
|
||||
self.assertTrue(np.all(self.mesh.nCx == 6))
|
||||
self.assertTrue(np.all(self.mesh.nCy == 2))
|
||||
self.assertTrue(np.all(self.mesh.nCz == 3))
|
||||
|
||||
def test_mesh_nf(self):
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [7, 2, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [6, 3, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnFz == [6, 2, 4]))
|
||||
|
||||
def test_mesh_ne(self):
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [6, 3, 4]))
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [7, 2, 4]))
|
||||
self.assertTrue(np.all(self.mesh.vnEz == [7, 3, 3]))
|
||||
|
||||
def test_mesh_numbers(self):
|
||||
self.assertTrue(self.mesh.nC == 36)
|
||||
self.assertTrue(np.all(self.mesh.vnF == [42, 54, 48]))
|
||||
self.assertTrue(np.all(self.mesh.vnE == [72, 56, 63]))
|
||||
self.assertTrue(np.all(self.mesh.nF == np.sum([42, 54, 48])))
|
||||
self.assertTrue(np.all(self.mesh.nE == np.sum([72, 56, 63])))
|
||||
|
||||
def test_mesh_r_E_V(self):
|
||||
ex = np.ones(self.mesh.nEx)
|
||||
ey = np.ones(self.mesh.nEy)*2
|
||||
ez = np.ones(self.mesh.nEz)*3
|
||||
e = np.r_[ex, ey, ez]
|
||||
tex = self.mesh.r(e, 'E', 'Ex', 'V')
|
||||
tey = self.mesh.r(e, 'E', 'Ey', 'V')
|
||||
tez = self.mesh.r(e, 'E', 'Ez', 'V')
|
||||
self.assertTrue(np.all(tex == ex))
|
||||
self.assertTrue(np.all(tey == ey))
|
||||
self.assertTrue(np.all(tez == ez))
|
||||
tex, tey, tez = self.mesh.r(e, 'E', 'E', 'V')
|
||||
self.assertTrue(np.all(tex == ex))
|
||||
self.assertTrue(np.all(tey == ey))
|
||||
self.assertTrue(np.all(tez == ez))
|
||||
|
||||
def test_mesh_r_F_V(self):
|
||||
fx = np.ones(self.mesh.nFx)
|
||||
fy = np.ones(self.mesh.nFy)*2
|
||||
fz = np.ones(self.mesh.nFz)*3
|
||||
f = np.r_[fx, fy, fz]
|
||||
tfx = self.mesh.r(f, 'F', 'Fx', 'V')
|
||||
tfy = self.mesh.r(f, 'F', 'Fy', 'V')
|
||||
tfz = self.mesh.r(f, 'F', 'Fz', 'V')
|
||||
self.assertTrue(np.all(tfx == fx))
|
||||
self.assertTrue(np.all(tfy == fy))
|
||||
self.assertTrue(np.all(tfz == fz))
|
||||
tfx, tfy, tfz = self.mesh.r(f, 'F', 'F', 'V')
|
||||
self.assertTrue(np.all(tfx == fx))
|
||||
self.assertTrue(np.all(tfy == fy))
|
||||
self.assertTrue(np.all(tfz == fz))
|
||||
|
||||
def test_mesh_r_E_M(self):
|
||||
g = np.ones((np.prod(self.mesh.vnEx), 3))
|
||||
g[:, 1] = 2
|
||||
g[:, 2] = 3
|
||||
Xex, Yex, Zex = self.mesh.r(g, 'Ex', 'Ex', 'M')
|
||||
self.assertTrue(np.all(Xex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Yex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Zex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Xex == 1))
|
||||
self.assertTrue(np.all(Yex == 2))
|
||||
self.assertTrue(np.all(Zex == 3))
|
||||
|
||||
def test_mesh_r_F_M(self):
|
||||
g = np.ones((np.prod(self.mesh.vnFx), 3))
|
||||
g[:, 1] = 2
|
||||
g[:, 2] = 3
|
||||
Xfx, Yfx, Zfx = self.mesh.r(g, 'Fx', 'Fx', 'M')
|
||||
self.assertTrue(np.all(Xfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Yfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Zfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Xfx == 1))
|
||||
self.assertTrue(np.all(Yfx == 2))
|
||||
self.assertTrue(np.all(Zfx == 3))
|
||||
|
||||
def test_mesh_r_CC_M(self):
|
||||
g = np.ones((self.mesh.nC, 3))
|
||||
g[:, 1] = 2
|
||||
g[:, 2] = 3
|
||||
Xc, Yc, Zc = self.mesh.r(g, 'CC', 'CC', 'M')
|
||||
self.assertTrue(np.all(Xc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Yc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Zc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Xc == 1))
|
||||
self.assertTrue(np.all(Yc == 2))
|
||||
self.assertTrue(np.all(Zc == 3))
|
||||
|
||||
|
||||
class TestMeshNumbers2D(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.mesh = BaseRectangularMesh([6, 2])
|
||||
|
||||
def test_meshDimensions(self):
|
||||
self.assertTrue(self.mesh.dim, 2)
|
||||
|
||||
def test_mesh_nc(self):
|
||||
self.assertTrue(np.all(self.mesh.vnC == [6, 2]))
|
||||
|
||||
def test_mesh_nc_xyz(self):
|
||||
self.assertTrue(np.all(self.mesh.nCx == 6))
|
||||
self.assertTrue(np.all(self.mesh.nCy == 2))
|
||||
self.assertTrue(self.mesh.nCz is None)
|
||||
|
||||
def test_mesh_nf(self):
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [7, 2]))
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [6, 3]))
|
||||
self.assertTrue(self.mesh.vnFz is None)
|
||||
|
||||
def test_mesh_ne(self):
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [6, 3]))
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [7, 2]))
|
||||
self.assertTrue(self.mesh.vnEz is None)
|
||||
|
||||
def test_mesh_numbers(self):
|
||||
c = self.mesh.nC == 12
|
||||
self.assertTrue(np.all(self.mesh.vnF == [14, 18]))
|
||||
self.assertTrue(np.all(self.mesh.nFx == 14))
|
||||
self.assertTrue(np.all(self.mesh.nFy == 18))
|
||||
self.assertTrue(np.all(self.mesh.nEx == 18))
|
||||
self.assertTrue(np.all(self.mesh.nEy == 14))
|
||||
self.assertTrue(np.all(self.mesh.vnE == [18, 14]))
|
||||
self.assertTrue(np.all(self.mesh.vnE == [18, 14]))
|
||||
self.assertTrue(np.all(self.mesh.nF == np.sum([14, 18])))
|
||||
self.assertTrue(np.all(self.mesh.nE == np.sum([18, 14])))
|
||||
|
||||
def test_mesh_r_E_V(self):
|
||||
ex = np.ones(self.mesh.nEx)
|
||||
ey = np.ones(self.mesh.nEy)*2
|
||||
e = np.r_[ex, ey]
|
||||
tex = self.mesh.r(e, 'E', 'Ex', 'V')
|
||||
tey = self.mesh.r(e, 'E', 'Ey', 'V')
|
||||
self.assertTrue(np.all(tex == ex))
|
||||
self.assertTrue(np.all(tey == ey))
|
||||
tex, tey = self.mesh.r(e, 'E', 'E', 'V')
|
||||
self.assertTrue(np.all(tex == ex))
|
||||
self.assertTrue(np.all(tey == ey))
|
||||
self.assertRaises(AssertionError, self.mesh.r, e, 'E', 'Ez', 'V')
|
||||
|
||||
def test_mesh_r_F_V(self):
|
||||
fx = np.ones(self.mesh.nFx)
|
||||
fy = np.ones(self.mesh.nFy)*2
|
||||
f = np.r_[fx, fy]
|
||||
tfx = self.mesh.r(f, 'F', 'Fx', 'V')
|
||||
tfy = self.mesh.r(f, 'F', 'Fy', 'V')
|
||||
self.assertTrue(np.all(tfx == fx))
|
||||
self.assertTrue(np.all(tfy == fy))
|
||||
tfx, tfy = self.mesh.r(f, 'F', 'F', 'V')
|
||||
self.assertTrue(np.all(tfx == fx))
|
||||
self.assertTrue(np.all(tfy == fy))
|
||||
self.assertRaises(AssertionError, self.mesh.r, f, 'F', 'Fz', 'V')
|
||||
|
||||
def test_mesh_r_E_M(self):
|
||||
g = np.ones((np.prod(self.mesh.vnEx), 2))
|
||||
g[:, 1] = 2
|
||||
Xex, Yex = self.mesh.r(g, 'Ex', 'Ex', 'M')
|
||||
self.assertTrue(np.all(Xex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Yex.shape == self.mesh.vnEx))
|
||||
self.assertTrue(np.all(Xex == 1))
|
||||
self.assertTrue(np.all(Yex == 2))
|
||||
|
||||
def test_mesh_r_F_M(self):
|
||||
g = np.ones((np.prod(self.mesh.vnFx), 2))
|
||||
g[:, 1] = 2
|
||||
Xfx, Yfx = self.mesh.r(g, 'Fx', 'Fx', 'M')
|
||||
self.assertTrue(np.all(Xfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Yfx.shape == self.mesh.vnFx))
|
||||
self.assertTrue(np.all(Xfx == 1))
|
||||
self.assertTrue(np.all(Yfx == 2))
|
||||
|
||||
def test_mesh_r_CC_M(self):
|
||||
g = np.ones((self.mesh.nC, 2))
|
||||
g[:, 1] = 2
|
||||
Xc, Yc = self.mesh.r(g, 'CC', 'CC', 'M')
|
||||
self.assertTrue(np.all(Xc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Yc.shape == self.mesh.vnC))
|
||||
self.assertTrue(np.all(Xc == 1))
|
||||
self.assertTrue(np.all(Yc == 2))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,500 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG import *
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh']
|
||||
|
||||
class Test1D_InhomogeneousDirichlet(Tests.OrderTest):
|
||||
name = "1D - Dirichlet"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 1
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.cos(np.pi*x)
|
||||
j_fun = lambda x: -np.pi*np.sin(np.pi*x)
|
||||
q_fun = lambda x: -(np.pi**2)*np.cos(np.pi*x)
|
||||
|
||||
xc_ana = phi(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
j_ana = j_fun(self.M.gridFx)
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
# vec = self.M.vectorNx
|
||||
vec = self.M.vectorCCx
|
||||
|
||||
phi_bc = phi(vec[[0,-1]])
|
||||
j_bc = j_fun(vec[[0,-1]])
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'dirichlet']])
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
V = Utils.sdiag(self.M.vol)
|
||||
G = -Pin.T*Pin*self.M.faceDiv.T * V
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_ana + P*phi_bc)
|
||||
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = V*D*Pin.T*Pin*McI*G
|
||||
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
|
||||
# A = D*McI*G
|
||||
# rhs = q_ana - D*McI*P*phi_bc
|
||||
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((j-j_ana), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-V*q_ana), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
solver = SolverCG(A, maxiter=1000)
|
||||
xc = solver * (rhs)
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc = Solver(A) * (rhs)
|
||||
print np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
j = McI*(G*xc + P*phi_bc)
|
||||
err = np.linalg.norm((j-j_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderX(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Inverse"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "1D - InhomogeneousDirichlet_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class Test2D_InhomogeneousDirichlet(Tests.OrderTest):
|
||||
name = "2D - Dirichlet"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
j_funX = lambda x: -np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
j_funY = lambda x: -np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
|
||||
q_fun = lambda x: -2*(np.pi**2)*phi(x)
|
||||
|
||||
xc_ana = phi(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
jX_ana = j_funX(self.M.gridFx)
|
||||
jY_ana = j_funY(self.M.gridFy)
|
||||
j_ana = np.r_[jX_ana,jY_ana]
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
# fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
# gBFx = self.M.gridFx[(fxm|fxp),:]
|
||||
# gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
fxm,fxp,fym,fyp = self.M.cellBoundaryInd
|
||||
gBFx = self.M.gridCC[(fxm|fxp),:]
|
||||
gBFy = self.M.gridCC[(fym|fyp),:]
|
||||
|
||||
bc = phi(np.r_[gBFx,gBFy])
|
||||
|
||||
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF('dirichlet')
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol)
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_ana + P*bc)
|
||||
q = D*j
|
||||
|
||||
# self.M.plotImage(j, 'FxFy', showIt=True)
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = D*McI*G
|
||||
rhs = q_ana - D*McI*P*bc
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((j-j_ana), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-q_ana), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
xc = Solver(A) * (rhs)
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
elif self.myTest == 'xcJ':
|
||||
xc = Solver(A) * (rhs)
|
||||
j = McI*(G*xc + P*bc)
|
||||
err = np.linalg.norm((j-j_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderX(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Inverse"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "2D - InhomogeneousDirichlet_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
class Test1D_InhomogeneousNeumann(Tests.OrderTest):
|
||||
name = "1D - Neumann"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 1
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.sin(np.pi*x)
|
||||
j_fun = lambda x: np.pi*np.cos(np.pi*x)
|
||||
q_fun = lambda x: -(np.pi**2)*np.sin(np.pi*x)
|
||||
|
||||
xc_ana = phi(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
j_ana = j_fun(self.M.gridFx)
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
vecN = self.M.vectorNx
|
||||
vecC = self.M.vectorCCx
|
||||
|
||||
phi_bc = phi(vecC[[0,-1]])
|
||||
j_bc = j_fun(vecN[[0,-1]])
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF([['neumann', 'neumann']])
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
V = Utils.sdiag(self.M.vol)
|
||||
G = -Pin.T*Pin*self.M.faceDiv.T * V
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_ana + P*phi_bc)
|
||||
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = V*D*Pin.T*Pin*McI*G
|
||||
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
|
||||
# A = D*McI*G
|
||||
# rhs = q_ana - D*McI*P*phi_bc
|
||||
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-V*q_ana), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
j = McI*(G*xc + P*phi_bc)
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "1D - InhomogeneousNeumann_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "1D - InhomogeneousNeumann_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "1D - InhomogeneousNeumann_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
class Test2D_InhomogeneousNeumann(Tests.OrderTest):
|
||||
name = "2D - Neumann"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
# meshSizes = [4]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.sin(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
|
||||
j_funX = lambda x: np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
|
||||
j_funY = lambda x: np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
q_fun = lambda x: -2*(np.pi**2)*phi(x)
|
||||
|
||||
xc_ana = phi(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
jX_ana = j_funX(self.M.gridFx)
|
||||
jY_ana = j_funY(self.M.gridFy)
|
||||
j_ana = np.r_[jX_ana,jY_ana]
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
|
||||
cxm,cxp,cym,cyp = self.M.cellBoundaryInd
|
||||
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
|
||||
gBFx = self.M.gridFx[(fxm|fxp),:]
|
||||
gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
|
||||
gBCx = self.M.gridCC[(cxm|cxp),:]
|
||||
gBCy = self.M.gridCC[(cym|cyp),:]
|
||||
|
||||
phi_bc = phi(np.r_[gBFx,gBFy])
|
||||
j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)]
|
||||
|
||||
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF('neumann')
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
V = Utils.sdiag(self.M.vol)
|
||||
G = -Pin.T*Pin*self.M.faceDiv.T * V
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_ana + P*phi_bc)
|
||||
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = V*D*Pin.T*Pin*McI*G
|
||||
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-V*q_ana), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
j = McI*(G*xc + P*phi_bc)
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "2D - InhomogeneousNeumann_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "2D - InhomogeneousNeumann_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "2D - InhomogeneousNeumann_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
class Test1D_InhomogeneousMixed(Tests.OrderTest):
|
||||
name = "1D - Mixed"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 1
|
||||
expectedOrders = 2
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.cos(0.5*np.pi*x)
|
||||
j_fun = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x)
|
||||
q_fun = lambda x: -0.25*(np.pi**2)*np.cos(0.5*np.pi*x)
|
||||
|
||||
xc_ana = phi(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
j_ana = j_fun(self.M.gridFx)
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
vecN = self.M.vectorNx
|
||||
vecC = self.M.vectorCCx
|
||||
|
||||
phi_bc = phi(vecC[[0,-1]])
|
||||
j_bc = j_fun(vecN[[0,-1]])
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann']])
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
V = Utils.sdiag(self.M.vol)
|
||||
G = -Pin.T*Pin*self.M.faceDiv.T * V
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_ana + P*phi_bc)
|
||||
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = V*D*Pin.T*Pin*McI*G
|
||||
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
|
||||
# A = D*McI*G
|
||||
# rhs = q_ana - D*McI*P*phi_bc
|
||||
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-V*q_ana), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
j = McI*(G*xc + P*phi_bc)
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "1D - InhomogeneousMixed_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "1D - InhomogeneousMixed_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "1D - InhomogeneousMixed_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
class Test2D_InhomogeneousMixed(Tests.OrderTest):
|
||||
name = "2D - Mixed"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
expectedOrders = 2
|
||||
meshSizes = [2, 4, 8, 16]
|
||||
# meshSizes = [4]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi = lambda x: np.cos(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1])
|
||||
j_funX = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1])
|
||||
j_funY = lambda x: -0.5*np.pi*np.cos(0.5*np.pi*x[:,0])*np.sin(0.5*np.pi*x[:,1])
|
||||
q_fun = lambda x: -2*((0.5*np.pi)**2)*phi(x)
|
||||
|
||||
xc_ana = phi(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
jX_ana = j_funX(self.M.gridFx)
|
||||
jY_ana = j_funY(self.M.gridFy)
|
||||
j_ana = np.r_[jX_ana,jY_ana]
|
||||
|
||||
#TODO: Check where our boundary conditions are CCx or Nx
|
||||
|
||||
cxm,cxp,cym,cyp = self.M.cellBoundaryInd
|
||||
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
|
||||
gBFx = self.M.gridFx[(fxm|fxp),:]
|
||||
gBFy = self.M.gridFy[(fym|fyp),:]
|
||||
|
||||
gBCx = self.M.gridCC[(cxm|cxp),:]
|
||||
gBCy = self.M.gridCC[(cym|cyp),:]
|
||||
|
||||
phi_bc = phi(np.r_[gBCx,gBCy])
|
||||
j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)]
|
||||
|
||||
# P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
|
||||
|
||||
P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann'], ['dirichlet', 'neumann']])
|
||||
|
||||
Mc = self.M.getFaceInnerProduct()
|
||||
McI = Utils.sdInv(self.M.getFaceInnerProduct())
|
||||
V = Utils.sdiag(self.M.vol)
|
||||
G = -Pin.T*Pin*self.M.faceDiv.T * V
|
||||
D = self.M.faceDiv
|
||||
j = McI*(G*xc_ana + P*phi_bc)
|
||||
q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
|
||||
|
||||
# Rearrange if we know q to solve for x
|
||||
A = V*D*Pin.T*Pin*McI*G
|
||||
rhs = V*q_ana - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
|
||||
|
||||
if self.myTest == 'j':
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
elif self.myTest == 'q':
|
||||
err = np.linalg.norm((q-V*q_ana), np.inf)
|
||||
elif self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
elif self.myTest == 'xcJ':
|
||||
#TODO: fix the null space
|
||||
xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
|
||||
j = McI*(G*xc + P*phi_bc)
|
||||
err = np.linalg.norm((Pin*j-Pin*j_ana), np.inf)
|
||||
if info > 0:
|
||||
print 'Solve does not work well'
|
||||
print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
|
||||
return err
|
||||
|
||||
def test_orderJ(self):
|
||||
self.name = "2D - InhomogeneousMixed_Forward j"
|
||||
self.myTest = 'j'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderQ(self):
|
||||
self.name = "2D - InhomogeneousMixed_Forward q"
|
||||
self.myTest = 'q'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderXJ(self):
|
||||
self.name = "2D - InhomogeneousMixed_Inverse J"
|
||||
self.myTest = 'xcJ'
|
||||
self.orderTest()
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,400 @@
|
||||
import unittest
|
||||
import sys
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class TestCyl2DMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
hx = np.r_[1,1,0.5]
|
||||
hz = np.r_[2,1]
|
||||
self.mesh = Mesh.CylMesh([hx, 1,hz])
|
||||
|
||||
def test_dim(self):
|
||||
self.assertTrue(self.mesh.dim == 3)
|
||||
|
||||
def test_nC(self):
|
||||
self.assertTrue(self.mesh.nC == 6)
|
||||
self.assertTrue(self.mesh.nCx == 3)
|
||||
self.assertTrue(self.mesh.nCy == 1)
|
||||
self.assertTrue(self.mesh.nCz == 2)
|
||||
self.assertTrue(np.all(self.mesh.vnC == [3, 1, 2]))
|
||||
|
||||
def test_nN(self):
|
||||
self.assertTrue(self.mesh.nN == 0)
|
||||
self.assertTrue(self.mesh.nNx == 3)
|
||||
self.assertTrue(self.mesh.nNy == 0)
|
||||
self.assertTrue(self.mesh.nNz == 3)
|
||||
self.assertTrue(np.all(self.mesh.vnN == [3, 0, 3]))
|
||||
|
||||
def test_nF(self):
|
||||
self.assertTrue(self.mesh.nFx == 6)
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [3, 1, 2]))
|
||||
self.assertTrue(self.mesh.nFy == 0)
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [3, 0, 2]))
|
||||
self.assertTrue(self.mesh.nFz == 9)
|
||||
self.assertTrue(np.all(self.mesh.vnFz == [3, 1, 3]))
|
||||
self.assertTrue(self.mesh.nF == 15)
|
||||
self.assertTrue(np.all(self.mesh.vnF == [6, 0, 9]))
|
||||
|
||||
def test_nE(self):
|
||||
self.assertTrue(self.mesh.nEx == 0)
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [3, 0, 3]))
|
||||
self.assertTrue(self.mesh.nEy == 9)
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [3, 1, 3]))
|
||||
self.assertTrue(self.mesh.nEz == 0)
|
||||
self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2]))
|
||||
self.assertTrue(self.mesh.nE == 9)
|
||||
self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0]))
|
||||
|
||||
def test_vectorsCC(self):
|
||||
v = np.r_[0.5, 1.5, 2.25]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0)
|
||||
v = np.r_[0]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0)
|
||||
v = np.r_[1, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0)
|
||||
|
||||
def test_vectorsN(self):
|
||||
v = np.r_[1, 2, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0)
|
||||
v = np.r_[0]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0)
|
||||
v = np.r_[0, 2, 3.]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0)
|
||||
|
||||
def test_edge(self):
|
||||
edge = np.r_[1, 2, 2.5, 1, 2, 2.5, 1, 2, 2.5] * 2 * np.pi
|
||||
self.assertTrue(np.linalg.norm((edge-self.mesh.edge)) == 0)
|
||||
|
||||
def test_area(self):
|
||||
r = np.r_[0, 1, 2, 2.5]
|
||||
a = r[1:]*2*np.pi
|
||||
areaX = np.r_[2*a,a]
|
||||
a = (r[1:]**2 - r[:-1]**2)*np.pi
|
||||
areaZ = np.r_[a,a,a]
|
||||
area = np.r_[areaX, areaZ]
|
||||
self.assertTrue(np.linalg.norm((area-self.mesh.area)) == 0)
|
||||
|
||||
def test_vol(self):
|
||||
r = np.r_[0, 1, 2, 2.5]
|
||||
a = (r[1:]**2 - r[:-1]**2)*np.pi
|
||||
vol = np.r_[2*a,a]
|
||||
self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0)
|
||||
|
||||
def test_gridSizes(self):
|
||||
self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3))
|
||||
self.assertTrue(self.mesh.gridN.shape == (9, 3))
|
||||
|
||||
self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3))
|
||||
self.assertTrue(self.mesh.gridFy is None)
|
||||
self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3))
|
||||
|
||||
self.assertTrue(self.mesh.gridEx is None)
|
||||
self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3))
|
||||
self.assertTrue(self.mesh.gridEz is None)
|
||||
|
||||
def test_gridCC(self):
|
||||
x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25]
|
||||
y = np.zeros(6)
|
||||
z = np.r_[1,1,1,2.5,2.5,2.5]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridCC).ravel()) == 0)
|
||||
|
||||
def test_gridN(self):
|
||||
x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5]
|
||||
y = np.zeros(9)
|
||||
z = np.r_[0,0,0,2,2,2,3,3,3.]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridN).ravel()) == 0)
|
||||
|
||||
def test_gridFx(self):
|
||||
x = np.r_[1,2,2.5,1,2,2.5]
|
||||
y = np.zeros(6)
|
||||
z = np.r_[1,1,1,2.5,2.5,2.5]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridFx).ravel()) == 0)
|
||||
|
||||
def test_gridFz(self):
|
||||
x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25,0.5,1.5,2.25]
|
||||
y = np.zeros(9)
|
||||
z = np.r_[0,0,0,2,2,2,3,3,3.]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridFz).ravel()) == 0)
|
||||
|
||||
def test_gridEy(self):
|
||||
x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5]
|
||||
y = np.zeros(9)
|
||||
z = np.r_[0,0,0,2,2,2,3,3,3.]
|
||||
G = np.c_[x,y,z]
|
||||
self.assertTrue(np.linalg.norm((G-self.mesh.gridEy).ravel()) == 0)
|
||||
|
||||
def test_lightOperators(self):
|
||||
self.assertTrue(self.mesh.nodalGrad is None)
|
||||
|
||||
def test_getInterpMatCartMesh_Cells(self):
|
||||
|
||||
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
|
||||
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
|
||||
|
||||
mc = np.arange(Mc.nC)
|
||||
xr = np.linspace(0,0.4,50)
|
||||
xc = np.linspace(0,0.4,50) + 0.2
|
||||
Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'CC')
|
||||
Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
|
||||
Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC')
|
||||
|
||||
assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
|
||||
|
||||
def test_getInterpMatCartMesh_Faces(self):
|
||||
|
||||
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
|
||||
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
|
||||
|
||||
Pf = Mc.getInterpolationMatCartMesh(Mr, 'F')
|
||||
mf = np.ones(Mc.nF)
|
||||
|
||||
frect = Pf * mf
|
||||
|
||||
fxcc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
|
||||
fycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
|
||||
fzcc = Mr.r(frect, 'F', 'Fz')
|
||||
|
||||
indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
|
||||
indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
|
||||
|
||||
TOL = 1e-2
|
||||
assert np.abs(float(fxcc[indX]) - 1) < TOL
|
||||
assert np.abs(float(fxcc[indY]) - 0) < TOL
|
||||
assert np.abs(float(fycc[indX]) - 0) < TOL
|
||||
assert np.abs(float(fycc[indY]) - 1) < TOL
|
||||
assert np.abs((fzcc - 1).sum()) < TOL
|
||||
|
||||
mag = (fxcc**2 + fycc**2)**0.5
|
||||
dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5
|
||||
|
||||
assert np.abs(mag[dist > 0.1].max() - 1) < TOL
|
||||
assert np.abs(mag[dist > 0.1].min() - 1) < TOL
|
||||
|
||||
|
||||
def test_getInterpMatCartMesh_Edges(self):
|
||||
|
||||
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
|
||||
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
|
||||
|
||||
Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
|
||||
me = np.ones(Mc.nE)
|
||||
|
||||
erect = Pe * me
|
||||
|
||||
excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
|
||||
eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
|
||||
ezcc = Mr.r(erect, 'E', 'Ez')
|
||||
|
||||
indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
|
||||
indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
|
||||
|
||||
TOL = 1e-2
|
||||
assert np.abs(float(excc[indX]) - 0) < TOL
|
||||
assert np.abs(float(excc[indY]) + 1) < TOL
|
||||
assert np.abs(float(eycc[indX]) - 1) < TOL
|
||||
assert np.abs(float(eycc[indY]) - 0) < TOL
|
||||
assert np.abs(ezcc.sum()) < TOL
|
||||
|
||||
mag = (excc**2 + eycc**2)**0.5
|
||||
dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5
|
||||
|
||||
assert np.abs(mag[dist > 0.1].max() - 1) < TOL
|
||||
assert np.abs(mag[dist > 0.1].min() - 1) < TOL
|
||||
|
||||
|
||||
MESHTYPES = ['uniformCylMesh']
|
||||
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 2])
|
||||
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
cyl_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
|
||||
cyl_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
|
||||
cylF2 = lambda M, fx, fy: np.vstack((cyl_row2(M.gridFx, fx, fy), cyl_row2(M.gridFz, fx, fy)))
|
||||
|
||||
|
||||
class TestFaceDiv2D(Tests.OrderTest):
|
||||
name = "FaceDiv"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
|
||||
funR = lambda r, z: np.sin(2.*np.pi*r)
|
||||
funZ = lambda r, z: np.sin(2.*np.pi*z)
|
||||
|
||||
sol = lambda r, t, z: (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r + 2*np.pi*np.cos(2*np.pi*z)
|
||||
|
||||
Fc = cylF2(self.M, funR, funZ)
|
||||
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
divF = self.M.faceDiv.dot(F)
|
||||
divF_ana = call3(sol, self.M.gridCC)
|
||||
|
||||
err = np.linalg.norm((divF-divF_ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestEdgeCurl2D(Tests.OrderTest):
|
||||
name = "EdgeCurl"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
# To Recreate or change the functions:
|
||||
|
||||
# import sympy
|
||||
# r,t,z = sympy.symbols('r,t,z')
|
||||
|
||||
# fR = 0
|
||||
# fZ = 0
|
||||
# fT = sympy.sin(2.*sympy.pi*z)
|
||||
|
||||
# print 1/r*sympy.diff(fZ,t) - sympy.diff(fT,z)
|
||||
# print sympy.diff(fR,z) - sympy.diff(fZ,r)
|
||||
# print 1/r*(sympy.diff(r*fT,r) - sympy.diff(fR,t))
|
||||
|
||||
funT = lambda r, t, z: np.sin(2.*np.pi*z)
|
||||
|
||||
solR = lambda r, z: -2.0*np.pi*np.cos(2.0*np.pi*z)
|
||||
solZ = lambda r, z: np.sin(2.0*np.pi*z)/r
|
||||
|
||||
E = call3(funT, self.M.gridEy)
|
||||
|
||||
curlE = self.M.edgeCurl.dot(E)
|
||||
|
||||
Fc = cylF2(self.M, solR, solZ)
|
||||
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
|
||||
curlE_ana = self.M.projectFaceVector(Fc)
|
||||
|
||||
err = np.linalg.norm((curlE-curlE_ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
# class TestInnerProducts2D(Tests.OrderTest):
|
||||
# """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
# meshTypes = MESHTYPES
|
||||
# meshDimension = 2
|
||||
# meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
# def getError(self):
|
||||
|
||||
# funR = lambda r, t, z: np.cos(2.0*np.pi*z)
|
||||
# funT = lambda r, t, z: 0*t
|
||||
# funZ = lambda r, t, z: np.sin(2.0*np.pi*r)
|
||||
|
||||
# call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
|
||||
# sigma1 = lambda r, t, z: z+1
|
||||
# sigma2 = lambda r, t, z: r*z+50
|
||||
# sigma3 = lambda r, t, z: 3+t*r
|
||||
# sigma4 = lambda r, t, z: 0.1*r*t*z
|
||||
# sigma5 = lambda r, t, z: 0.2*z*r*t
|
||||
# sigma6 = lambda r, t, z: 0.1*t
|
||||
|
||||
# Gc = self.M.gridCC
|
||||
# if self.sigmaTest == 1:
|
||||
# sigma = np.c_[call(sigma1, Gc)]
|
||||
# analytic = 144877./360 # Found using sympy. z=5
|
||||
# elif self.sigmaTest == 2:
|
||||
# sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
||||
# analytic = 189959./120 # Found using sympy. z=5
|
||||
# elif self.sigmaTest == 3:
|
||||
# sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
# analytic = 781427./360 # Found using sympy. z=5
|
||||
|
||||
# if self.location == 'edges':
|
||||
# E = call(funT, self.M.gridEy)
|
||||
# A = self.M.getEdgeInnerProduct(sigma)
|
||||
# numeric = E.T.dot(A.dot(E))
|
||||
# elif self.location == 'faces':
|
||||
# Fr = call(funR, self.M.gridFx)
|
||||
# Fz = call(funZ, self.M.gridFz)
|
||||
# A = self.M.getFaceInnerProduct(sigma)
|
||||
# F = np.r_[Fr,Fz]
|
||||
# numeric = F.T.dot(A.dot(F))
|
||||
|
||||
# print numeric
|
||||
# err = np.abs(numeric - analytic)
|
||||
# return err
|
||||
|
||||
# def test_order1_faces(self):
|
||||
# self.name = "2D Face Inner Product - Isotropic"
|
||||
# self.location = 'faces'
|
||||
# self.sigmaTest = 1
|
||||
# self.orderTest()
|
||||
|
||||
|
||||
class TestCyl3DMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
hx = np.r_[1,1,0.5]
|
||||
hy = np.r_[np.pi, np.pi]
|
||||
hz = np.r_[2,1]
|
||||
self.mesh = Mesh.CylMesh([hx, hy,hz])
|
||||
|
||||
def test_dim(self):
|
||||
self.assertTrue(self.mesh.dim == 3)
|
||||
|
||||
def test_nC(self):
|
||||
self.assertTrue(self.mesh.nCx == 3)
|
||||
self.assertTrue(self.mesh.nCy == 2)
|
||||
self.assertTrue(self.mesh.nCz == 2)
|
||||
self.assertTrue(np.all(self.mesh.vnC == [3, 2, 2]))
|
||||
|
||||
def test_nN(self):
|
||||
self.assertTrue(self.mesh.nN == 24)
|
||||
self.assertTrue(self.mesh.nNx == 4)
|
||||
self.assertTrue(self.mesh.nNy == 2)
|
||||
self.assertTrue(self.mesh.nNz == 3)
|
||||
self.assertTrue(np.all(self.mesh.vnN == [4, 2, 3]))
|
||||
|
||||
def test_nF(self):
|
||||
self.assertTrue(self.mesh.nFx == 12)
|
||||
self.assertTrue(np.all(self.mesh.vnFx == [3, 2, 2]))
|
||||
self.assertTrue(self.mesh.nFy == 12)
|
||||
self.assertTrue(np.all(self.mesh.vnFy == [3, 2, 2]))
|
||||
self.assertTrue(self.mesh.nFz == 18)
|
||||
self.assertTrue(np.all(self.mesh.vnFz == [3, 2, 3]))
|
||||
self.assertTrue(self.mesh.nF == 42)
|
||||
self.assertTrue(np.all(self.mesh.vnF == [12, 12, 18]))
|
||||
|
||||
def test_nE(self):
|
||||
self.assertTrue(self.mesh.nEx == 18)
|
||||
self.assertTrue(np.all(self.mesh.vnEx == [3, 2, 3]))
|
||||
self.assertTrue(self.mesh.nEy == 18)
|
||||
self.assertTrue(np.all(self.mesh.vnEy == [3, 2, 3]))
|
||||
self.assertTrue(self.mesh.nEz == 12 + 2)
|
||||
self.assertTrue(self.mesh.vnEz is None)
|
||||
self.assertTrue(self.mesh.nE == 50)
|
||||
self.assertTrue(np.all(self.mesh.vnE == [18, 18, 14]))
|
||||
|
||||
def test_vectorsCC(self):
|
||||
v = np.r_[0.5, 1.5, 2.25]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0)
|
||||
v = np.r_[0, np.pi]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0)
|
||||
v = np.r_[1, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0)
|
||||
|
||||
def test_vectorsN(self):
|
||||
v = np.r_[0, 1, 2, 2.5]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0)
|
||||
v = np.r_[np.pi/2, 1.5*np.pi]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0)
|
||||
v = np.r_[0, 2, 3]
|
||||
self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,402 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG import Utils, Tests
|
||||
|
||||
|
||||
class TestInnerProducts(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
|
||||
meshDimension = 3
|
||||
meshSizes = [16, 32]
|
||||
|
||||
def getError(self):
|
||||
|
||||
call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
|
||||
ex = lambda x, y, z: x**2+y*z
|
||||
ey = lambda x, y, z: (z**2)*x+y*z
|
||||
ez = lambda x, y, z: y**2+x*z
|
||||
|
||||
sigma1 = lambda x, y, z: x*y+1
|
||||
sigma2 = lambda x, y, z: x*z+2
|
||||
sigma3 = lambda x, y, z: 3+z*y
|
||||
sigma4 = lambda x, y, z: 0.1*x*y*z
|
||||
sigma5 = lambda x, y, z: 0.2*x*y
|
||||
sigma6 = lambda x, y, z: 0.1*z
|
||||
|
||||
Gc = self.M.gridCC
|
||||
if self.sigmaTest == 1:
|
||||
sigma = np.c_[call(sigma1, Gc)]
|
||||
analytic = 647./360 # Found using sympy.
|
||||
elif self.sigmaTest == 3:
|
||||
sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
analytic = 37./12 # Found using sympy.
|
||||
elif self.sigmaTest == 6:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc),
|
||||
call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)]
|
||||
analytic = 69881./21600 # Found using sympy.
|
||||
|
||||
if self.location == 'edges':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
||||
Ec = np.vstack((cart(self.M.gridEx),
|
||||
cart(self.M.gridEy),
|
||||
cart(self.M.gridEz)))
|
||||
E = self.M.projectEdgeVector(Ec)
|
||||
|
||||
if self.invProp:
|
||||
A = self.M.getEdgeInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
||||
else:
|
||||
A = self.M.getEdgeInnerProduct(sigma)
|
||||
numeric = E.T.dot(A.dot(E))
|
||||
elif self.location == 'faces':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
|
||||
Fc = np.vstack((cart(self.M.gridFx),
|
||||
cart(self.M.gridFy),
|
||||
cart(self.M.gridFz)))
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
if self.invProp:
|
||||
A = self.M.getFaceInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
||||
else:
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T.dot(A.dot(F))
|
||||
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order1_edges(self):
|
||||
self.name = "Edge Inner Product - Isotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_edges_invProp(self):
|
||||
self.name = "Edge Inner Product - Isotropic - invProp"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_edges(self):
|
||||
self.name = "Edge Inner Product - Anisotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_edges_invProp(self):
|
||||
self.name = "Edge Inner Product - Anisotropic - invProp"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_edges(self):
|
||||
self.name = "Edge Inner Product - Full Tensor"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 6
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_edges_invProp(self):
|
||||
self.name = "Edge Inner Product - Full Tensor - invProp"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 6
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces(self):
|
||||
self.name = "Face Inner Product - Isotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces_invProp(self):
|
||||
self.name = "Face Inner Product - Isotropic - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_faces(self):
|
||||
self.name = "Face Inner Product - Anisotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_faces_invProp(self):
|
||||
self.name = "Face Inner Product - Anisotropic - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_faces(self):
|
||||
self.name = "Face Inner Product - Full Tensor"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 6
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_faces_invProp(self):
|
||||
self.name = "Face Inner Product - Full Tensor - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 6
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestInnerProducts2D(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
|
||||
meshDimension = 2
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
|
||||
z = 5 # Because 5 is just such a great number.
|
||||
|
||||
call = lambda fun, xy: fun(xy[:, 0], xy[:, 1])
|
||||
|
||||
ex = lambda x, y: x**2+y*z
|
||||
ey = lambda x, y: (z**2)*x+y*z
|
||||
|
||||
sigma1 = lambda x, y: x*y+1
|
||||
sigma2 = lambda x, y: x*z+2
|
||||
sigma3 = lambda x, y: 3+z*y
|
||||
|
||||
Gc = self.M.gridCC
|
||||
if self.sigmaTest == 1:
|
||||
sigma = np.c_[call(sigma1, Gc)]
|
||||
analytic = 144877./360 # Found using sympy. z=5
|
||||
elif self.sigmaTest == 2:
|
||||
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
||||
analytic = 189959./120 # Found using sympy. z=5
|
||||
elif self.sigmaTest == 3:
|
||||
sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
analytic = 781427./360 # Found using sympy. z=5
|
||||
|
||||
if self.location == 'edges':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
||||
Ec = np.vstack((cart(self.M.gridEx),
|
||||
cart(self.M.gridEy)))
|
||||
E = self.M.projectEdgeVector(Ec)
|
||||
if self.invProp:
|
||||
A = self.M.getEdgeInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
||||
else:
|
||||
A = self.M.getEdgeInnerProduct(sigma)
|
||||
numeric = E.T.dot(A.dot(E))
|
||||
elif self.location == 'faces':
|
||||
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
|
||||
Fc = np.vstack((cart(self.M.gridFx),
|
||||
cart(self.M.gridFy)))
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
if self.invProp:
|
||||
A = self.M.getFaceInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
|
||||
else:
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T.dot(A.dot(F))
|
||||
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order1_edges(self):
|
||||
self.name = "2D Edge Inner Product - Isotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_edges_invProp(self):
|
||||
self.name = "2D Edge Inner Product - Isotropic - invProp"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_edges(self):
|
||||
self.name = "2D Edge Inner Product - Anisotropic"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 2
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_edges_invProp(self):
|
||||
self.name = "2D Edge Inner Product - Anisotropic - invProp"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 2
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_edges(self):
|
||||
self.name = "2D Edge Inner Product - Full Tensor"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order6_edges_invProp(self):
|
||||
self.name = "2D Edge Inner Product - Full Tensor - invProp"
|
||||
self.location = 'edges'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces(self):
|
||||
self.name = "2D Face Inner Product - Isotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces_invProp(self):
|
||||
self.name = "2D Face Inner Product - Isotropic - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order2_faces(self):
|
||||
self.name = "2D Face Inner Product - Anisotropic"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 2
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order2_faces_invProp(self):
|
||||
self.name = "2D Face Inner Product - Anisotropic - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 2
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_faces(self):
|
||||
self.name = "2D Face Inner Product - Full Tensor"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order3_faces_invProp(self):
|
||||
self.name = "2D Face Inner Product - Full Tensor - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 3
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
|
||||
|
||||
class TestInnerProducts1D(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 1
|
||||
meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
|
||||
y = 12 # Because 12 is just such a great number.
|
||||
z = 5 # Because 5 is just such a great number as well!
|
||||
|
||||
call = lambda fun, x: fun(x)
|
||||
|
||||
ex = lambda x: x**2+y*z
|
||||
|
||||
sigma1 = lambda x: x*y+1
|
||||
|
||||
Gc = self.M.gridCC
|
||||
sigma = call(sigma1, Gc)
|
||||
analytic = 128011./5 # Found using sympy. y=12, z=5
|
||||
|
||||
if self.location == 'faces':
|
||||
F = call(ex, self.M.gridFx)
|
||||
if self.invProp:
|
||||
A = self.M.getFaceInnerProduct(1/sigma, invProp=True)
|
||||
else:
|
||||
A = self.M.getFaceInnerProduct(sigma)
|
||||
numeric = F.T.dot(A.dot(F))
|
||||
|
||||
err = np.abs(numeric - analytic)
|
||||
return err
|
||||
|
||||
def test_order1_faces(self):
|
||||
self.name = "1D Face Inner Product"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = False
|
||||
self.orderTest()
|
||||
|
||||
def test_order1_faces_invProp(self):
|
||||
self.name = "1D Face Inner Product - invProp"
|
||||
self.location = 'faces'
|
||||
self.sigmaTest = 1
|
||||
self.invProp = True
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
###################################################
|
||||
#### Uncomment to Reevaluate the InnerProducts ####
|
||||
###################################################
|
||||
|
||||
# if __name__ == '__main__':
|
||||
# import sympy
|
||||
|
||||
# x,y,z = sympy.symbols(['x','y','z'])
|
||||
# ex = x**2+y*z
|
||||
# ey = (z**2)*x+y*z
|
||||
# ez = y**2+x*z
|
||||
# e = sympy.Matrix([ex,ey,ez])
|
||||
|
||||
# sigma1 = x*y+1
|
||||
# sigma2 = x*z+2
|
||||
# sigma3 = 3+z*y
|
||||
# sigma4 = 0.1*x*y*z
|
||||
# sigma5 = 0.2*x*y
|
||||
# sigma6 = 0.1*z
|
||||
|
||||
# S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
|
||||
# S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
|
||||
# S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
|
||||
|
||||
# print '3D'
|
||||
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
|
||||
|
||||
# z = 5
|
||||
# ex = x**2+y*z
|
||||
# ey = (z**2)*x+y*z
|
||||
# e = sympy.Matrix([ex,ey])
|
||||
|
||||
# sigma1 = x*y+1
|
||||
# sigma2 = x*z+2
|
||||
# sigma3 = 3+z*y
|
||||
|
||||
# S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
|
||||
# S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
|
||||
# S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
|
||||
|
||||
# print '2D'
|
||||
# print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
|
||||
# print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
|
||||
# print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
|
||||
|
||||
# y = 12
|
||||
# z = 5
|
||||
# ex = x**2+y*z
|
||||
# e = ex
|
||||
|
||||
# sigma1 = x*y+1
|
||||
|
||||
# print '1D'
|
||||
# print sympy.integrate(e*sigma1*e, (x,0,1))
|
||||
@@ -0,0 +1,263 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class TestInnerProductsDerivs(unittest.TestCase):
|
||||
|
||||
def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False):
|
||||
if meshType == 'Curv':
|
||||
hRect = Utils.exampleLrmGrid(h,'rotate')
|
||||
mesh = Mesh.CurvilinearMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
mesh = Mesh.TreeMesh(h)
|
||||
elif meshType == 'Tensor':
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nF)
|
||||
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
|
||||
def fun(sig):
|
||||
M = mesh.getFaceInnerProduct(sig, invProp=invProp, invMat=invMat)
|
||||
Md = mesh.getFaceInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
|
||||
return M*v, Md(v)
|
||||
print meshType, 'Face', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
|
||||
return Tests.checkDerivative(fun, sig, num=5, plotIt=False)
|
||||
|
||||
def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
|
||||
if meshType == 'Curv':
|
||||
hRect = Utils.exampleLrmGrid(h,'rotate')
|
||||
mesh = Mesh.CurvilinearMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
mesh = Mesh.TreeMesh(h)
|
||||
elif meshType == 'Tensor':
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nE)
|
||||
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
|
||||
def fun(sig):
|
||||
M = mesh.getEdgeInnerProduct(sig, invProp=invProp, invMat=invMat)
|
||||
Md = mesh.getEdgeInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
|
||||
return M*v, Md(v)
|
||||
print meshType, 'Edge', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
|
||||
return Tests.checkDerivative(fun, sig, num=5, plotIt=False)
|
||||
|
||||
def test_FaceIP_1D_float(self):
|
||||
self.assertTrue(self.doTestFace([10],0, False, 'Tensor'))
|
||||
def test_FaceIP_2D_float(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, 'Tensor'))
|
||||
def test_FaceIP_3D_float(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tensor'))
|
||||
def test_FaceIP_1D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10],1, False, 'Tensor'))
|
||||
def test_FaceIP_2D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, 'Tensor'))
|
||||
def test_FaceIP_3D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tensor'))
|
||||
def test_FaceIP_2D_anisotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, 'Tensor'))
|
||||
def test_FaceIP_3D_anisotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tensor'))
|
||||
def test_FaceIP_2D_tensor(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False, 'Tensor'))
|
||||
def test_FaceIP_3D_tensor(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tensor'))
|
||||
|
||||
def test_FaceIP_1D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10],0, True, 'Tensor'))
|
||||
def test_FaceIP_2D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor'))
|
||||
def test_FaceIP_3D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor'))
|
||||
def test_FaceIP_1D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10],1, True, 'Tensor'))
|
||||
def test_FaceIP_2D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor'))
|
||||
def test_FaceIP_3D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor'))
|
||||
def test_FaceIP_2D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor'))
|
||||
def test_FaceIP_3D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor'))
|
||||
|
||||
def test_EdgeIP_1D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10],0, False, 'Tensor'))
|
||||
def test_EdgeIP_2D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tensor'))
|
||||
def test_EdgeIP_3D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tensor'))
|
||||
def test_EdgeIP_1D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10],1, False, 'Tensor'))
|
||||
def test_EdgeIP_2D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tensor'))
|
||||
def test_EdgeIP_3D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tensor'))
|
||||
def test_EdgeIP_2D_anisotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tensor'))
|
||||
def test_EdgeIP_3D_anisotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tensor'))
|
||||
def test_EdgeIP_2D_tensor(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tensor'))
|
||||
def test_EdgeIP_3D_tensor(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tensor'))
|
||||
|
||||
def test_EdgeIP_1D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10],0, True, 'Tensor'))
|
||||
def test_EdgeIP_2D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tensor'))
|
||||
def test_EdgeIP_3D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tensor'))
|
||||
def test_EdgeIP_1D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10],1, True, 'Tensor'))
|
||||
def test_EdgeIP_2D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tensor'))
|
||||
def test_EdgeIP_3D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tensor'))
|
||||
def test_EdgeIP_2D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tensor'))
|
||||
def test_EdgeIP_3D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tensor'))
|
||||
|
||||
|
||||
|
||||
def test_FaceIP_1D_float_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_2D_float_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_3D_float_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_1D_isotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10],1, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_2D_isotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_3D_isotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_2D_anisotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_3D_anisotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))
|
||||
|
||||
|
||||
|
||||
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_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_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_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'))
|
||||
|
||||
|
||||
|
||||
|
||||
def test_FaceIP_2D_float_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, 'Tree'))
|
||||
def test_FaceIP_3D_float_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tree'))
|
||||
def test_FaceIP_2D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, 'Tree'))
|
||||
def test_FaceIP_3D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tree'))
|
||||
def test_FaceIP_2D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, 'Tree'))
|
||||
def test_FaceIP_3D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tree'))
|
||||
def test_FaceIP_2D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False, 'Tree'))
|
||||
def test_FaceIP_3D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tree'))
|
||||
|
||||
def test_FaceIP_2D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tree'))
|
||||
def test_FaceIP_3D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tree'))
|
||||
def test_FaceIP_2D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tree'))
|
||||
def test_FaceIP_3D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tree'))
|
||||
def test_FaceIP_2D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tree'))
|
||||
def test_FaceIP_3D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tree'))
|
||||
|
||||
def test_EdgeIP_2D_float_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tree'))
|
||||
def test_EdgeIP_3D_float_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tree'))
|
||||
def test_EdgeIP_2D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tree'))
|
||||
def test_EdgeIP_3D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tree'))
|
||||
def test_EdgeIP_2D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tree'))
|
||||
def test_EdgeIP_3D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tree'))
|
||||
def test_EdgeIP_2D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tree'))
|
||||
def test_EdgeIP_3D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tree'))
|
||||
|
||||
def test_EdgeIP_2D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tree'))
|
||||
def test_EdgeIP_3D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tree'))
|
||||
def test_EdgeIP_2D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tree'))
|
||||
def test_EdgeIP_3D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tree'))
|
||||
def test_EdgeIP_2D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tree'))
|
||||
def test_EdgeIP_3D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tree'))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,302 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG.Utils import mkvc
|
||||
from SimPEG import Mesh, Tests
|
||||
import unittest
|
||||
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
|
||||
TOLERANCES = [0.9, 0.5, 0.5]
|
||||
call1 = lambda fun, xyz: fun(xyz)
|
||||
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)]
|
||||
cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
|
||||
cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
|
||||
cartF2Cyl = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFz, fx, fy)))
|
||||
cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
|
||||
cartE2Cyl = lambda M, ex, ey: cart_row2(M.gridEy, ex, ey)
|
||||
cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
|
||||
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
|
||||
|
||||
TOL = 1e-7
|
||||
|
||||
|
||||
class TestInterpolation1D(Tests.OrderTest):
|
||||
LOCS = np.random.rand(50)*0.6+0.2
|
||||
name = "Interpolation 1D"
|
||||
meshTypes = MESHTYPES
|
||||
tolerance = TOLERANCES
|
||||
meshDimension = 1
|
||||
meshSizes = [8, 16, 32, 64, 128]
|
||||
|
||||
def getError(self):
|
||||
funX = lambda x: np.cos(2*np.pi*x)
|
||||
|
||||
ana = call1(funX, self.LOCS)
|
||||
|
||||
if 'CC' == self.type:
|
||||
grid = call1(funX, self.M.gridCC)
|
||||
elif 'N' == self.type:
|
||||
grid = call1(funX, self.M.gridN)
|
||||
|
||||
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
|
||||
|
||||
err = np.linalg.norm((comp - ana), 2)
|
||||
return err
|
||||
|
||||
def test_orderCC(self):
|
||||
self.type = 'CC'
|
||||
self.name = 'Interpolation 1D: CC'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN(self):
|
||||
self.type = 'N'
|
||||
self.name = 'Interpolation 1D: N'
|
||||
self.orderTest()
|
||||
|
||||
class TestOutliersInterp1D(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
|
||||
def test_outliers(self):
|
||||
M = Mesh.TensorMesh([4])
|
||||
Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True)
|
||||
x = np.arange(4)+1
|
||||
self.assertTrue(np.linalg.norm(Q*x - np.r_[1,1.004,1.008]) < TOL)
|
||||
Q = M.getInterpolationMat(np.array([[-1],[0.126],[0.127]]),'CC',zerosOutside=True)
|
||||
self.assertTrue(np.linalg.norm(Q*x - np.r_[0,1.004,1.008]) < TOL)
|
||||
|
||||
class TestInterpolation2d(Tests.OrderTest):
|
||||
name = "Interpolation 2D"
|
||||
LOCS = np.random.rand(50,2)*0.6+0.2
|
||||
meshTypes = MESHTYPES
|
||||
tolerance = TOLERANCES
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
funX = lambda x, y: np.cos(2*np.pi*y)
|
||||
funY = lambda x, y: np.cos(2*np.pi*x)
|
||||
|
||||
if 'x' in self.type:
|
||||
ana = call2(funX, self.LOCS)
|
||||
elif 'y' in self.type:
|
||||
ana = call2(funY, self.LOCS)
|
||||
else:
|
||||
ana = call2(funX, self.LOCS)
|
||||
|
||||
if 'F' in self.type:
|
||||
Fc = cartF2(self.M, funX, funY)
|
||||
grid = self.M.projectFaceVector(Fc)
|
||||
elif 'E' in self.type:
|
||||
Ec = cartE2(self.M, funX, funY)
|
||||
grid = self.M.projectEdgeVector(Ec)
|
||||
elif 'CC' == self.type:
|
||||
grid = call2(funX, self.M.gridCC)
|
||||
elif 'N' == self.type:
|
||||
grid = call2(funX, self.M.gridN)
|
||||
|
||||
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
|
||||
|
||||
err = np.linalg.norm((comp - ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderCC(self):
|
||||
self.type = 'CC'
|
||||
self.name = 'Interpolation 2D: CC'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN(self):
|
||||
self.type = 'N'
|
||||
self.name = 'Interpolation 2D: N'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFx(self):
|
||||
self.type = 'Fx'
|
||||
self.name = 'Interpolation 2D: Fx'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFy(self):
|
||||
self.type = 'Fy'
|
||||
self.name = 'Interpolation 2D: Fy'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEx(self):
|
||||
self.type = 'Ex'
|
||||
self.name = 'Interpolation 2D: Ex'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEy(self):
|
||||
self.type = 'Ey'
|
||||
self.name = 'Interpolation 2D: Ey'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestInterpolation2dCyl_Simple(unittest.TestCase):
|
||||
def test_simpleInter(self):
|
||||
M = Mesh.CylMesh([4,1,1])
|
||||
locs = np.r_[0,0,0.5]
|
||||
fx = np.array([[ 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
|
||||
self.assertTrue( np.all(fx == M.getInterpolationMat(locs, 'Fx').todense()) )
|
||||
fz = np.array([[ 0., 0., 0., 0., 0.5, 0., 0., 0., 0.5, 0., 0., 0.]])
|
||||
self.assertTrue( np.all(fz == M.getInterpolationMat(locs, 'Fz').todense()) )
|
||||
|
||||
def test_exceptions(self):
|
||||
M = Mesh.CylMesh([4,1,1])
|
||||
locs = np.r_[0,0,0.5]
|
||||
self.assertRaises(Exception,lambda:M.getInterpolationMat(locs, 'Fy'))
|
||||
self.assertRaises(Exception,lambda:M.getInterpolationMat(locs, 'Ex'))
|
||||
self.assertRaises(Exception,lambda:M.getInterpolationMat(locs, 'Ez'))
|
||||
|
||||
|
||||
class TestInterpolation2dCyl(Tests.OrderTest):
|
||||
name = "Interpolation 2D"
|
||||
LOCS = np.c_[np.random.rand(4)*0.6+0.2, np.zeros(4), np.random.rand(4)*0.6+0.2]
|
||||
meshTypes = ['uniformCylMesh'] # MESHTYPES +
|
||||
tolerance = 0.6
|
||||
meshDimension = 2
|
||||
meshSizes = [32, 64, 128, 256]
|
||||
|
||||
def getError(self):
|
||||
funX = lambda x, y: np.cos(2*np.pi*y)
|
||||
funY = lambda x, y: np.cos(2*np.pi*x)
|
||||
|
||||
if 'x' in self.type:
|
||||
ana = call2(funX, self.LOCS)
|
||||
elif 'y' in self.type:
|
||||
ana = call2(funY, self.LOCS)
|
||||
elif 'z' in self.type:
|
||||
ana = call2(funY, self.LOCS)
|
||||
else:
|
||||
ana = call2(funX, self.LOCS)
|
||||
|
||||
if 'Fx' == self.type:
|
||||
Fc = cartF2Cyl(self.M, funX, funY)
|
||||
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
|
||||
grid = self.M.projectFaceVector(Fc)
|
||||
elif 'Fz' == self.type:
|
||||
Fc = cartF2Cyl(self.M, funX, funY)
|
||||
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
|
||||
|
||||
grid = self.M.projectFaceVector(Fc)
|
||||
elif 'E' in self.type:
|
||||
Ec = cartE2Cyl(self.M, funX, funY)
|
||||
grid = Ec[:,1]
|
||||
elif 'CC' == self.type:
|
||||
grid = call2(funX, self.M.gridCC)
|
||||
elif 'N' == self.type:
|
||||
grid = call2(funX, self.M.gridN)
|
||||
|
||||
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
|
||||
|
||||
err = np.linalg.norm((comp - ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderCC(self):
|
||||
self.type = 'CC'
|
||||
self.name = 'Interpolation 2D CYLMESH: CC'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN(self):
|
||||
self.type = 'N'
|
||||
self.name = 'Interpolation 2D CYLMESH: N'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFx(self):
|
||||
self.type = 'Fx'
|
||||
self.name = 'Interpolation 2D CYLMESH: Fx'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFz(self):
|
||||
self.type = 'Fz'
|
||||
self.name = 'Interpolation 2D CYLMESH: Fz'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEy(self):
|
||||
self.type = 'Ey'
|
||||
self.name = 'Interpolation 2D CYLMESH: Ey'
|
||||
self.orderTest()
|
||||
|
||||
class TestInterpolation3D(Tests.OrderTest):
|
||||
name = "Interpolation"
|
||||
LOCS = np.random.rand(50,3)*0.6+0.2
|
||||
meshTypes = MESHTYPES
|
||||
tolerance = TOLERANCES
|
||||
meshDimension = 3
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
funX = lambda x, y, z: np.cos(2*np.pi*y)
|
||||
funY = lambda x, y, z: np.cos(2*np.pi*z)
|
||||
funZ = lambda x, y, z: np.cos(2*np.pi*x)
|
||||
|
||||
if 'x' in self.type:
|
||||
ana = call3(funX, self.LOCS)
|
||||
elif 'y' in self.type:
|
||||
ana = call3(funY, self.LOCS)
|
||||
elif 'z' in self.type:
|
||||
ana = call3(funZ, self.LOCS)
|
||||
else:
|
||||
ana = call3(funX, self.LOCS)
|
||||
|
||||
if 'F' in self.type:
|
||||
Fc = cartF3(self.M, funX, funY, funZ)
|
||||
grid = self.M.projectFaceVector(Fc)
|
||||
elif 'E' in self.type:
|
||||
Ec = cartE3(self.M, funX, funY, funZ)
|
||||
grid = self.M.projectEdgeVector(Ec)
|
||||
elif 'CC' == self.type:
|
||||
grid = call3(funX, self.M.gridCC)
|
||||
elif 'N' == self.type:
|
||||
grid = call3(funX, self.M.gridN)
|
||||
|
||||
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
|
||||
|
||||
err = np.linalg.norm((comp - ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderCC(self):
|
||||
self.type = 'CC'
|
||||
self.name = 'Interpolation 3D: CC'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN(self):
|
||||
self.type = 'N'
|
||||
self.name = 'Interpolation 3D: N'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFx(self):
|
||||
self.type = 'Fx'
|
||||
self.name = 'Interpolation 3D: Fx'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFy(self):
|
||||
self.type = 'Fy'
|
||||
self.name = 'Interpolation 3D: Fy'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFz(self):
|
||||
self.type = 'Fz'
|
||||
self.name = 'Interpolation 3D: Fz'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEx(self):
|
||||
self.type = 'Ex'
|
||||
self.name = 'Interpolation 3D: Ex'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEy(self):
|
||||
self.type = 'Ey'
|
||||
self.name = 'Interpolation 3D: Ey'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEz(self):
|
||||
self.type = 'Ez'
|
||||
self.name = 'Interpolation 3D: Ez'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,483 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG.Tests import OrderTest
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
#TODO: 'randomTensorMesh'
|
||||
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)]
|
||||
cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
|
||||
cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
|
||||
cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
|
||||
cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
|
||||
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
|
||||
|
||||
|
||||
class TestCurl(OrderTest):
|
||||
name = "Curl"
|
||||
meshTypes = MESHTYPES
|
||||
|
||||
def getError(self):
|
||||
# fun: i (cos(y)) + j (cos(z)) + k (cos(x))
|
||||
# sol: i (sin(z)) + j (sin(x)) + k (sin(y))
|
||||
|
||||
funX = lambda x, y, z: np.cos(2*np.pi*y)
|
||||
funY = lambda x, y, z: np.cos(2*np.pi*z)
|
||||
funZ = lambda x, y, z: np.cos(2*np.pi*x)
|
||||
|
||||
solX = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*z)
|
||||
solY = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)
|
||||
solZ = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*y)
|
||||
|
||||
Ec = cartE3(self.M, funX, funY, funZ)
|
||||
E = self.M.projectEdgeVector(Ec)
|
||||
|
||||
Fc = cartF3(self.M, solX, solY, solZ)
|
||||
curlE_ana = self.M.projectFaceVector(Fc)
|
||||
|
||||
curlE = self.M.edgeCurl.dot(E)
|
||||
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)
|
||||
else:
|
||||
err = np.linalg.norm((curlE - curlE_ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestCurl2D(OrderTest):
|
||||
name = "Cell Grad 2D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
ex = lambda x, y: np.cos(y)
|
||||
ey = lambda x, y: np.cos(x)
|
||||
sol = lambda x, y: -np.sin(x)+np.sin(y)
|
||||
|
||||
sol_curl2d = call2(sol, self.M.gridCC)
|
||||
Ec = cartE2(self.M, ex, ey)
|
||||
sol_ana = self.M.edgeCurl*self.M.projectFaceVector(Ec)
|
||||
err = np.linalg.norm((sol_curl2d-sol_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
|
||||
name = "Cell Grad 1D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 1
|
||||
expectedOrders = 1 # because of the averaging involved in the ghost point. u_b = (u_n + u_g)/2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x: -2*np.pi*np.sin(2*np.pi*x)
|
||||
sol = lambda x: np.cos(2*np.pi*x)
|
||||
|
||||
|
||||
xc = sol(self.M.gridCC)
|
||||
|
||||
gradX_ana = fx(self.M.gridFx)
|
||||
|
||||
bc = np.array([1,1])
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
gradX = self.M.cellGrad.dot(xc) + self.M.cellGradBC*bc
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestCellGrad2D_Dirichlet(OrderTest):
|
||||
name = "Cell Grad 2D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)
|
||||
fy = lambda x, y: 2*np.pi*np.cos(2*np.pi*y)*np.sin(2*np.pi*x)
|
||||
sol = lambda x, y: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)
|
||||
|
||||
xc = call2(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF2(self.M, fx, fy)
|
||||
gradX_ana = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestCellGrad3D_Dirichlet(OrderTest):
|
||||
name = "Cell Grad 3D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 3
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y, z: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
fy = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
fz = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
sol = lambda x, y, z: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
|
||||
xc = call3(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF3(self.M, fx, fy, fz)
|
||||
gradX_ana = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestCellGrad2D_Neumann(OrderTest):
|
||||
name = "Cell Grad 2D - Neumann"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)
|
||||
fy = lambda x, y: -2*np.pi*np.sin(2*np.pi*y)*np.cos(2*np.pi*x)
|
||||
sol = lambda x, y: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)
|
||||
|
||||
xc = call2(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF2(self.M, fx, fy)
|
||||
gradX_ana = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('neumann')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestCellGrad3D_Neumann(OrderTest):
|
||||
name = "Cell Grad 3D - Neumann"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 3
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y, z: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
fy = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
fz = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
sol = lambda x, y, z: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
|
||||
xc = call3(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF3(self.M, fx, fy, fz)
|
||||
gradX_ana = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('neumann')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestFaceDiv3D(OrderTest):
|
||||
name = "Face Divergence 3D"
|
||||
meshTypes = MESHTYPES
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y, z: np.sin(2*np.pi*x)
|
||||
fy = lambda x, y, z: np.sin(2*np.pi*y)
|
||||
fz = lambda x, y, z: np.sin(2*np.pi*z)
|
||||
sol = lambda x, y, z: (2*np.pi*np.cos(2*np.pi*x)+2*np.pi*np.cos(2*np.pi*y)+2*np.pi*np.cos(2*np.pi*z))
|
||||
|
||||
Fc = cartF3(self.M, fx, fy, fz)
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
divF = self.M.faceDiv.dot(F)
|
||||
divF_ana = call3(sol, self.M.gridCC)
|
||||
|
||||
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)
|
||||
else:
|
||||
err = np.linalg.norm((divF-divF_ana), np.inf)
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestFaceDiv2D(OrderTest):
|
||||
name = "Face Divergence 2D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y: np.sin(2*np.pi*x)
|
||||
fy = lambda x, y: np.sin(2*np.pi*y)
|
||||
sol = lambda x, y: 2*np.pi*(np.cos(2*np.pi*x)+np.cos(2*np.pi*y))
|
||||
|
||||
Fc = cartF2(self.M, fx, fy)
|
||||
F = self.M.projectFaceVector(Fc)
|
||||
|
||||
divF = self.M.faceDiv.dot(F)
|
||||
divF_ana = call2(sol, self.M.gridCC)
|
||||
|
||||
err = np.linalg.norm((divF-divF_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad(OrderTest):
|
||||
name = "Nodal Gradient"
|
||||
meshTypes = MESHTYPES
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))
|
||||
# i (sin(x)) + j (sin(y)) + k (sin(z))
|
||||
solX = lambda x, y, z: -np.sin(x)
|
||||
solY = lambda x, y, z: -np.sin(y)
|
||||
solZ = lambda x, y, z: -np.sin(z)
|
||||
|
||||
phi = call3(fun, self.M.gridN)
|
||||
gradE = self.M.nodalGrad.dot(phi)
|
||||
|
||||
Ec = cartE3(self.M, solX, solY, solZ)
|
||||
gradE_ana = self.M.projectEdgeVector(Ec)
|
||||
|
||||
err = np.linalg.norm((gradE-gradE_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad2D(OrderTest):
|
||||
name = "Nodal Gradient 2D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fun = lambda x, y: (np.cos(x)+np.cos(y))
|
||||
# i (sin(x)) + j (sin(y)) + k (sin(z))
|
||||
solX = lambda x, y: -np.sin(x)
|
||||
solY = lambda x, y: -np.sin(y)
|
||||
|
||||
phi = call2(fun, self.M.gridN)
|
||||
gradE = self.M.nodalGrad.dot(phi)
|
||||
|
||||
Ec = cartE2(self.M, solX, solY)
|
||||
gradE_ana = self.M.projectEdgeVector(Ec)
|
||||
|
||||
err = np.linalg.norm((gradE-gradE_ana), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestAveraging2D(OrderTest):
|
||||
name = "Averaging 2D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
|
||||
def getError(self):
|
||||
num = self.getAve(self.M) * self.getHere(self.M)
|
||||
err = np.linalg.norm((self.getThere(self.M)-num), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderN2CC(self):
|
||||
self.name = "Averaging 2D: N2CC"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: call2(fun, M.gridN)
|
||||
self.getThere = lambda M: call2(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveN2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN2F(self):
|
||||
self.name = "Averaging 2D: N2F"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: call2(fun, M.gridN)
|
||||
self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
|
||||
self.getAve = lambda M: M.aveN2F
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN2E(self):
|
||||
self.name = "Averaging 2D: N2E"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: call2(fun, M.gridN)
|
||||
self.getThere = lambda M: np.r_[call2(fun, M.gridEx), call2(fun, M.gridEy)]
|
||||
self.getAve = lambda M: M.aveN2E
|
||||
self.orderTest()
|
||||
|
||||
def test_orderF2CC(self):
|
||||
self.name = "Averaging 2D: F2CC"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
|
||||
self.getThere = lambda M: call2(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveF2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderF2CCV(self):
|
||||
self.name = "Averaging 2D: F2CCV"
|
||||
funX = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
funY = lambda x, y: (np.cos(y)*np.sin(x))
|
||||
self.getHere = lambda M: np.r_[call2(funX, M.gridFx), call2(funY, M.gridFy)]
|
||||
self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveF2CCV
|
||||
self.orderTest()
|
||||
|
||||
def test_orderCC2F(self):
|
||||
self.name = "Averaging 2D: CC2F"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: call2(fun, M.gridCC)
|
||||
self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
|
||||
self.getAve = lambda M: M.aveCC2F
|
||||
self.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
def test_orderE2CC(self):
|
||||
self.name = "Averaging 2D: E2CC"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: np.r_[call2(fun, M.gridEx), call2(fun, M.gridEy)]
|
||||
self.getThere = lambda M: call2(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveE2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderE2CCV(self):
|
||||
self.name = "Averaging 2D: E2CCV"
|
||||
funX = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
funY = lambda x, y: (np.cos(y)*np.sin(x))
|
||||
self.getHere = lambda M: np.r_[call2(funX, M.gridEx), call2(funY, M.gridEy)]
|
||||
self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveE2CCV
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestAveraging3D(OrderTest):
|
||||
name = "Averaging 3D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 3
|
||||
|
||||
def getError(self):
|
||||
num = self.getAve(self.M) * self.getHere(self.M)
|
||||
err = np.linalg.norm((self.getThere(self.M)-num), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderN2CC(self):
|
||||
self.name = "Averaging 3D: N2CC"
|
||||
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: call3(fun, M.gridN)
|
||||
self.getThere = lambda M: call3(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveN2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN2F(self):
|
||||
self.name = "Averaging 3D: N2F"
|
||||
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: call3(fun, M.gridN)
|
||||
self.getThere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
|
||||
self.getAve = lambda M: M.aveN2F
|
||||
self.orderTest()
|
||||
|
||||
def test_orderN2E(self):
|
||||
self.name = "Averaging 3D: N2E"
|
||||
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: call3(fun, M.gridN)
|
||||
self.getThere = lambda M: np.r_[call3(fun, M.gridEx), call3(fun, M.gridEy), call3(fun, M.gridEz)]
|
||||
self.getAve = lambda M: M.aveN2E
|
||||
self.orderTest()
|
||||
|
||||
def test_orderF2CC(self):
|
||||
self.name = "Averaging 3D: F2CC"
|
||||
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
|
||||
self.getThere = lambda M: call3(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveF2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderF2CCV(self):
|
||||
self.name = "Averaging 3D: F2CCV"
|
||||
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
|
||||
funZ = lambda x, y, z: (np.cos(x)*np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funX, M.gridFx), call3(funY, M.gridFy), call3(funZ, M.gridFz)]
|
||||
self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveF2CCV
|
||||
self.orderTest()
|
||||
|
||||
def test_orderE2CC(self):
|
||||
self.name = "Averaging 3D: E2CC"
|
||||
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(fun, M.gridEx), call3(fun, M.gridEy), call3(fun, M.gridEz)]
|
||||
self.getThere = lambda M: call3(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveE2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderE2CCV(self):
|
||||
self.name = "Averaging 3D: E2CCV"
|
||||
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
|
||||
funZ = lambda x, y, z: (np.cos(x)*np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funX, M.gridEx), call3(funY, M.gridEy), call3(funZ, M.gridEz)]
|
||||
self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveE2CCV
|
||||
self.orderTest()
|
||||
|
||||
def test_orderCC2F(self):
|
||||
self.name = "Averaging 3D: CC2F"
|
||||
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: call3(fun, M.gridCC)
|
||||
self.getThere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
|
||||
self.getAve = lambda M: M.aveCC2F
|
||||
self.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,128 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG.Mesh import TensorMesh
|
||||
from SimPEG import Solver, Tests
|
||||
|
||||
TOL = 1e-10
|
||||
|
||||
class BasicTensorMeshTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
a = np.array([1, 1, 1])
|
||||
b = np.array([1, 2])
|
||||
c = np.array([1, 4])
|
||||
self.mesh2 = TensorMesh([a, b], [3, 5])
|
||||
self.mesh3 = TensorMesh([a, b, c])
|
||||
|
||||
def test_vectorN_2D(self):
|
||||
testNx = np.array([3, 4, 5, 6])
|
||||
testNy = np.array([5, 6, 8])
|
||||
|
||||
xtest = np.all(self.mesh2.vectorNx == testNx)
|
||||
ytest = np.all(self.mesh2.vectorNy == testNy)
|
||||
self.assertTrue(xtest and ytest)
|
||||
|
||||
def test_vectorCC_2D(self):
|
||||
testNx = np.array([3.5, 4.5, 5.5])
|
||||
testNy = np.array([5.5, 7])
|
||||
|
||||
xtest = np.all(self.mesh2.vectorCCx == testNx)
|
||||
ytest = np.all(self.mesh2.vectorCCy == testNy)
|
||||
self.assertTrue(xtest and ytest)
|
||||
|
||||
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])
|
||||
t1 = np.all(self.mesh3.area == test_area)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_vol_3D(self):
|
||||
test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
|
||||
t1 = np.all(self.mesh3.vol == test_vol)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_vol_2D(self):
|
||||
test_vol = np.array([1, 1, 1, 2, 2, 2])
|
||||
t1 = np.all(self.mesh2.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.mesh3.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.mesh2.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_oneCell(self):
|
||||
hx = np.array([1e-5])
|
||||
M = TensorMesh([hx])
|
||||
self.assertTrue(M.nC == 1)
|
||||
|
||||
def test_printing(self):
|
||||
print TensorMesh([10])
|
||||
print TensorMesh([10,10])
|
||||
print TensorMesh([10,10,10])
|
||||
|
||||
def test_centering(self):
|
||||
M1d = TensorMesh([10], 'C')
|
||||
M2d = TensorMesh([10,10], 'CC')
|
||||
M3d = TensorMesh([10,10,10], 'CCC')
|
||||
self.assertLess(np.abs(M1d.x0 + 0.5).sum(), TOL)
|
||||
self.assertLess(np.abs(M2d.x0 + 0.5).sum(), TOL)
|
||||
self.assertLess(np.abs(M3d.x0 + 0.5).sum(), TOL)
|
||||
|
||||
def test_negative(self):
|
||||
M1d = TensorMesh([10], 'N')
|
||||
self.assertRaises(Exception, TensorMesh, [10], 'F')
|
||||
M2d = TensorMesh([10,10], 'NN')
|
||||
M3d = TensorMesh([10,10,10], 'NNN')
|
||||
self.assertLess(np.abs(M1d.x0 + 1.0).sum(), TOL)
|
||||
self.assertLess(np.abs(M2d.x0 + 1.0).sum(), TOL)
|
||||
self.assertLess(np.abs(M3d.x0 + 1.0).sum(), TOL)
|
||||
|
||||
def test_cent_neg(self):
|
||||
M3d = TensorMesh([10,10,10], 'C0N')
|
||||
self.assertLess(np.abs(M3d.x0 + np.r_[0.5,0,1.0]).sum(), TOL)
|
||||
|
||||
def test_tensor(self):
|
||||
M = TensorMesh([[(10.,2)]])
|
||||
self.assertLess(np.abs(M.hx - np.r_[10.,10.]).sum(), TOL)
|
||||
|
||||
class TestPoissonEqn(Tests.OrderTest):
|
||||
name = "Poisson Equation"
|
||||
meshSizes = [10, 16, 20]
|
||||
|
||||
def getError(self):
|
||||
# Create some functions to integrate
|
||||
fun = lambda x: np.sin(2*np.pi*x[:, 0])*np.sin(2*np.pi*x[:, 1])*np.sin(2*np.pi*x[:, 2])
|
||||
sol = lambda x: -3.*((2*np.pi)**2)*fun(x)
|
||||
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
|
||||
D = self.M.faceDiv
|
||||
G = self.M.cellGrad
|
||||
if self.forward:
|
||||
sA = sol(self.M.gridCC)
|
||||
sN = D*G*fun(self.M.gridCC)
|
||||
err = np.linalg.norm((sA - sN), np.inf)
|
||||
else:
|
||||
fA = fun(self.M.gridCC)
|
||||
fN = Solver(D*G) * (sol(self.M.gridCC))
|
||||
err = np.linalg.norm((fA - fN), np.inf)
|
||||
return err
|
||||
|
||||
def test_orderForward(self):
|
||||
self.name = "Poisson Equation - Forward"
|
||||
self.forward = True
|
||||
self.orderTest()
|
||||
|
||||
def test_orderBackward(self):
|
||||
self.name = "Poisson Equation - Backward"
|
||||
self.forward = False
|
||||
self.orderTest()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user