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

started testing of faceDiv

Browse Source
This commit is contained in:
rowanc1
2014-02-16 18:05:11 -08:00
parent d5c09ce036
commit 8b5ca46ecf
3 changed files with 117 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/CylMesh.py
+46
View File
@@ -136,6 +136,45 @@ class CylMesh(BaseTensorMesh):
"""Nodal grid vector (1D) in the y direction."""
return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5
def getTensor(self, locType):
""" Returns a tensor list.
:param str locType: What tensor (see below)
:rtype: list
:return: list of the tensors that make up the mesh.
locType can be::
'Ex' -> x-component of field defined on edges
'Ey' -> y-component of field defined on edges
'Ez' -> z-component of field defined on edges
'Fx' -> x-component of field defined on faces
'Fy' -> y-component of field defined on faces
'Fz' -> z-component of field defined on faces
'N' -> scalar field defined on nodes
'CC' -> scalar field defined on cell centers
"""
if locType is 'Fx':
ten = [self.vectorNx , self.vectorCCy, self.vectorCCz]
elif locType is 'Fy':
ten = [self.vectorCCx, self.vectorNy , self.vectorCCz]
elif locType is 'Fz':
ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ]
elif locType is 'Ex':
ten = [self.vectorCCx, self.vectorNy , self.vectorNz ]
elif locType is 'Ey':
ten = [self.vectorNx , self.vectorCCy, self.vectorNz ]
elif locType is 'Ez':
ten = [self.vectorNx , self.vectorNy , self.vectorCCz]
elif locType is 'CC':
ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz]
elif locType is 'N':
ten = [self.vectorNx , self.vectorNy , self.vectorNz ]
return [t for t in ten if t is not None]
@property
def edge(self):
"""Edge lengths"""
@@ -219,6 +258,13 @@ class CylMesh(BaseTensorMesh):
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
return self._faceDivz
@property
def cellGrad(self):
"""The cell centered Gradient, takes you to cell faces."""
raise NotImplementedError('Cell Grad is not yet implemented.')
@property
def nodalGrad(self):
"""Construct gradient operator (nodes to edges)."""
SimPEG/Tests/TestUtils.py
+16 -5
View File
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
from numpy.linalg import norm
from SimPEG.Utils import mkvc, sdiag
from SimPEG import Utils
from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh
from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh, CylMesh
import numpy as np
import scipy.sparse as sp
import unittest
@@ -88,10 +88,7 @@ class OrderTest(unittest.TestCase):
"""
if 'TensorMesh' in self._meshType:
if 'uniform' in self._meshType:
h1 = np.ones(nc)/nc
h2 = np.ones(nc)/nc
h3 = np.ones(nc)/nc
h = [h1, h2, h3]
h = [nc, nc, nc]
elif 'random' in self._meshType:
h1 = np.random.rand(nc)
h2 = np.random.rand(nc)
@@ -104,6 +101,20 @@ class OrderTest(unittest.TestCase):
max_h = max([np.max(hi) for hi in self.M.h])
return max_h
elif 'CylMesh' in self._meshType:
if 'uniform' in self._meshType:
h = [nc, nc, nc]
else:
raise Exception('Unexpected meshType')
if self.meshDimension == 2:
self.M = CylMesh([h[0], 1, h[2]])
max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]])
elif self.meshDimension == 3:
self.M = CylMesh(h)
max_h = max([np.max(hi) for hi in self.M.h])
return max_h
elif 'LOM' in self._meshType:
if 'uniform' in self._meshType:
kwrd = 'rect'
SimPEG/Tests/test_cylMesh.py
+55 -2
View File
@@ -1,6 +1,7 @@
import unittest
import sys
from SimPEG import *
from TestUtils import OrderTest
class TestCyl2DMesh(unittest.TestCase):
@@ -75,8 +76,60 @@ class TestCyl2DMesh(unittest.TestCase):
vol = np.r_[2*a,a]
self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0)
def test_faceDiv(self):
print self.mesh.faceDiv
def test_gridSizes(self):
self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3))
# self.assertTrue(self.mesh.gridN.shape == (self.mesh.nN, 3))
self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3))
# self.assertTrue(self.mesh.gridFy.shape == (self.mesh.nFy, 3))
self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3))
# self.assertTrue(self.mesh.gridEx.shape == (self.mesh.nEx, 3))
self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3))
# self.assertTrue(self.mesh.gridEz.shape == (self.mesh.nEz, 3))
MESHTYPES = ['uniformCylMesh']
MESHDIMENSION = 2
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
class TestFaceDiv(OrderTest):
name = "FaceDiv"
meshTypes = MESHTYPES
meshDimension = MESHDIMENSION
def getError(self):
funX = lambda x, y, z: np.cos(2*np.pi*y)
funY = lambda x, y, z: 1+x*0
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)
print Ec.shape, self.M.nE
print self.M
E = self.M.projectEdgeVector(Ec)
Fc = cartF3(self.M, solX, solY, solZ)
curlE_anal = self.M.projectFaceVector(Fc)
curlE = self.M.edgeCurl.dot(E)
err = np.linalg.norm((curlE - curlE_anal), np.inf)
return err
# def test_order(self):
# self.orderTest()
class TestCyl3DMesh(unittest.TestCase):