mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-12 12:30:37 +08:00
Merge branch 'dev' of https://github.com/simpeg/simpeg into mt/dev
Conflicts: SimPEG/Utils/meshutils.py docs/index.rst tests/mesh/test_MeshIO.py
This commit is contained in:
Rowan Cockett
@@ -4,11 +4,17 @@ from SimPEG import *
|
||||
from scipy.sparse.linalg import dsolve
|
||||
import inspect
|
||||
|
||||
TOL = 1e-20
|
||||
|
||||
class RegularizationTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.mesh2 = Mesh.TensorMesh([3, 2])
|
||||
hx, hy, hz = np.random.rand(10), np.random.rand(9), np.random.rand(8)
|
||||
hx, hy, hz = hx/hx.sum(), hy/hy.sum(), hz/hz.sum()
|
||||
mesh1 = Mesh.TensorMesh([hx])
|
||||
mesh2 = Mesh.TensorMesh([hx, hy])
|
||||
mesh3 = Mesh.TensorMesh([hx, hy, hz])
|
||||
self.meshlist = [mesh1,mesh2, mesh3]
|
||||
|
||||
def test_regularization(self):
|
||||
for R in dir(Regularization):
|
||||
@@ -16,18 +22,63 @@ class RegularizationTests(unittest.TestCase):
|
||||
if not inspect.isclass(r): continue
|
||||
if not issubclass(r, Regularization.BaseRegularization):
|
||||
continue
|
||||
# if 'Regularization' not in R: continue
|
||||
mapping = r.mapPair(self.mesh2)
|
||||
reg = r(self.mesh2, mapping=mapping)
|
||||
m = np.random.rand(mapping.nP)
|
||||
reg.mref = m[:]*np.mean(m)
|
||||
|
||||
print 'Check:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
print 'Check 2 Deriv:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
print 'Testing %iD'%mesh.dim
|
||||
|
||||
mapping = r.mapPair(mesh)
|
||||
reg = r(mesh, mapping=mapping)
|
||||
m = np.random.rand(mapping.nP)
|
||||
reg.mref = np.ones_like(m)*np.mean(m)
|
||||
|
||||
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
|
||||
passed = reg.eval(reg.mref) < TOL
|
||||
self.assertTrue(passed)
|
||||
|
||||
print 'Check:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
print 'Check 2 Deriv:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_regularization_ActiveCells(self):
|
||||
for R in dir(Regularization):
|
||||
r = getattr(Regularization, R)
|
||||
if not inspect.isclass(r): continue
|
||||
if not issubclass(r, Regularization.BaseRegularization):
|
||||
continue
|
||||
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
print 'Testing Active Cells %iD'%(mesh.dim)
|
||||
|
||||
if mesh.dim == 1:
|
||||
indAct = Utils.mkvc(mesh.gridCC <= 0.8)
|
||||
elif mesh.dim == 2:
|
||||
indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
|
||||
elif mesh.dim == 3:
|
||||
indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
|
||||
|
||||
mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
|
||||
|
||||
reg = r(mesh, mapping=mapping, indActive=indAct)
|
||||
m = np.random.rand(mesh.nC)[indAct]
|
||||
reg.mref = np.ones_like(m)*np.mean(m)
|
||||
|
||||
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
|
||||
passed = reg.eval(reg.mref) < TOL
|
||||
self.assertTrue(passed)
|
||||
|
||||
print 'Check:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
print 'Check 2 Deriv:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
+61
-11
@@ -4,11 +4,65 @@ import SimPEG as simpeg
|
||||
from SimPEG.Mesh import TensorMesh, TreeMesh
|
||||
|
||||
|
||||
class TestOcTreeIO(unittest.TestCase):
|
||||
class TestTensorMeshIO(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
h = np.ones(16)
|
||||
mesh = simpeg.Mesh.TreeMesh([h,2*h,3*h])
|
||||
mesh = TensorMesh([h,2*h,3*h])
|
||||
self.mesh = mesh
|
||||
|
||||
def test_UBCfiles(self):
|
||||
|
||||
mesh = self.mesh
|
||||
# Make a vector
|
||||
vec = np.arange(mesh.nC)
|
||||
# Write and read
|
||||
mesh.writeUBC('temp.msh', {'arange.txt':vec})
|
||||
meshUBC = TensorMesh.readUBC('temp.msh')
|
||||
vecUBC = meshUBC.readModelUBC('arange.txt')
|
||||
|
||||
# The mesh
|
||||
assert mesh.__str__() == meshUBC.__str__()
|
||||
assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0
|
||||
assert np.sum(vec - vecUBC) == 0
|
||||
assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0)
|
||||
|
||||
|
||||
vecUBC = mesh.readModelUBC('arange.txt')
|
||||
assert np.sum(vec - vecUBC) == 0
|
||||
|
||||
mesh.writeModelUBC('arange2.txt', vec + 1)
|
||||
vec2UBC = mesh.readModelUBC('arange2.txt')
|
||||
assert np.sum(vec + 1 - vec2UBC) == 0
|
||||
|
||||
print 'IO of UBC tensor mesh files is working'
|
||||
os.remove('temp.msh')
|
||||
os.remove('arange.txt')
|
||||
os.remove('arange2.txt')
|
||||
|
||||
def test_VTKfiles(self):
|
||||
mesh = self.mesh
|
||||
vec = np.arange(mesh.nC)
|
||||
|
||||
mesh.writeVTK('temp.vtr', {'arange.txt':vec})
|
||||
meshVTR, models = TensorMesh.readVTK('temp.vtr')
|
||||
|
||||
assert mesh.__str__() == meshVTR.__str__()
|
||||
assert np.all(np.array(mesh.h) - np.array(meshVTR.h) == 0)
|
||||
|
||||
assert 'arange.txt' in models
|
||||
vecVTK = models['arange.txt']
|
||||
assert np.sum(vec - vecVTK) == 0
|
||||
|
||||
print 'IO of VTR tensor mesh files is working'
|
||||
os.remove('temp.vtr')
|
||||
|
||||
|
||||
class TestOcTreeMeshIO(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
h = np.ones(16)
|
||||
mesh = TreeMesh([h,2*h,3*h])
|
||||
mesh.refine(3)
|
||||
mesh._refineCell([0,0,0,3])
|
||||
mesh._refineCell([0,2,0,3])
|
||||
@@ -19,9 +73,10 @@ class TestOcTreeIO(unittest.TestCase):
|
||||
mesh = self.mesh
|
||||
# Make a vector
|
||||
vec = np.arange(mesh.nC)
|
||||
# Write aand read
|
||||
simpeg.Utils.meshutils.writeUBCocTreeFiles('temp.msh',mesh,{'arange.txt':vec})
|
||||
meshUBC, vecUBC = simpeg.Utils.meshutils.readUBCocTreeFiles('temp.msh',['arange.txt'])
|
||||
# Write and read
|
||||
mesh.writeUBC('temp.msh', {'arange.txt':vec})
|
||||
meshUBC = TreeMesh.readUBC('temp.msh')
|
||||
vecUBC = meshUBC.readModelUBC('arange.txt')
|
||||
|
||||
# The mesh
|
||||
assert mesh.__str__() == meshUBC.__str__()
|
||||
@@ -35,12 +90,7 @@ class TestOcTreeIO(unittest.TestCase):
|
||||
def test_VTUfiles(self):
|
||||
mesh = self.mesh
|
||||
vec = np.arange(mesh.nC)
|
||||
try:
|
||||
simpeg.Utils.meshutils.writeVTUFile('temp.vtu',mesh,{'arange':vec})
|
||||
run = True
|
||||
except:
|
||||
run = False
|
||||
assert run
|
||||
mesh.writeVTK('temp.vtu',{'arange':vec})
|
||||
print 'Writing of VTU files is working'
|
||||
os.remove('temp.vtu')
|
||||
|
||||
|
||||
@@ -26,6 +26,27 @@ class TestSimpleQuadTree(unittest.TestCase):
|
||||
|
||||
assert np.allclose(np.r_[M._areaFxFull, M._areaFyFull], M._deflationMatrix('F') * M.area)
|
||||
|
||||
def test_getitem(self):
|
||||
M = Mesh.TreeMesh([4,4])
|
||||
M.refine(1)
|
||||
assert M.nC == 4
|
||||
assert len(M) == M.nC
|
||||
assert np.allclose(M[0].center, [0.25,0.25])
|
||||
actual = [[0,0],[0.5,0],[0,0.5],[0.5,0.5]]
|
||||
for i, n in enumerate(M[0].nodes):
|
||||
assert np.allclose(M._gridN[n,:], actual[i])
|
||||
|
||||
def test_getitem3D(self):
|
||||
M = Mesh.TreeMesh([4,4,4])
|
||||
M.refine(1)
|
||||
assert M.nC == 8
|
||||
assert len(M) == M.nC
|
||||
assert np.allclose(M[0].center, [0.25,0.25,0.25])
|
||||
actual = [[0,0,0],[0.5,0,0],[0,0.5,0],[0.5,0.5,0],
|
||||
[0,0,0.5],[0.5,0,0.5],[0,0.5,0.5],[0.5,0.5,0.5]]
|
||||
for i, n in enumerate(M[0].nodes):
|
||||
assert np.allclose(M._gridN[n,:], actual[i])
|
||||
|
||||
def test_refine(self):
|
||||
M = Mesh.TreeMesh([4,4,4])
|
||||
M.refine(1)
|
||||
|
||||
Reference in New Issue
Block a user