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

changes to CC at middle of first h_r

Browse Source
This commit is contained in:
rowanc1
2014-02-21 15:31:31 -08:00
parent 897b43bf80
commit fee512dda3
3 changed files with 24 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/CylMesh.py
+1 -2
View File
@@ -89,8 +89,7 @@ class CylMesh(BaseTensorMesh):
@property
def vectorCCx(self):
"""Cell-centered grid vector (1D) in the x direction."""
firstEl = -self.hx[0]*0.5 if self.nCy == 1 else 0
return np.r_[firstEl, self.hx[:-1].cumsum()] + self.hx*0.5
return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5
@property
def vectorCCy(self):
SimPEG/Mesh/TensorMesh.py
+1 -1
View File
@@ -22,7 +22,7 @@ class BaseTensorMesh(BaseRectangularMesh):
assert type(h_in) is list, 'h_in must be a list'
h = range(len(h_in))
for i, h_i in enumerate(h_in):
if type(h_i) in [int, long, float]:
if type(h_i) in [int, long, float, np.int_]:
# This gives you something over the unit cube.
h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i)
assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
SimPEG/Tests/test_cylMesh.py
+22 -26
View File
@@ -49,7 +49,7 @@ class TestCyl2DMesh(unittest.TestCase):
self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0]))
def test_vectorsCC(self):
v = np.r_[0, 1.5, 2.25]
v = np.r_[0.5, 1.5, 2.25]
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0)
v = np.r_[0]
self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0)
@@ -96,7 +96,7 @@ class TestCyl2DMesh(unittest.TestCase):
self.assertTrue(self.mesh.gridEz is None)
def test_gridCC(self):
x = np.r_[0,1.5,2.25,0,1.5,2.25]
x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25]
y = np.zeros(6)
z = np.r_[1,1,1,2.5,2.5,2.5]
G = np.c_[x,y,z]
@@ -117,7 +117,7 @@ class TestCyl2DMesh(unittest.TestCase):
self.assertTrue(np.linalg.norm((G-self.mesh.gridFx).ravel()) == 0)
def test_gridFz(self):
x = np.r_[0,1.5,2.25,0,1.5,2.25,0,1.5,2.25]
x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25,0.5,1.5,2.25]
y = np.zeros(9)
z = np.r_[0,0,0,2,2,2,3,3,3.]
G = np.c_[x,y,z]
@@ -137,14 +137,14 @@ class TestCyl2DMesh(unittest.TestCase):
MESHTYPES = ['uniformCylMesh']
MESHDIMENSION = 2
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 2])
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)))
cyl_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
cyl_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
cylF2 = lambda M, fx, fy: np.vstack((cyl_row2(M.gridFx, fx, fy), cyl_row2(M.gridFz, fx, fy)))
# cylE2 = lambda M, ex, ey: np.vstack((cyl_row2(M.gridEx, ex, ey), cyl_row2(M.gridEy, ex, ey)))
# cylF3 = lambda M, fx, fy, fz: np.vstack((cyl_row3(M.gridFx, fx, fy, fz), cyl_row3(M.gridFy, fx, fy, fz), cyl_row3(M.gridFz, fx, fy, fz)))
# cylE3 = lambda M, ex, ey, ez: np.vstack((cyl_row3(M.gridEx, ex, ey, ez), cyl_row3(M.gridEy, ex, ey, ez), cyl_row3(M.gridEz, ex, ey, ez)))
class TestFaceDiv(OrderTest):
@@ -154,28 +154,24 @@ class TestFaceDiv(OrderTest):
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)
funR = lambda r, z: np.sin(2.*np.pi*r)
funZ = lambda r, z: np.sin(2.*np.pi*z)
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)
sol = lambda r, t, z: (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r + 2*np.pi*np.cos(2*np.pi*z)
Ec = cartE3(self.M, funX, funY, funZ)
print Ec.shape, self.M.nE
print self.M
E = self.M.projectEdgeVector(Ec)
Fc = cylF2(self.M, funR, funZ)
Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
F = self.M.projectFaceVector(Fc)
Fc = cartF3(self.M, solX, solY, solZ)
curlE_anal = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
tm = Mesh.TensorMesh([self.M.nCx, self.M.nCz])
divF_anal = call3(sol, self.M.gridCC)
curlE = self.M.edgeCurl.dot(E)
err = np.linalg.norm((curlE - curlE_anal), np.inf)
err = np.linalg.norm((divF-divF_anal), np.inf)
return err
# def test_order(self):
# self.orderTest()
def test_order(self):
self.orderTest()
class TestCyl3DMesh(unittest.TestCase):