mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-15 11:26:09 +08:00
Rowan Cockett
@@ -38,5 +38,4 @@ nosetests.xml
|
||||
*.sublime-project
|
||||
*.sublime-workspace
|
||||
docs/_build/
|
||||
*_cython.c
|
||||
Makefile
|
||||
|
||||
@@ -30,6 +30,7 @@ install:
|
||||
- pip install nose-cov python-coveralls
|
||||
# - pip install -r requirements.txt
|
||||
- python setup.py install
|
||||
- python setup.py build_ext --inplace
|
||||
|
||||
# Run test
|
||||
script:
|
||||
|
||||
@@ -27,6 +27,7 @@ class BaseMesh(object):
|
||||
# Ensure x0 & n are 1D vectors
|
||||
self._n = np.array(n, dtype=int).ravel()
|
||||
self._x0 = np.array(x0, dtype=float).ravel()
|
||||
self._dim = len(self._x0)
|
||||
|
||||
@property
|
||||
def x0(self):
|
||||
@@ -46,7 +47,7 @@ class BaseMesh(object):
|
||||
:rtype: int
|
||||
:return: dim
|
||||
"""
|
||||
return len(self._n)
|
||||
return self._dim
|
||||
|
||||
@property
|
||||
def nC(self):
|
||||
|
||||
@@ -2,12 +2,12 @@ import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from scipy.constants import pi
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
|
||||
from TensorMesh import BaseTensorMesh
|
||||
from TensorMesh import BaseTensorMesh, BaseRectangularMesh
|
||||
from InnerProducts import InnerProducts
|
||||
from View import CylView
|
||||
|
||||
|
||||
class CylMesh(BaseTensorMesh, InnerProducts, CylView):
|
||||
class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
|
||||
"""
|
||||
CylMesh is a mesh class for cylindrical problems
|
||||
|
||||
|
||||
+27
-24
@@ -1,10 +1,10 @@
|
||||
from SimPEG import Utils, np, sp
|
||||
from BaseMesh import BaseRectangularMesh
|
||||
from BaseMesh import BaseMesh, BaseRectangularMesh
|
||||
from View import TensorView
|
||||
from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
|
||||
class BaseTensorMesh(BaseRectangularMesh):
|
||||
class BaseTensorMesh(BaseMesh):
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
@@ -42,7 +42,10 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
else:
|
||||
raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i)
|
||||
|
||||
BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
|
||||
if isinstance(self, BaseRectangularMesh):
|
||||
BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
|
||||
else:
|
||||
BaseMesh.__init__(self, np.array([x.size for x in h]), x0)
|
||||
|
||||
# Ensure h contains 1D vectors
|
||||
self._h = [Utils.mkvc(x.astype(float)) for x in h]
|
||||
@@ -356,7 +359,7 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
|
||||
|
||||
|
||||
class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts):
|
||||
class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
|
||||
"""
|
||||
TensorMesh is a mesh class that deals with tensor product meshes.
|
||||
|
||||
@@ -413,34 +416,34 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts):
|
||||
break
|
||||
|
||||
if n == 1:
|
||||
outStr = outStr + ' {0:.2f},'.format(h)
|
||||
outStr += ' {0:.2f},'.format(h)
|
||||
else:
|
||||
outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h)
|
||||
outStr += ' {0:d}*{1:.2f},'.format(n,h)
|
||||
|
||||
return outStr[:-1]
|
||||
|
||||
if self.dim == 1:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr += '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr += '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr += printH(self.hx, outStr='\n hx:')
|
||||
pass
|
||||
elif self.dim == 2:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr = outStr + printH(self.hy, outStr='\n hy:')
|
||||
outStr += '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr += '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr += '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr += '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr += printH(self.hx, outStr='\n hx:')
|
||||
outStr += printH(self.hy, outStr='\n hy:')
|
||||
elif self.dim == 3:
|
||||
outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2])
|
||||
outStr = outStr + '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr = outStr + '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr = outStr + '\n nCz: {0:d}'.format(self.nCz)
|
||||
outStr = outStr + printH(self.hx, outStr='\n hx:')
|
||||
outStr = outStr + printH(self.hy, outStr='\n hy:')
|
||||
outStr = outStr + printH(self.hz, outStr='\n hz:')
|
||||
outStr += '\n x0: {0:.2f}'.format(self.x0[0])
|
||||
outStr += '\n y0: {0:.2f}'.format(self.x0[1])
|
||||
outStr += '\n z0: {0:.2f}'.format(self.x0[2])
|
||||
outStr += '\n nCx: {0:d}'.format(self.nCx)
|
||||
outStr += '\n nCy: {0:d}'.format(self.nCy)
|
||||
outStr += '\n nCz: {0:d}'.format(self.nCz)
|
||||
outStr += printH(self.hx, outStr='\n hx:')
|
||||
outStr += printH(self.hy, outStr='\n hy:')
|
||||
outStr += printH(self.hz, outStr='\n hz:')
|
||||
|
||||
return outStr
|
||||
|
||||
|
||||
+2320
-1105
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,85 @@
|
||||
# from __future__ import division
|
||||
# import numpy as np
|
||||
# cimport numpy as np
|
||||
# from libcpp.vector cimport vector
|
||||
|
||||
|
||||
"""
|
||||
The Z-order curve is generated by interleaving the bits of an offset.
|
||||
|
||||
See also:
|
||||
|
||||
https://github.com/cortesi/scurve
|
||||
Aldo Cortesi <aldo@corte.si>
|
||||
|
||||
"""
|
||||
|
||||
def bitrange(long x, int width, int start, int end):
|
||||
"""
|
||||
Extract a bit range as an integer.
|
||||
(start, end) is inclusive lower bound, exclusive upper bound.
|
||||
"""
|
||||
return x >> (width-end) & ((2**(end-start))-1)
|
||||
|
||||
def index(int dimension, int bits, int levelBits, list p, int level):
|
||||
cdef long idx = 0
|
||||
cdef int iwidth
|
||||
cdef int i
|
||||
cdef long b
|
||||
cdef int bitoff
|
||||
|
||||
p = [_ for _ in p]
|
||||
|
||||
p.reverse()
|
||||
iwidth = bits * dimension
|
||||
for i in range(iwidth):
|
||||
bitoff = bits-(i/dimension)-1
|
||||
poff = dimension-(i%dimension)-1
|
||||
b = bitrange(p[poff], bits, bitoff, bitoff+1) << i
|
||||
idx |= b
|
||||
|
||||
return (idx << levelBits) + level
|
||||
|
||||
def point(int dimension, int bits, int levelBits, long idx):
|
||||
cdef list p
|
||||
cdef int iwidth
|
||||
cdef int i, n
|
||||
cdef long b
|
||||
|
||||
n = idx & (2**levelBits-1)
|
||||
idx = idx >> levelBits
|
||||
|
||||
p = [0]*dimension
|
||||
iwidth = bits * dimension
|
||||
for i in range(iwidth):
|
||||
b = bitrange(idx, iwidth, i, i+1) << (iwidth-i-1)/dimension
|
||||
p[i%dimension] |= b
|
||||
p.reverse()
|
||||
return p + [n]
|
||||
|
||||
|
||||
# def _refineCell(int dimension, int bits, self, pointer):
|
||||
# self._structureChange()
|
||||
# pointer = self._asPointer(pointer)
|
||||
# ind = self._asIndex(pointer)
|
||||
# assert ind in self
|
||||
# h = self._levelWidth(pointer[-1])/2 # halfWidth
|
||||
# nL = pointer[-1] + 1 # new level
|
||||
# add = lambda p:p[0]+p[1]
|
||||
# added = []
|
||||
# def addCell(p):
|
||||
# i = self._index(p+[nL])
|
||||
# self._treeInds.add(i)
|
||||
# added.append(i)
|
||||
|
||||
# addCell(map(add, zip(pointer[:-1], [0,0,0])))
|
||||
# addCell(map(add, zip(pointer[:-1], [h,0,0])))
|
||||
# addCell(map(add, zip(pointer[:-1], [0,h,0])))
|
||||
# addCell(map(add, zip(pointer[:-1], [h,h,0])))
|
||||
# if self.dim == 3:
|
||||
# addCell(map(add, zip(pointer[:-1], [0,0,h])))
|
||||
# addCell(map(add, zip(pointer[:-1], [h,0,h])))
|
||||
# addCell(map(add, zip(pointer[:-1], [0,h,h])))
|
||||
# addCell(map(add, zip(pointer[:-1], [h,h,h])))
|
||||
# self._treeInds.remove(ind)
|
||||
# return added
|
||||
+2
-1
@@ -173,7 +173,7 @@ class TensorView(object):
|
||||
ax=None, clim=None, showIt=False,
|
||||
pcolorOpts={},
|
||||
streamOpts={'color':'k'},
|
||||
gridOpts={'color':'k'}
|
||||
gridOpts={'color':'k', 'alpha':0.5}
|
||||
):
|
||||
|
||||
"""
|
||||
@@ -216,6 +216,7 @@ class TensorView(object):
|
||||
if ind is None: ind = int(szSliceDim/2)
|
||||
assert type(ind) in [int, long], 'ind must be an integer'
|
||||
|
||||
assert not (v.dtype == complex and view == 'vec'), 'Can not plot a complex vector.'
|
||||
# The slicing and plotting code!!
|
||||
|
||||
def getIndSlice(v):
|
||||
|
||||
+1
-1
@@ -367,7 +367,7 @@ class BaseSurvey(object):
|
||||
|
||||
"""
|
||||
if getattr(self, 'dobs', None) is not None and not force:
|
||||
raise Exception('Survey already has dobs.')
|
||||
raise Exception('Survey already has dobs. You can use force=True to override this exception.')
|
||||
self.mtrue = m
|
||||
self.dtrue = self.dpred(m, u=u)
|
||||
noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
|
||||
|
||||
@@ -4,6 +4,7 @@ from numpy.linalg import norm
|
||||
from SimPEG.Utils import mkvc, sdiag, diagEst
|
||||
from SimPEG import Utils
|
||||
from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
|
||||
from SimPEG.Mesh.TreeMesh import TreeMesh as Tree
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
@@ -132,6 +133,34 @@ class OrderTest(unittest.TestCase):
|
||||
self.M = CurvilinearMesh([X, Y, Z])
|
||||
return 1./nc
|
||||
|
||||
elif 'Tree' in self._meshType:
|
||||
nc *= 2
|
||||
if 'uniform' in self._meshType or 'notatree' in self._meshType:
|
||||
h = [nc, nc, nc]
|
||||
elif 'random' in self._meshType:
|
||||
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
|
||||
else:
|
||||
raise Exception('Unexpected meshType')
|
||||
|
||||
levels = int(np.log(nc)/np.log(2))
|
||||
self.M = Tree(h[:self.meshDimension], levels=levels)
|
||||
def function(cell):
|
||||
if 'notatree' in self._meshType:
|
||||
return levels - 1
|
||||
r = cell.center - np.array([0.5]*len(cell.center))
|
||||
dist = np.sqrt(r.dot(r))
|
||||
if dist < 0.2:
|
||||
return levels
|
||||
return levels - 1
|
||||
self.M.refine(function,balance=False)
|
||||
self.M.number(balance=False)
|
||||
# self.M.plotGrid(showIt=True)
|
||||
max_h = max([np.max(hi) for hi in self.M.h])
|
||||
return max_h
|
||||
|
||||
def getError(self):
|
||||
"""For given h, generate A[h], f and A(f) and return norm of error."""
|
||||
return 1.
|
||||
|
||||
@@ -124,13 +124,13 @@ if not _interpCython:
|
||||
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
|
||||
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
|
||||
|
||||
inds += [( ind_x1, ind_y2),
|
||||
( ind_x1, ind_y1),
|
||||
inds += [( ind_x1, ind_y1),
|
||||
( ind_x1, ind_y2),
|
||||
( ind_x2, ind_y1),
|
||||
( ind_x2, ind_y2)]
|
||||
|
||||
vals += [wx1*wy2,
|
||||
wx1*wy1,
|
||||
vals += [wx1*wy1,
|
||||
wx1*wy2,
|
||||
wx2*wy1,
|
||||
wx2*wy2]
|
||||
|
||||
@@ -152,8 +152,8 @@ if not _interpCython:
|
||||
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
|
||||
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
|
||||
|
||||
inds += [( ind_x1, ind_y2, ind_z1),
|
||||
( ind_x1, ind_y1, ind_z1),
|
||||
inds += [( ind_x1, ind_y1, ind_z1),
|
||||
( ind_x1, ind_y2, ind_z1),
|
||||
( ind_x2, ind_y1, ind_z1),
|
||||
( ind_x2, ind_y2, ind_z1),
|
||||
( ind_x1, ind_y1, ind_z2),
|
||||
@@ -161,8 +161,8 @@ if not _interpCython:
|
||||
( ind_x2, ind_y1, ind_z2),
|
||||
( ind_x2, ind_y2, ind_z2)]
|
||||
|
||||
vals += [wx1*wy2*wz1,
|
||||
wx1*wy1*wz1,
|
||||
vals += [wx1*wy1*wz1,
|
||||
wx1*wy2*wz1,
|
||||
wx2*wy1*wz1,
|
||||
wx2*wy2*wz1,
|
||||
wx1*wy1*wz2,
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -71,12 +71,12 @@ def _interpmat2D(np.ndarray[np.float64_t, ndim=2] locs,
|
||||
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
|
||||
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
|
||||
|
||||
inds += [( ind_x1, ind_y2),
|
||||
( ind_x1, ind_y1),
|
||||
inds += [( ind_x1, ind_y1),
|
||||
( ind_x1, ind_y2),
|
||||
( ind_x2, ind_y1),
|
||||
( ind_x2, ind_y2)]
|
||||
|
||||
vals += [wx1*wy2, wx1*wy1, wx2*wy1, wx2*wy2]
|
||||
vals += [wx1*wy1, wx1*wy2, wx2*wy1, wx2*wy2]
|
||||
|
||||
return inds, vals
|
||||
|
||||
@@ -98,8 +98,8 @@ def _interpmat3D(np.ndarray[np.float64_t, ndim=2] locs,
|
||||
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
|
||||
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
|
||||
|
||||
inds += [( ind_x1, ind_y2, ind_z1),
|
||||
( ind_x1, ind_y1, ind_z1),
|
||||
inds += [( ind_x1, ind_y1, ind_z1),
|
||||
( ind_x1, ind_y2, ind_z1),
|
||||
( ind_x2, ind_y1, ind_z1),
|
||||
( ind_x2, ind_y2, ind_z1),
|
||||
( ind_x1, ind_y1, ind_z2),
|
||||
@@ -107,8 +107,8 @@ def _interpmat3D(np.ndarray[np.float64_t, ndim=2] locs,
|
||||
( ind_x2, ind_y1, ind_z2),
|
||||
( ind_x2, ind_y2, ind_z2)]
|
||||
|
||||
vals += [wx1*wy2*wz1,
|
||||
wx1*wy1*wz1,
|
||||
vals += [wx1*wy1*wz1,
|
||||
wx1*wy2*wz1,
|
||||
wx2*wy1*wz1,
|
||||
wx2*wy2*wz1,
|
||||
wx1*wy1*wz2,
|
||||
|
||||
+25
-25
@@ -149,7 +149,7 @@ def readUBCTensorModel(fileName, mesh):
|
||||
|
||||
Input:
|
||||
:param fileName, path to the UBC GIF mesh file to read
|
||||
:param mesh, TensorMesh object, mesh that coresponds to the model
|
||||
:param mesh, TensorMesh object, mesh that coresponds to the model
|
||||
|
||||
Output:
|
||||
:return numpy array, model with TensorMesh ordered
|
||||
@@ -170,7 +170,7 @@ def writeUBCTensorMesh(fileName, mesh):
|
||||
|
||||
:param str fileName: File to write to
|
||||
:param simpeg.Mesh.TensorMesh mesh: The mesh
|
||||
|
||||
|
||||
"""
|
||||
assert mesh.dim == 3
|
||||
s = ''
|
||||
@@ -216,7 +216,7 @@ def readVTRFile(fileName):
|
||||
Output:
|
||||
:return SimPEG TensorMesh object
|
||||
:return SimPEG model dictionary
|
||||
|
||||
|
||||
"""
|
||||
# Import
|
||||
from vtk import vtkXMLRectilinearGridReader as vtrFileReader
|
||||
@@ -324,56 +324,56 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
|
||||
Extracts Core Mesh from Global mesh
|
||||
xyzlim: 2D array [ndim x 2]
|
||||
mesh: SimPEG mesh
|
||||
This function ouputs:
|
||||
This function ouputs:
|
||||
- actind: corresponding boolean index from global to core
|
||||
- meshcore: core SimPEG mesh
|
||||
- meshcore: core SimPEG mesh
|
||||
Warning: 1D and 2D has not been tested
|
||||
"""
|
||||
from SimPEG import Mesh
|
||||
if mesh.dim ==1:
|
||||
xyzlim = xyzlim.flatten()
|
||||
xmin, xmax = xyzlim[0], xyzlim[1]
|
||||
|
||||
xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
|
||||
|
||||
|
||||
xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
|
||||
|
||||
xc = mesh.vectorCCx[xind]
|
||||
|
||||
hx = mesh.hx[xind]
|
||||
|
||||
|
||||
x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
|
||||
|
||||
|
||||
meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
|
||||
|
||||
|
||||
actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
|
||||
|
||||
|
||||
elif mesh.dim ==2:
|
||||
xmin, xmax = xyzlim[0,0], xyzlim[0,1]
|
||||
ymin, ymax = xyzlim[1,0], xyzlim[1,1]
|
||||
|
||||
yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
|
||||
zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
|
||||
zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
|
||||
|
||||
xc = mesh.vectorCCx[xind]
|
||||
yc = mesh.vectorCCy[yind]
|
||||
|
||||
hx = mesh.hx[xind]
|
||||
hy = mesh.hy[yind]
|
||||
|
||||
|
||||
x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
|
||||
|
||||
|
||||
meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
|
||||
|
||||
|
||||
actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
|
||||
& (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
|
||||
|
||||
|
||||
elif mesh.dim==3:
|
||||
xmin, xmax = xyzlim[0,0], xyzlim[0,1]
|
||||
ymin, ymax = xyzlim[1,0], xyzlim[1,1]
|
||||
zmin, zmax = xyzlim[2,0], xyzlim[2,1]
|
||||
|
||||
|
||||
xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
|
||||
yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
|
||||
zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
|
||||
zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
|
||||
|
||||
xc = mesh.vectorCCx[xind]
|
||||
yc = mesh.vectorCCy[yind]
|
||||
@@ -382,19 +382,19 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
|
||||
hx = mesh.hx[xind]
|
||||
hy = mesh.hy[yind]
|
||||
hz = mesh.hz[zind]
|
||||
|
||||
|
||||
x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
|
||||
|
||||
|
||||
meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
|
||||
|
||||
|
||||
actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
|
||||
& (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
|
||||
& (mesh.gridCC[:,2]>zmin) & (mesh.gridCC[:,2]<zmax)
|
||||
|
||||
|
||||
else:
|
||||
raise(Exception("Not implemented!"))
|
||||
|
||||
|
||||
|
||||
|
||||
return actind, meshCore
|
||||
|
||||
|
||||
|
||||
@@ -5,10 +5,15 @@ SimPEG is a python package for simulation and gradient based
|
||||
parameter estimation in the context of geophysical applications.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
import os
|
||||
import sys
|
||||
import subprocess
|
||||
|
||||
from distutils.core import setup
|
||||
from setuptools import find_packages
|
||||
from Cython.Build import cythonize
|
||||
import numpy as np
|
||||
from distutils.extension import Extension
|
||||
|
||||
CLASSIFIERS = [
|
||||
'Development Status :: 4 - Beta',
|
||||
@@ -26,6 +31,44 @@ CLASSIFIERS = [
|
||||
'Natural Language :: English',
|
||||
]
|
||||
|
||||
args = sys.argv[1:]
|
||||
|
||||
# Make a `cleanall` rule to get rid of intermediate and library files
|
||||
if "cleanall" in args:
|
||||
print "Deleting cython files..."
|
||||
# Just in case the build directory was created by accident,
|
||||
# note that shell=True should be OK here because the command is constant.
|
||||
subprocess.Popen("rm -rf build", shell=True, executable="/bin/bash")
|
||||
subprocess.Popen("find . -name \*.c -type f -delete", shell=True, executable="/bin/bash")
|
||||
subprocess.Popen("find . -name \*.so -type f -delete", shell=True, executable="/bin/bash")
|
||||
# Now do a normal clean
|
||||
sys.argv[sys.argv.index('cleanall')] = "clean"
|
||||
|
||||
# We want to always use build_ext --inplace
|
||||
if args.count("build_ext") > 0 and args.count("--inplace") == 0:
|
||||
sys.argv.insert(sys.argv.index("build_ext")+1, "--inplace")
|
||||
|
||||
try:
|
||||
from Cython.Build import cythonize
|
||||
from Cython.Distutils import build_ext
|
||||
cythonKwargs = dict(cmdclass={'build_ext': build_ext})
|
||||
USE_CYTHON = True
|
||||
except Exception, e:
|
||||
USE_CYTHON = False
|
||||
cythonKwargs = dict()
|
||||
|
||||
ext = '.pyx' if USE_CYTHON else '.c'
|
||||
|
||||
cython_files = [
|
||||
"SimPEG/Utils/interputils_cython",
|
||||
"SimPEG/Mesh/TreeUtils"
|
||||
]
|
||||
extensions = [Extension(f, [f+ext]) for f in cython_files]
|
||||
|
||||
if USE_CYTHON and "cleanall" not in args:
|
||||
from Cython.Build import cythonize
|
||||
extensions = cythonize(extensions)
|
||||
|
||||
import os, os.path
|
||||
|
||||
with open("README.rst") as f:
|
||||
@@ -51,5 +94,6 @@ setup(
|
||||
platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],
|
||||
use_2to3 = False,
|
||||
include_dirs=[np.get_include()],
|
||||
ext_modules = cythonize('SimPEG/Utils/interputils_cython.pyx')
|
||||
ext_modules = extensions,
|
||||
**cythonKwargs
|
||||
)
|
||||
|
||||
@@ -0,0 +1,182 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG import Utils, Tests
|
||||
|
||||
MESHTYPES = ['uniformTree'] #['randomTree', 'uniformTree']
|
||||
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)))
|
||||
|
||||
|
||||
plotIt = False
|
||||
|
||||
|
||||
MESHTYPES = ['uniformTree','notatreeTree']
|
||||
|
||||
|
||||
"""
|
||||
|
||||
Face interpolation is O(h)
|
||||
Edge interpolation is O(h^2)
|
||||
|
||||
"""
|
||||
|
||||
class TestInterpolation2d(Tests.OrderTest):
|
||||
name = "Interpolation 2D"
|
||||
np.random.seed(1)
|
||||
LOCS = np.random.rand(50,2)*0.6+0.2
|
||||
# LOCS = np.c_[np.ones(100)*0.51, np.linspace(0.3,0.7,100)]
|
||||
meshTypes = MESHTYPES
|
||||
# tolerance = TOLERANCES
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32]
|
||||
expectedOrders = 1
|
||||
|
||||
def getError(self):
|
||||
funX = lambda x, y: np.cos(2.*np.pi*y)*np.cos(2.*np.pi*x) + x
|
||||
funY = lambda x, y: np.cos(2.*np.pi*x)*np.cos(2.*np.pi*y) + y
|
||||
|
||||
# self.LOCS = self.M.gridCC
|
||||
|
||||
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)
|
||||
|
||||
P = self.M.getInterpolationMat(self.LOCS, self.type)
|
||||
# print P
|
||||
comp = P*grid
|
||||
|
||||
err = np.linalg.norm((comp - ana), np.inf)
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
ax = plt.subplot(211)
|
||||
self.M.plotGrid(ax=ax)
|
||||
plt.plot(self.LOCS[:,0],self.LOCS[:,1], 'mx')
|
||||
# ax = plt.subplot(111)
|
||||
# self.M.plotImage(call2(funX, self.M.gridCC),ax=ax)
|
||||
ax = plt.subplot(212)
|
||||
plt.plot(self.LOCS[:,1],comp, 'bx')
|
||||
plt.plot(self.LOCS[:,1],ana, 'ro')
|
||||
plt.show()
|
||||
return err
|
||||
|
||||
def test_orderFx(self):
|
||||
self.type = 'Fx'
|
||||
self.name = 'TreeMesh Interpolation 2D: Fx'
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFy(self):
|
||||
self.type = 'Fy'
|
||||
self.name = 'TreeMesh Interpolation 2D: Fy'
|
||||
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]
|
||||
|
||||
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.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
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.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
def test_orderFy(self):
|
||||
self.type = 'Fy'
|
||||
self.name = 'Interpolation 3D: Fy'
|
||||
self.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
def test_orderFz(self):
|
||||
self.type = 'Fz'
|
||||
self.name = 'Interpolation 3D: Fz'
|
||||
self.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
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()
|
||||
+280
-479
@@ -1,504 +1,305 @@
|
||||
from SimPEG.Mesh import TensorMesh
|
||||
from SimPEG.Mesh.TreeMesh import TreeMesh, TreeFace, TreeCell
|
||||
from SimPEG import Mesh, Tests
|
||||
from SimPEG.Mesh.TreeMesh import CellLookUpException
|
||||
import numpy as np
|
||||
import unittest
|
||||
import matplotlib.pyplot as plt
|
||||
import unittest
|
||||
|
||||
TOL = 1e-10
|
||||
TOL = 1e-8
|
||||
|
||||
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
|
||||
class TestSimpleQuadTree(unittest.TestCase):
|
||||
|
||||
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
|
||||
nc = 8
|
||||
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h = [hi/np.sum(hi) for hi in [h1, h2]] # normalize
|
||||
M = Mesh.TreeMesh(h)
|
||||
M._refineCell([0,0,0])
|
||||
M._refineCell([0,0,1])
|
||||
M.number()
|
||||
# M.plotGrid(showIt=True)
|
||||
print M
|
||||
assert M.nhFx == 2
|
||||
assert M.nFx == 9
|
||||
|
||||
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)
|
||||
assert np.allclose(M.vol.sum(), 1.0)
|
||||
|
||||
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)
|
||||
assert np.allclose(np.r_[M._areaFxFull, M._areaFyFull], M._deflationMatrix('F') * M.area)
|
||||
|
||||
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_refine(self):
|
||||
M = Mesh.TreeMesh([4,4,4])
|
||||
M.refine(1)
|
||||
assert M.nC == 8
|
||||
M.refine(0)
|
||||
assert M.nC == 8
|
||||
M.corsen(0)
|
||||
assert M.nC == 1
|
||||
|
||||
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_corsen(self):
|
||||
nc = 8
|
||||
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h = [hi/np.sum(hi) for hi in [h1, h2]] # normalize
|
||||
M = Mesh.TreeMesh(h)
|
||||
M._refineCell([0,0,0])
|
||||
M._refineCell([0,0,1])
|
||||
self.assertRaises(CellLookUpException, M._refineCell, [0,0,1])
|
||||
assert M._index([0,0,1]) not in M
|
||||
assert M._index([0,0,2]) in M
|
||||
assert M._index([2,0,2]) in M
|
||||
assert M._index([0,2,2]) in M
|
||||
assert M._index([2,2,2]) in M
|
||||
|
||||
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)
|
||||
self.assertRaises(CellLookUpException, M._corsenCell, [0,0,1])
|
||||
M._corsenCell([0,0,2])
|
||||
assert M._index([0,0,1]) in M
|
||||
assert M._index([0,0,2]) not in M
|
||||
assert M._index([2,0,2]) not in M
|
||||
assert M._index([0,2,2]) not in M
|
||||
assert M._index([2,2,2]) not in M
|
||||
M._refineCell([0,0,1])
|
||||
|
||||
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)
|
||||
self.assertRaises(CellLookUpException, M._corsenCell, [0,0,1])
|
||||
M._corsenCell([2,0,2])
|
||||
assert M._index([0,0,1]) in M
|
||||
assert M._index([0,0,2]) not in M
|
||||
assert M._index([2,0,2]) not in M
|
||||
assert M._index([0,2,2]) not in M
|
||||
assert M._index([2,2,2]) not in M
|
||||
M._refineCell([0,0,1])
|
||||
|
||||
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)
|
||||
self.assertRaises(CellLookUpException, M._corsenCell, [0,0,1])
|
||||
M._corsenCell([0,2,2])
|
||||
assert M._index([0,0,1]) in M
|
||||
assert M._index([0,0,2]) not in M
|
||||
assert M._index([2,0,2]) not in M
|
||||
assert M._index([0,2,2]) not in M
|
||||
assert M._index([2,2,2]) not in M
|
||||
M._refineCell([0,0,1])
|
||||
|
||||
|
||||
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])
|
||||
self.assertRaises(CellLookUpException, M._corsenCell, [0,0,1])
|
||||
M._corsenCell([2,2,2])
|
||||
assert M._index([0,0,1]) in M
|
||||
assert M._index([0,0,2]) not in M
|
||||
assert M._index([2,0,2]) not in M
|
||||
assert M._index([0,2,2]) not in M
|
||||
assert M._index([2,2,2]) not in M
|
||||
|
||||
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)
|
||||
hx, hy = np.r_[1.,2,3,4], np.r_[5.,6,7,8]
|
||||
T = Mesh.TreeMesh([hx, hy], levels=2)
|
||||
T.refine(lambda xc:2)
|
||||
# T.plotGrid(showIt=True)
|
||||
M = Mesh.TensorMesh([hx, hy])
|
||||
assert M.nC == T.nC
|
||||
assert M.nF == T.nF
|
||||
assert M.nFx == T.nFx
|
||||
assert M.nFy == T.nFy
|
||||
assert M.nE == T.nE
|
||||
assert M.nEx == T.nEx
|
||||
assert M.nEy == T.nEy
|
||||
assert np.allclose(M.area, T.permuteF*T.area)
|
||||
assert np.allclose(M.edge, T.permuteE*T.edge)
|
||||
assert np.allclose(M.vol, T.permuteCC*T.vol)
|
||||
|
||||
# plt.subplot(211).spy(M.faceDiv)
|
||||
# plt.subplot(212).spy(T.permuteCC*T.faceDiv*T.permuteF.T)
|
||||
# plt.show()
|
||||
|
||||
assert (M.faceDiv - T.permuteCC*T.faceDiv*T.permuteF.T).nnz == 0
|
||||
|
||||
|
||||
class TestOcTree(unittest.TestCase):
|
||||
|
||||
def test_counts(self):
|
||||
nc = 8
|
||||
h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
|
||||
h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize
|
||||
M = Mesh.TreeMesh(h, levels=3)
|
||||
M._refineCell([0,0,0,0])
|
||||
M._refineCell([0,0,0,1])
|
||||
M.number()
|
||||
# M.plotGrid(showIt=True)
|
||||
# assert M.nhFx == 2
|
||||
# assert M.nFx == 9
|
||||
|
||||
assert np.allclose(M.vol.sum(), 1.0)
|
||||
|
||||
# assert np.allclose(M._areaFxFull, (M._deflationMatrix('F') * M.area)[:M.ntFx])
|
||||
# assert np.allclose(M._areaFyFull, (M._deflationMatrix('F') * M.area)[M.ntFx:(M.ntFx+M.ntFy)])
|
||||
# assert np.allclose(M._areaFzFull, (M._deflationMatrix('F') * M.area)[(M.ntFx+M.ntFy):])
|
||||
|
||||
# assert np.allclose(M._edgeExFull, (M._deflationMatrix('E') * M.edge)[:M.ntEx])
|
||||
# assert np.allclose(M._edgeEyFull, (M._deflationMatrix('E') * M.edge)[M.ntEx:(M.ntEx+M.ntEy)])
|
||||
# assert np.allclose(M._edgeEzFull, (M._deflationMatrix('E') * M.edge)[(M.ntEx+M.ntEy):])
|
||||
|
||||
def test_faceDiv(self):
|
||||
|
||||
hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
|
||||
M = Mesh.TreeMesh([hx, hy, hz], levels=2)
|
||||
M.refine(lambda xc:2)
|
||||
# M.plotGrid(showIt=True)
|
||||
Mr = Mesh.TensorMesh([hx, hy, hz])
|
||||
assert M.nC == Mr.nC
|
||||
assert M.nF == Mr.nF
|
||||
assert M.nFx == Mr.nFx
|
||||
assert M.nFy == Mr.nFy
|
||||
assert M.nE == Mr.nE
|
||||
assert M.nEx == Mr.nEx
|
||||
assert M.nEy == Mr.nEy
|
||||
assert np.allclose(Mr.area, M.permuteF*M.area)
|
||||
assert np.allclose(Mr.edge, M.permuteE*M.edge)
|
||||
assert np.allclose(Mr.vol, M.permuteCC*M.vol)
|
||||
|
||||
# plt.subplot(211).spy(Mr.faceDiv)
|
||||
# plt.subplot(212).spy(M.permuteCC*M.faceDiv*M.permuteF.T)
|
||||
# plt.show()
|
||||
|
||||
assert (Mr.faceDiv - M.permuteCC*M.faceDiv*M.permuteF.T).nnz == 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)
|
||||
hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
|
||||
M = Mesh.TreeMesh([hx, hy, hz], levels=2)
|
||||
M.refine(lambda xc:2)
|
||||
# M.plotGrid(showIt=True)
|
||||
Mr = Mesh.TensorMesh([hx, hy, hz])
|
||||
|
||||
# plt.subplot(211).spy(Mr.faceDiv)
|
||||
# plt.subplot(212).spy(M.permuteCC.T*M.faceDiv*M.permuteF)
|
||||
# plt.show()
|
||||
|
||||
assert (Mr.edgeCurl - M.permuteF*M.edgeCurl*M.permuteE.T).nnz == 0
|
||||
|
||||
def test_faceInnerProduct(self):
|
||||
|
||||
hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
|
||||
# hx, hy, hz = [[(1,4)], [(1,4)], [(1,4)]]
|
||||
|
||||
M = Mesh.TreeMesh([hx, hy, hz], levels=2)
|
||||
M.refine(lambda xc:2)
|
||||
# M.plotGrid(showIt=True)
|
||||
Mr = Mesh.TensorMesh([hx, hy, hz])
|
||||
|
||||
# plt.subplot(211).spy(Mr.getFaceInnerProduct())
|
||||
# plt.subplot(212).spy(M.getFaceInnerProduct())
|
||||
# plt.show()
|
||||
|
||||
# print M.nC, M.nF, M.getFaceInnerProduct().shape, M.permuteF.shape
|
||||
|
||||
assert np.allclose(Mr.getFaceInnerProduct().todense(), (M.permuteF * M.getFaceInnerProduct() * M.permuteF.T).todense())
|
||||
assert np.allclose(Mr.getEdgeInnerProduct().todense(), (M.permuteE * M.getEdgeInnerProduct() * M.permuteE.T).todense())
|
||||
|
||||
def test_VectorIdenties(self):
|
||||
hx, hy, hz = [[(1,4)], [(1,4)], [(1,4)]]
|
||||
|
||||
M = Mesh.TreeMesh([hx, hy, hz], levels=2)
|
||||
Mr = Mesh.TensorMesh([hx, hy, hz])
|
||||
|
||||
assert (M.faceDiv * M.edgeCurl).nnz == 0
|
||||
assert (Mr.faceDiv * Mr.edgeCurl).nnz == 0
|
||||
|
||||
hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
|
||||
|
||||
M = Mesh.TreeMesh([hx, hy, hz], levels=2)
|
||||
Mr = Mesh.TensorMesh([hx, hy, hz])
|
||||
|
||||
assert np.max(np.abs((M.faceDiv * M.edgeCurl).todense().flatten())) < TOL
|
||||
assert np.max(np.abs((Mr.faceDiv * Mr.edgeCurl).todense().flatten())) < TOL
|
||||
|
||||
class Test2DInterpolation(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
def topo(x):
|
||||
return np.sin(x*(2.*np.pi))*0.3 + 0.5
|
||||
|
||||
def function(cell):
|
||||
r = cell.center - np.array([0.5]*len(cell.center))
|
||||
dist1 = np.sqrt(r.dot(r)) - 0.08
|
||||
dist2 = np.abs(cell.center[-1] - topo(cell.center[0]))
|
||||
|
||||
dist = min([dist1,dist2])
|
||||
# if dist < 0.05:
|
||||
# return 5
|
||||
if dist < 0.05:
|
||||
return 6
|
||||
if dist < 0.2:
|
||||
return 5
|
||||
if dist < 0.3:
|
||||
return 4
|
||||
if dist < 1.0:
|
||||
return 3
|
||||
else:
|
||||
return 0
|
||||
|
||||
M = Mesh.TreeMesh([64,64],levels=6)
|
||||
M.refine(function)
|
||||
self.M = M
|
||||
|
||||
def test_fx(self):
|
||||
r = np.random.rand(self.M.nFx)
|
||||
P = self.M.getInterpolationMat(self.M.gridFx, 'Fx')
|
||||
assert np.abs(P[:,:self.M.nFx]*r - r).max() < TOL
|
||||
|
||||
def test_fy(self):
|
||||
r = np.random.rand(self.M.nFy)
|
||||
P = self.M.getInterpolationMat(self.M.gridFy, 'Fy')
|
||||
assert np.abs(P[:,self.M.nFx:]*r - r).max() < TOL
|
||||
|
||||
|
||||
class Test3DInterpolation(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
def function(cell):
|
||||
r = cell.center - np.array([0.5]*len(cell.center))
|
||||
dist = np.sqrt(r.dot(r))
|
||||
if dist < 0.2:
|
||||
return 4
|
||||
if dist < 0.3:
|
||||
return 3
|
||||
if dist < 1.0:
|
||||
return 2
|
||||
else:
|
||||
return 0
|
||||
|
||||
M = Mesh.TreeMesh([16,16,16],levels=4)
|
||||
M.refine(function)
|
||||
# M.plotGrid(showIt=True)
|
||||
self.M = M
|
||||
|
||||
def test_Fx(self):
|
||||
r = np.random.rand(self.M.nFx)
|
||||
P = self.M.getInterpolationMat(self.M.gridFx, 'Fx')
|
||||
assert np.abs(P[:,:self.M.nFx]*r - r).max() < TOL
|
||||
|
||||
def test_Fy(self):
|
||||
r = np.random.rand(self.M.nFy)
|
||||
P = self.M.getInterpolationMat(self.M.gridFy, 'Fy')
|
||||
assert np.abs(P[:,self.M.nFx:(self.M.nFx+self.M.nFy)]*r - r).max() < TOL
|
||||
|
||||
def test_Fz(self):
|
||||
r = np.random.rand(self.M.nFz)
|
||||
P = self.M.getInterpolationMat(self.M.gridFz, 'Fz')
|
||||
assert np.abs(P[:,(self.M.nFx+self.M.nFy):]*r - r).max() < TOL
|
||||
|
||||
def test_Ex(self):
|
||||
r = np.random.rand(self.M.nEx)
|
||||
P = self.M.getInterpolationMat(self.M.gridEx, 'Ex')
|
||||
assert np.abs(P[:,:self.M.nEx]*r - r).max() < TOL
|
||||
|
||||
def test_Ey(self):
|
||||
r = np.random.rand(self.M.nEy)
|
||||
P = self.M.getInterpolationMat(self.M.gridEy, 'Ey')
|
||||
assert np.abs(P[:,self.M.nEx:(self.M.nEx+self.M.nEy)]*r - r).max() < TOL
|
||||
|
||||
def test_Ez(self):
|
||||
r = np.random.rand(self.M.nEz)
|
||||
P = self.M.getInterpolationMat(self.M.gridEz, 'Ez')
|
||||
assert np.abs(P[:,(self.M.nEx+self.M.nEy):]*r - r).max() < TOL
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -0,0 +1,675 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG import Utils, Tests
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
MESHTYPES = ['uniformTree'] #['randomTree', 'uniformTree']
|
||||
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)))
|
||||
|
||||
|
||||
plotIt = False
|
||||
|
||||
class TestFaceDiv2D(Tests.OrderTest):
|
||||
name = "Face Divergence 2D"
|
||||
meshTypes = MESHTYPES
|
||||
meshDimension = 2
|
||||
meshSizes = [16, 32]
|
||||
|
||||
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)
|
||||
|
||||
# self.M.plotImage(divF-divF_ana, showIt=True)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestFaceDiv3D(Tests.OrderTest):
|
||||
name = "Face Divergence 3D"
|
||||
meshTypes = MESHTYPES
|
||||
meshSizes = [8, 16]
|
||||
|
||||
def getError(self):
|
||||
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)
|
||||
|
||||
return np.linalg.norm((divF-divF_ana), np.inf)
|
||||
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestCurl(Tests.OrderTest):
|
||||
name = "Curl"
|
||||
meshTypes = ['notatreeTree', 'uniformTree'] #, 'randomTree']#, 'uniformTree']
|
||||
meshSizes = [8, 16]#, 32]
|
||||
expectedOrders = [2,1] # This is due to linear interpolation in the Re projection
|
||||
|
||||
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)
|
||||
|
||||
err = np.linalg.norm((curlE - curlE_ana), np.inf)
|
||||
# err = np.linalg.norm((curlE - curlE_ana)*self.M.area, 2)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestNodalGrad(Tests.OrderTest):
|
||||
name = "Nodal Gradient"
|
||||
meshTypes = ['notatreeTree', 'uniformTree'] #['randomTree', 'uniformTree']
|
||||
meshSizes = [8, 16]#, 32]
|
||||
expectedOrders = [2,1]
|
||||
|
||||
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(Tests.OrderTest):
|
||||
name = "Nodal Gradient 2D"
|
||||
meshTypes = ['notatreeTree', 'uniformTree'] #['randomTree', 'uniformTree']
|
||||
meshSizes = [8, 16]#, 32]
|
||||
expectedOrders = [2,1]
|
||||
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 TestTreeInnerProducts(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTree', 'notatreeTree'] #['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
|
||||
meshDimension = 3
|
||||
meshSizes = [4, 8]
|
||||
|
||||
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 TestTreeInnerProducts2D(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTree']
|
||||
meshDimension = 2
|
||||
meshSizes = [4, 8]
|
||||
|
||||
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 TestTreeAveraging2D(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['notatreeTree', 'uniformTree']#, 'randomTree']
|
||||
meshDimension = 2
|
||||
meshSizes = [4,8,16]
|
||||
expectedOrders = [2,1]
|
||||
|
||||
def getError(self):
|
||||
if plotIt:
|
||||
plt.spy(self.getAve(self.M))
|
||||
plt.show()
|
||||
|
||||
num = self.getAve(self.M) * self.getHere(self.M)
|
||||
err = np.linalg.norm((self.getThere(self.M)-num), np.inf)
|
||||
|
||||
if plotIt:
|
||||
self.M.plotImage(self.getThere(self.M)-num)
|
||||
plt.show()
|
||||
plt.tight_layout
|
||||
|
||||
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, np.r_[M.gridFx, M.gridFy])]
|
||||
self.getThere = lambda M: call2(fun, M.gridCC)
|
||||
self.getAve = lambda M: M.aveF2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFx2CC(self):
|
||||
self.name = "Averaging 2D: Fx2CC"
|
||||
funX = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: np.r_[call2(funX, M.gridFx)]
|
||||
self.getThere = lambda M: np.r_[call2(funX, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveFx2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFy2CC(self):
|
||||
self.name = "Averaging 2D: Fy2CC"
|
||||
funY = lambda x, y: (np.cos(y)*np.sin(x))
|
||||
self.getHere = lambda M: np.r_[call2(funY, M.gridFy)]
|
||||
self.getThere = lambda M: np.r_[call2(funY, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveFy2CC
|
||||
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
|
||||
|
||||
class TestAveraging3D(Tests.OrderTest):
|
||||
name = "Averaging 3D"
|
||||
meshTypes = ['notatreeTree', 'uniformTree']#, 'randomTree']
|
||||
meshDimension = 3
|
||||
meshSizes = [8,16]
|
||||
expectedOrders = [2,1]
|
||||
|
||||
def getError(self):
|
||||
if plotIt:
|
||||
plt.spy(self.getAve(self.M))
|
||||
plt.show()
|
||||
|
||||
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_orderFx2CC(self):
|
||||
self.name = "Averaging 3D: Fx2CC"
|
||||
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funX, M.gridFx)]
|
||||
self.getThere = lambda M: np.r_[call3(funX, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveFx2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFy2CC(self):
|
||||
self.name = "Averaging 3D: Fy2CC"
|
||||
funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funY, M.gridFy)]
|
||||
self.getThere = lambda M: np.r_[call3(funY, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveFy2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderFz2CC(self):
|
||||
self.name = "Averaging 3D: Fz2CC"
|
||||
funZ = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funZ, M.gridFz)]
|
||||
self.getThere = lambda M: np.r_[call3(funZ, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveFz2CC
|
||||
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_orderEx2CC(self):
|
||||
self.name = "Averaging 3D: Ex2CC"
|
||||
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funX, M.gridEx)]
|
||||
self.getThere = lambda M: np.r_[call3(funX, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveEx2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEy2CC(self):
|
||||
self.name = "Averaging 3D: Ey2CC"
|
||||
funY = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funY, M.gridEy)]
|
||||
self.getThere = lambda M: np.r_[call3(funY, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveEy2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderEz2CC(self):
|
||||
self.name = "Averaging 3D: Ez2CC"
|
||||
funZ = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
self.getHere = lambda M: np.r_[call3(funZ, M.gridEz)]
|
||||
self.getThere = lambda M: np.r_[call3(funZ, M.gridCC)]
|
||||
self.getAve = lambda M: M.aveEz2CC
|
||||
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()
|
||||
@@ -10,7 +10,9 @@ class TestInnerProductsDerivs(unittest.TestCase):
|
||||
hRect = Utils.exampleLrmGrid(h,'rotate')
|
||||
mesh = Mesh.CurvilinearMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
mesh = Mesh.TreeMesh(h)
|
||||
mesh = Mesh.TreeMesh(h, levels=3)
|
||||
mesh.refine(lambda xc: 3)
|
||||
mesh.number(balance=False)
|
||||
elif meshType == 'Tensor':
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nF)
|
||||
@@ -27,7 +29,9 @@ class TestInnerProductsDerivs(unittest.TestCase):
|
||||
hRect = Utils.exampleLrmGrid(h,'rotate')
|
||||
mesh = Mesh.CurvilinearMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
mesh = Mesh.TreeMesh(h)
|
||||
mesh = Mesh.TreeMesh(h, levels=3)
|
||||
mesh.refine(lambda xc: 3)
|
||||
mesh.number(balance=False)
|
||||
elif meshType == 'Tensor':
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nE)
|
||||
@@ -197,67 +201,65 @@ class TestInnerProductsDerivs(unittest.TestCase):
|
||||
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'))
|
||||
self.assertTrue(self.doTestFace([8, 8],0, False, 'Tree'))
|
||||
def test_FaceIP_3D_float_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],0, False, 'Tree'))
|
||||
def test_FaceIP_2D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8],1, False, 'Tree'))
|
||||
def test_FaceIP_3D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],1, False, 'Tree'))
|
||||
def test_FaceIP_2D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8],2, False, 'Tree'))
|
||||
def test_FaceIP_3D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],3, False, 'Tree'))
|
||||
def test_FaceIP_2D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8],3, False, 'Tree'))
|
||||
def test_FaceIP_3D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],6, False, 'Tree'))
|
||||
|
||||
def test_FaceIP_2D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8],0, True, 'Tree'))
|
||||
def test_FaceIP_3D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],0, True, 'Tree'))
|
||||
def test_FaceIP_2D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8],1, True, 'Tree'))
|
||||
def test_FaceIP_3D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],1, True, 'Tree'))
|
||||
def test_FaceIP_2D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8],2, True, 'Tree'))
|
||||
def test_FaceIP_3D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tree'))
|
||||
self.assertTrue(self.doTestFace([8, 8, 8],3, True, 'Tree'))
|
||||
|
||||
def test_EdgeIP_2D_float_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tree'))
|
||||
# def test_EdgeIP_2D_float_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],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'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],0, False, 'Tree'))
|
||||
# def test_EdgeIP_2D_isotropic_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],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'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],1, False, 'Tree'))
|
||||
# def test_EdgeIP_2D_anisotropic_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],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'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],3, False, 'Tree'))
|
||||
# def test_EdgeIP_2D_tensor_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],3, False, 'Tree'))
|
||||
def test_EdgeIP_3D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tree'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],6, False, 'Tree'))
|
||||
|
||||
def test_EdgeIP_2D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tree'))
|
||||
# def test_EdgeIP_2D_float_fast_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],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'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],0, True, 'Tree'))
|
||||
# def test_EdgeIP_2D_isotropic_fast_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],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'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],1, True, 'Tree'))
|
||||
# def test_EdgeIP_2D_anisotropic_fast_Tree(self):
|
||||
# self.assertTrue(self.doTestEdge([8, 8],2, True, 'Tree'))
|
||||
def test_EdgeIP_3D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tree'))
|
||||
self.assertTrue(self.doTestEdge([8, 8, 8],3, True, 'Tree'))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user