mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-02 10:12:29 +08:00
Merge branch 'boundaryConditions' of https://bitbucket.org/rcockett/simpeg into richards
This commit is contained in:
Rowan Cockett
@@ -4,3 +4,8 @@ import mesh
|
||||
import inverse
|
||||
import forward
|
||||
import regularization
|
||||
|
||||
|
||||
import scipy.version as _v
|
||||
if _v.version < '0.13.0':
|
||||
print 'Warning: upgrade your scipy to 0.13.0'
|
||||
|
||||
+152
-21
@@ -1,10 +1,16 @@
|
||||
import numpy as np
|
||||
from scipy import sparse as sp
|
||||
from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av
|
||||
from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap
|
||||
|
||||
|
||||
def checkBC(bc):
|
||||
""" Checks if boundary condition 'bc' is valid. """
|
||||
"""
|
||||
|
||||
Checks if boundary condition 'bc' is valid.
|
||||
|
||||
Each bc must be either 'dirichlet' or 'neumann'
|
||||
|
||||
"""
|
||||
if(type(bc) is str):
|
||||
bc = [bc, bc]
|
||||
assert type(bc) is list, 'bc must be a list'
|
||||
@@ -17,7 +23,33 @@ def checkBC(bc):
|
||||
|
||||
|
||||
def ddxCellGrad(n, bc):
|
||||
"""Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1"""
|
||||
"""
|
||||
Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
|
||||
|
||||
For Cell-Centered **Dirichlet**, use a ghost point::
|
||||
|
||||
(u_1 - u_g)/hf = grad
|
||||
|
||||
u_g u_1 u_2
|
||||
* | * | * ...
|
||||
^
|
||||
0
|
||||
|
||||
u_g = - u_1
|
||||
grad = 2*u1/dx
|
||||
negitive on the other side.
|
||||
|
||||
For Cell-Centered **Neumann**, use a ghost point::
|
||||
|
||||
(u_1 - u_g)/hf = 0
|
||||
|
||||
u_g u_1 u_2
|
||||
* | * | * ...
|
||||
|
||||
u_g = u_1
|
||||
grad = 0; put a zero in.
|
||||
|
||||
"""
|
||||
bc = checkBC(bc)
|
||||
|
||||
D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
|
||||
@@ -33,6 +65,56 @@ def ddxCellGrad(n, bc):
|
||||
D[-1, -1] = 0
|
||||
return D
|
||||
|
||||
def ddxCellGradBC(n, bc):
|
||||
"""
|
||||
|
||||
Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
|
||||
|
||||
For Cell-Centered **Dirichlet**, use a ghost point::
|
||||
|
||||
(u_1 - u_g)/hf = grad
|
||||
|
||||
u_g u_1 u_2
|
||||
* | * | * ...
|
||||
^
|
||||
u_b
|
||||
|
||||
We know the value at the boundary (u_b)::
|
||||
|
||||
(u_g+u_1)/2 = u_b (the average)
|
||||
u_g = 2*u_b - u_1
|
||||
|
||||
So plug in to gradient:
|
||||
|
||||
(u_1 - (2*u_b - u_1))/hf = grad
|
||||
2*(u_1-u_b)/hf = grad
|
||||
|
||||
Separate, because BC are known (and can move to RHS later)::
|
||||
|
||||
( 2/hf )*u_1 + ( -2/hf )*u_b = grad
|
||||
|
||||
( ^ ) JUST RETURN THIS
|
||||
|
||||
|
||||
"""
|
||||
bc = checkBC(bc)
|
||||
|
||||
ij = (np.array([0, n]),np.array([0, 1]))
|
||||
vals = np.zeros(2)
|
||||
|
||||
# Set the first side
|
||||
if(bc[0] == 'dirichlet'):
|
||||
vals[0] = -2
|
||||
elif(bc[0] == 'neumann'):
|
||||
vals[0] = 0
|
||||
# Set the second side
|
||||
if(bc[1] == 'dirichlet'):
|
||||
vals[1] = 2
|
||||
elif(bc[1] == 'neumann'):
|
||||
vals[1] = 0
|
||||
D = sp.csr_matrix((vals, ij), shape=(n+1,2))
|
||||
return D
|
||||
|
||||
|
||||
class DiffOperators(object):
|
||||
"""
|
||||
@@ -151,17 +233,19 @@ class DiffOperators(object):
|
||||
for i, bc_i in enumerate(BC):
|
||||
BC[i] = checkBC(bc_i)
|
||||
|
||||
self._cellGrad = None # ensure we create a new gradient next time we call it
|
||||
self._cellGradBC = BC
|
||||
# ensure we create a new gradient next time we call it
|
||||
self._cellGrad = None
|
||||
self._cellGradBC = None
|
||||
self._cellGradBC_list = BC
|
||||
return BC
|
||||
_cellGradBC = 'neumann'
|
||||
_cellGradBC_list = 'neumann'
|
||||
|
||||
def cellGrad():
|
||||
doc = "The cell centered Gradient, takes you to cell faces."
|
||||
|
||||
def fget(self):
|
||||
if(self._cellGrad is None):
|
||||
BC = self.setCellGradBC(self._cellGradBC)
|
||||
BC = self.setCellGradBC(self._cellGradBC_list)
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
G = ddxCellGrad(n[0], BC[0])
|
||||
@@ -176,13 +260,40 @@ class DiffOperators(object):
|
||||
G = sp.vstack((G1, G2, G3), format="csr")
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.area
|
||||
V = self.vol
|
||||
self._cellGrad = sdiag(S)*G*sdiag(1/V)
|
||||
V = self.aveCC2F*self.vol # Average volume between adjacent cells
|
||||
self._cellGrad = sdiag(S/V)*G
|
||||
return self._cellGrad
|
||||
return locals()
|
||||
_cellGrad = None
|
||||
cellGrad = property(**cellGrad())
|
||||
|
||||
def cellGradBC():
|
||||
doc = "The cell centered Gradient boundary condition matrix"
|
||||
|
||||
def fget(self):
|
||||
if(self._cellGradBC is None):
|
||||
BC = self.setCellGradBC(self._cellGradBC_list)
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
G = ddxCellGradBC(n[0], BC[0])
|
||||
elif(self.dim == 2):
|
||||
G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
|
||||
G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
|
||||
G = sp.vstack((G1, G2), format="csr")
|
||||
elif(self.dim == 3):
|
||||
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
|
||||
G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
|
||||
G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
G = sp.vstack((G1, G2, G3), format="csr")
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.area
|
||||
V = self.aveCC2F*self.vol # Average volume between adjacent cells
|
||||
self._cellGradBC = sdiag(S/V)*G
|
||||
return self._cellGradBC
|
||||
return locals()
|
||||
_cellGradBC = None
|
||||
cellGradBC = property(**cellGradBC())
|
||||
|
||||
def cellGradx():
|
||||
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
|
||||
|
||||
@@ -197,9 +308,9 @@ class DiffOperators(object):
|
||||
elif(self.dim == 3):
|
||||
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.r(self.area, 'F','Fx', 'V')
|
||||
V = self.vol
|
||||
self._cellGradx = sdiag(S)*G1*sdiag(1/V)
|
||||
V = self.aveCC2F*self.vol
|
||||
L = self.r(self.area/V, 'F','Fx', 'V')
|
||||
self._cellGradx = sdiag(L)*G1
|
||||
return self._cellGradx
|
||||
return locals()
|
||||
cellGradx = property(**cellGradx())
|
||||
@@ -217,9 +328,9 @@ class DiffOperators(object):
|
||||
elif(self.dim == 3):
|
||||
G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.r(self.area, 'F','Fy', 'V')
|
||||
V = self.vol
|
||||
self._cellGrady = sdiag(S)*G2*sdiag(1/V)
|
||||
V = self.aveCC2F*self.vol
|
||||
L = self.r(self.area/V, 'F','Fy', 'V')
|
||||
self._cellGrady = sdiag(L)*G2
|
||||
return self._cellGrady
|
||||
return locals()
|
||||
cellGrady = property(**cellGrady())
|
||||
@@ -234,9 +345,9 @@ class DiffOperators(object):
|
||||
n = self.n
|
||||
G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.r(self.area, 'F','Fz', 'V')
|
||||
V = self.vol
|
||||
self._cellGradz = sdiag(S)*G3*sdiag(1/V)
|
||||
V = self.aveCC2F*self.vol
|
||||
L = self.r(self.area/V, 'F','Fz', 'V')
|
||||
self._cellGradz = sdiag(L)*G3
|
||||
return self._cellGradz
|
||||
return locals()
|
||||
cellGradz = property(**cellGradz())
|
||||
@@ -295,16 +406,36 @@ class DiffOperators(object):
|
||||
self._aveF2CC = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
|
||||
sp.kron(av(n[1]), speye(n[0]))), format="csr")
|
||||
sp.kron(av(n[1]), speye(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._aveF2CC
|
||||
return locals()
|
||||
_aveF2CC = None
|
||||
aveF2CC = property(**aveF2CC())
|
||||
|
||||
def aveCC2F():
|
||||
doc = "Construct the averaging operator on cell cell centers to faces."
|
||||
|
||||
def fget(self):
|
||||
if(self._aveCC2F is None):
|
||||
n = self.n
|
||||
if(self.dim == 1):
|
||||
self._aveCC2F = avExtrap(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
|
||||
sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
|
||||
kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
|
||||
kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._aveCC2F
|
||||
return locals()
|
||||
_aveCC2F = None
|
||||
aveCC2F = property(**aveCC2F())
|
||||
|
||||
def aveE2CC():
|
||||
doc = "Construct the averaging operator on cell edges to cell centers."
|
||||
|
||||
|
||||
@@ -26,18 +26,25 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
|
||||
|
||||
.. plot:: examples/mesh/plot_TensorMesh.py
|
||||
|
||||
For a quick tensor mesh on a (10x12x15) unit cube::
|
||||
|
||||
mesh = TensorMesh([10, 12, 15])
|
||||
|
||||
"""
|
||||
_meshType = 'TENSOR'
|
||||
|
||||
def __init__(self, h, x0=None):
|
||||
super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0)
|
||||
|
||||
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
|
||||
|
||||
for i, h_i in enumerate(h):
|
||||
if type(h_i) in [int, long, float]:
|
||||
# This gives you something over the unit cube.
|
||||
h_i = np.ones(int(h_i))/int(h_i)
|
||||
h[i] = h_i
|
||||
assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
|
||||
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
|
||||
|
||||
BaseMesh.__init__(self, np.array([x.size for x in h]), x0)
|
||||
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
|
||||
|
||||
# Ensure h contains 1D vectors
|
||||
self._h = [mkvc(x.astype(float)) for x in h]
|
||||
|
||||
|
||||
@@ -89,18 +89,18 @@ class TensorView(object):
|
||||
# Determine the subplot number: 131, 121
|
||||
numPlots = 130 if plotAll else len(imageType)/2*10+100
|
||||
pltNum = 1
|
||||
fx, fy, fz = self.r(I,'F','F','M')
|
||||
fxyz = self.r(I,'F','F','M')
|
||||
if plotAll or 'Fx' in imageType:
|
||||
ax_x = plt.subplot(numPlots+pltNum)
|
||||
self.plotImage(fx, imageType='Fx', ax=ax_x, **options)
|
||||
self.plotImage(fxyz[0], imageType='Fx', ax=ax_x, **options)
|
||||
pltNum +=1
|
||||
if plotAll or 'Fy' in imageType:
|
||||
ax_y = plt.subplot(numPlots+pltNum)
|
||||
self.plotImage(fy, imageType='Fy', ax=ax_y, **options)
|
||||
self.plotImage(fxyz[1], imageType='Fy', ax=ax_y, **options)
|
||||
pltNum +=1
|
||||
if plotAll or 'Fz' in imageType:
|
||||
ax_z = plt.subplot(numPlots+pltNum)
|
||||
self.plotImage(fz, imageType='Fz', ax=ax_z, **options)
|
||||
self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
|
||||
pltNum +=1
|
||||
return
|
||||
else:
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from TestUtils import OrderTest
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
|
||||
call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
|
||||
@@ -48,8 +49,148 @@ class TestCurl(OrderTest):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestFaceDiv(OrderTest):
|
||||
name = "Face Divergence"
|
||||
class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
|
||||
name = "Cell Grad 1D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 1
|
||||
expectedOrders = 1 # because of the averaging involved in the ghost point. u_b = (u_n + u_g)/2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x: -2*np.pi*np.sin(2*np.pi*x)
|
||||
sol = lambda x: np.cos(2*np.pi*x)
|
||||
|
||||
|
||||
xc = sol(self.M.gridCC)
|
||||
|
||||
gradX_anal = fx(self.M.gridFx)
|
||||
|
||||
bc = np.array([1,1])
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
gradX = self.M.cellGrad.dot(xc) + self.M.cellGradBC*bc
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestCellGrad2D_Dirichlet(OrderTest):
|
||||
name = "Cell Grad 2D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)
|
||||
fy = lambda x, y: 2*np.pi*np.cos(2*np.pi*y)*np.sin(2*np.pi*x)
|
||||
sol = lambda x, y: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)
|
||||
|
||||
xc = call2(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF2(self.M, fx, fy)
|
||||
gradX_anal = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestCellGrad3D_Dirichlet(OrderTest):
|
||||
name = "Cell Grad 3D - Dirichlet"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 3
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y, z: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
fy = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
fz = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
sol = lambda x, y, z: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
|
||||
xc = call3(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF3(self.M, fx, fy, fz)
|
||||
gradX_anal = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('dirichlet')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestCellGrad2D_Neumann(OrderTest):
|
||||
name = "Cell Grad 2D - Neumann"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 2
|
||||
meshSizes = [8, 16, 32, 64]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)
|
||||
fy = lambda x, y: -2*np.pi*np.sin(2*np.pi*y)*np.cos(2*np.pi*x)
|
||||
sol = lambda x, y: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)
|
||||
|
||||
xc = call2(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF2(self.M, fx, fy)
|
||||
gradX_anal = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('neumann')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestCellGrad3D_Neumann(OrderTest):
|
||||
name = "Cell Grad 3D - Neumann"
|
||||
meshTypes = ['uniformTensorMesh']
|
||||
meshDimension = 3
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
fx = lambda x, y, z: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
fy = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
fz = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z)
|
||||
sol = lambda x, y, z: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z)
|
||||
|
||||
xc = call3(sol, self.M.gridCC)
|
||||
|
||||
Fc = cartF3(self.M, fx, fy, fz)
|
||||
gradX_anal = self.M.projectFaceVector(Fc)
|
||||
|
||||
self.M.setCellGradBC('neumann')
|
||||
gradX = self.M.cellGrad.dot(xc)
|
||||
|
||||
err = np.linalg.norm((gradX-gradX_anal), np.inf)
|
||||
|
||||
return err
|
||||
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
class TestFaceDiv3D(OrderTest):
|
||||
name = "Face Divergence 3D"
|
||||
meshTypes = MESHTYPES
|
||||
meshSizes = [8, 16, 32]
|
||||
|
||||
@@ -197,6 +338,16 @@ class TestAveraging2D(OrderTest):
|
||||
self.getAve = lambda M: M.aveF2CC
|
||||
self.orderTest()
|
||||
|
||||
def test_orderCC2F(self):
|
||||
self.name = "Averaging 2D: CC2F"
|
||||
fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
self.getHere = lambda M: call2(fun, M.gridCC)
|
||||
self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
|
||||
self.getAve = lambda M: M.aveCC2F
|
||||
self.expectedOrders = 1
|
||||
self.orderTest()
|
||||
self.expectedOrders = 2
|
||||
|
||||
|
||||
def test_orderE2CC(self):
|
||||
self.name = "Averaging 2D: E2CC"
|
||||
@@ -258,6 +409,17 @@ class TestAveraging3D(OrderTest):
|
||||
self.getAve = lambda M: M.aveE2CC
|
||||
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()
|
||||
|
||||
@@ -4,7 +4,7 @@ import lomutils
|
||||
import interputils
|
||||
import ModelBuilder
|
||||
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
|
||||
from sputils import spzeros, kron3, speye, sdiag, ddx, av
|
||||
from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap
|
||||
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
|
||||
from interputils import interpmat
|
||||
from ipythonUtils import easyAnimate as animate
|
||||
|
||||
@@ -29,5 +29,10 @@ def ddx(n):
|
||||
|
||||
|
||||
def av(n):
|
||||
"""Define 1D averaging operator from cell-centres to nodes."""
|
||||
"""Define 1D averaging operator from nodes to cell-centers."""
|
||||
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
|
||||
|
||||
def avExtrap(n):
|
||||
"""Define 1D averaging operator from cell-centers to nodes."""
|
||||
Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
|
||||
return Av
|
||||
|
||||
Reference in New Issue
Block a user