"""Generate a diffusion stencil
Supports isotropic diffusion (FE,FD), anisotropic diffusion (FE, FD), and
rotated anisotropic diffusion (FD).
The stencils include redundancy to maintain readability for simple cases (e.g.
isotropic diffusion).
"""
import numpy as np
__docformat__ = "restructuredtext en"
__all__ = ['diffusion_stencil_2d']
[docs]def diffusion_stencil_2d(epsilon=1.0, theta=0.0, type='FE'):
"""
Rotated Anisotropic diffusion in 2d of the form:
-div Q A Q^T grad u
Q = [cos(theta) -sin(theta)]
[sin(theta) cos(theta)]
A = [1 0 ]
[0 eps ]
Parameters
----------
epsilon : float, optional
Anisotropic diffusion coefficient: -div A grad u,
where A = [1 0; 0 epsilon]. The default is isotropic, epsilon=1.0
theta : float, optional
Rotation angle `theta` in radians defines -div Q A Q^T grad,
where Q = [cos(`theta`) -sin(`theta`); sin(`theta`) cos(`theta`)].
type : {'FE','FD'}
Specifies the discretization as Q1 finite element (FE) or 2nd order
finite difference (FD)
The default is `theta` = 0.0
Returns
-------
stencil : numpy array
A 3x3 diffusion stencil
See Also
--------
stencil_grid, poisson
Notes
-----
Not all combinations are supported.
Examples
--------
>>> import scipy
>>> from pyamg.gallery.diffusion import diffusion_stencil_2d
>>> sten = diffusion_stencil_2d(epsilon=0.0001,theta=scipy.pi/6,type='FD')
>>> print sten
[[-0.2164847 -0.750025 0.2164847]
[-0.250075 2.0002 -0.250075 ]
[ 0.2164847 -0.750025 -0.2164847]]
"""
eps = float(epsilon) # for brevity
theta = float(theta)
C = np.cos(theta)
S = np.sin(theta)
CS = C*S
CC = C**2
SS = S**2
if(type == 'FE'):
"""FE approximation to::
- (eps c^2 + s^2) u_xx +
-2(eps - 1) c s u_xy +
- ( c^2 + eps s^2) u_yy
[ -c^2*eps-s^2+3*c*s*(eps-1)-c^2-s^2*eps,
2*c^2*eps+2*s^2-4*c^2-4*s^2*eps,
-c^2*eps-s^2-3*c*s*(eps-1)-c^2-s^2*eps]
[-4*c^2*eps-4*s^2+2*c^2+2*s^2*eps,
8*c^2*eps+8*s^2+8*c^2+8*s^2*eps,
-4*c^2*eps-4*s^2+2*c^2+2*s^2*eps]
[-c^2*eps-s^2-3*c*s*(eps-1)-c^2-s^2*eps,
2*c^2*eps+2*s^2-4*c^2-4*s^2*eps,
-c^2*eps-s^2+3*c*s*(eps-1)-c^2-s^2*eps]
c = cos(theta)
s = sin(theta)
"""
a = (-1*eps - 1)*CC + (-1*eps - 1)*SS + (3*eps - 3)*CS
b = (2*eps - 4)*CC + (-4*eps + 2)*SS
c = (-1*eps - 1)*CC + (-1*eps - 1)*SS + (-3*eps + 3)*CS
d = (-4*eps + 2)*CC + (2*eps - 4)*SS
e = (8*eps + 8)*CC + (8*eps + 8)*SS
stencil = np.array([[a, b, c],
[d, e, d],
[c, b, a]]) / 6.0
elif type == 'FD':
"""FD approximation to:
- (eps c^2 + s^2) u_xx +
-2(eps - 1) c s u_xy +
- ( c^2 + eps s^2) u_yy
c = cos(theta)
s = sin(theta)
A = [ 1/2(eps - 1) c s -(c^2 + eps s^2) -1/2(eps - 1) c s ]
[ ]
[ -(eps c^2 + s^2) 2 (eps + 1) -(eps c^2 + s^2) ]
[ ]
[ -1/2(eps - 1) c s -(c^2 + eps s^2) 1/2(eps - 1) c s ]
"""
a = 0.5*(eps - 1)*CS
b = -(eps*SS + CC)
c = -a
d = -(eps*CC + SS)
e = 2.0*(eps + 1)
stencil = np.array([[a, b, c],
[d, e, d],
[c, b, a]])
return stencil
def _symbolic_rotation_helper():
"""
Simple SymPy script to generate the 3D rotation matrix and products for
diffusion_stencil_3d.
"""
from sympy import symbols, Matrix
cpsi, spsi = symbols('cpsi, spsi')
cth, sth = symbols('cth, sth')
cphi, sphi = symbols('cphi, sphi')
Rpsi = Matrix([[cpsi, spsi, 0], [-spsi, cpsi, 0], [0, 0, 1]])
Rth = Matrix([[1, 0, 0], [0, cth, sth], [0, -sth, cth]])
Rphi = Matrix([[cphi, sphi, 0], [-sphi, cphi, 0], [0, 0, 1]])
Q = Rpsi * Rth * Rphi
epsy, epsz = symbols('epsy, epsz')
A = Matrix([[1, 0, 0], [0, epsy, 0], [0, 0, epsz]])
D = Q * A * Q.T
for i in range(3):
for j in range(3):
print 'D[%d, %d] = %s' % (i, j, D[i, j])
def _symbolic_product_helper():
"""
Simple SymPy script to generate the 3D products for diffusion_stencil_3d.
"""
from sympy import symbols, Matrix
D11, D12, D13, D21, D22, D23, D31, D32, D33 =\
symbols('D11, D12, D13, D21, D22, D23, D31, D32, D33')
D = Matrix([[D11, D12, D13], [D21, D22, D23], [D31, D32, D33]])
grad = Matrix([['dx', 'dy', 'dz']]).T
div = grad.T
a = div * D * grad
print a[0]
def diffusion_stencil_3d(epsilony=1.0, epsilonz=1.0, theta=0.0, phi=0.0,
psi=0.0, type='FD'):
"""
Rotated Anisotropic diffusion in 3d of the form:
-div Q A Q^T grad u
Q = Rpsi Rtheta Rphi
Rpsi = [ c s 0 ]
[-s c 0 ]
[ 0 0 1 ]
c = cos(psi)
s = sin(psi)
Rtheta = [ 1 0 0 ]
[ 0 c s ]
[ 0 -s c ]
c = cos(theta)
s = sin(theta)
Rphi = [ c s 0 ]
[-s c 0 ]
[ 0 0 1 ]
c = cos(phi)
s = sin(phi)
Here Euler Angles are used:
http://en.wikipedia.org/wiki/Euler_angles
This results in
Q = [ cphi*cpsi - cth*sphi*spsi, cpsi*sphi + cphi*cth*spsi, spsi*sth]
[ - cphi*spsi - cpsi*cth*sphi, cphi*cpsi*cth - sphi*spsi, cpsi*sth]
[ sphi*sth, -cphi*sth, cth]
A = [1 0 ]
[0 epsy ]
[0 0 epsz]
D = Q A Q^T
Parameters
----------
epsilony : float, optional
Anisotropic diffusion coefficient in the y-direction
where A = [1 0 0; 0 epsilon_y 0; 0 0 epsilon_z]. The default is
isotropic, epsilon=1.0
epsilonz : float, optional
Anisotropic diffusion coefficient in the z-direction
where A = [1 0 0; 0 epsilon_y 0; 0 0 epsilon_z]. The default is
isotropic, epsilon=1.0
theta : float, optional
Euler rotation angle `theta` in radians. The default is 0.0.
phi : float, optional
Euler rotation angle `phi` in radians. The default is 0.0.
psi : float, optional
Euler rotation angle `psi` in radians. The default is 0.0.
type : {'FE','FD'}
Specifies the discretization as Q1 finite element (FE) or 2nd order
finite difference (FD)
Returns
-------
stencil : numpy array
A 3x3 diffusion stencil
See Also
--------
stencil_grid, poisson, _symbolic_rotation_helper, _symbolic_product_helper
Notes
-----
Not all combinations are supported.
Examples
--------
>>> import scipy
>>> from pyamg.gallery.diffusion import diffusion_stencil_2d
>>> sten = diffusion_stencil_2d(epsilon=0.0001,theta=scipy.pi/6,type='FD')
>>> print sten
[[-0.2164847 -0.750025 0.2164847]
[-0.250075 2.0002 -0.250075 ]
[ 0.2164847 -0.750025 -0.2164847]]
"""
epsy = float(epsilony) # for brevity
epsz = float(epsilonz) # for brevity
theta = float(theta)
phi = float(phi)
psi = float(psi)
D = np.zeros((3, 3))
cphi = np.cos(phi)
sphi = np.sin(phi)
cth = np.cos(theta)
sth = np.sin(theta)
cpsi = np.cos(psi)
spsi = np.sin(psi)
# from _symbolic_rotation_helper
D[0, 0] = epsy*(cphi*cth*spsi + cpsi*sphi)**2 + epsz*spsi**2*sth**2 +\
(cphi*cpsi - cth*sphi*spsi)**2
D[0, 1] = cpsi*epsz*spsi*sth**2 +\
epsy*(cphi*cpsi*cth - sphi*spsi)*(cphi*cth*spsi + cpsi*sphi) +\
(cphi*cpsi - cth*sphi*spsi)*(-cphi*spsi - cpsi*cth*sphi)
D[0, 2] = -cphi*epsy*sth*(cphi*cth*spsi + cpsi*sphi) +\
cth*epsz*spsi*sth + sphi*sth*(cphi*cpsi - cth*sphi*spsi)
D[1, 0] = cpsi*epsz*spsi*sth**2 +\
epsy*(cphi*cpsi*cth - sphi*spsi)*(cphi*cth*spsi + cpsi*sphi) +\
(cphi*cpsi - cth*sphi*spsi)*(-cphi*spsi - cpsi*cth*sphi)
D[1, 1] = cpsi**2*epsz*sth**2 + epsy*(cphi*cpsi*cth - sphi*spsi)**2 +\
(-cphi*spsi - cpsi*cth*sphi)**2
D[1, 2] = -cphi*epsy*sth*(cphi*cpsi*cth - sphi*spsi) +\
cpsi*cth*epsz*sth + sphi*sth*(-cphi*spsi - cpsi*cth*sphi)
D[2, 0] = -cphi*epsy*sth*(cphi*cth*spsi + cpsi*sphi) + cth*epsz*spsi*sth +\
sphi*sth*(cphi*cpsi - cth*sphi*spsi)
D[2, 1] = -cphi*epsy*sth*(cphi*cpsi*cth - sphi*spsi) + cpsi*cth*epsz*sth +\
sphi*sth*(-cphi*spsi - cpsi*cth*sphi)
D[2, 2] = cphi**2*epsy*sth**2 + cth**2*epsz + sphi**2*sth**2
stencil = np.zeros((3, 3, 3))
if type == 'FD':
# from _symbolic_product_helper
# dx*(D11*dx + D21*dy + D31*dz) +
# dy*(D12*dx + D22*dy + D32*dz) +
# dz*(D13*dx + D23*dy + D33*dz)
#
# D00*dxx +
# (D10+D01)*dxy +
# (D20+D02)*dxz +
# D11*dyy +
# (D21+D12)*dyz +
# D22*dzz
i, j, k = (1, 1, 1)
# dxx
stencil[[i-1, i, i+1], j, k] += np.array([-1, 2, -1]) * D[0, 0]
# dyy
stencil[i, [j-1, j, j+1], k] += np.array([-1, 2, -1]) * D[1, 1]
# dzz
stencil[i, j, [k-1, k, k+1]] += np.array([-1, 2, -1]) * D[2, 2]
L = np.array([-1, -1, 1, 1])
M = np.array([-1, 1, -1, 1])
# dxy
stencil[i + L, j + M, k] \
+= 0.25 * np.array([1, -1, -1, 1]) * (D[1, 0] + D[0, 1])
# dxz
stencil[i + L, j, k + M] \
+= 0.25 * np.array([1, -1, -1, 1]) * (D[2, 0] + D[0, 2])
# dyz
stencil[i, j + L, k + M] \
+= 0.25 * np.array([1, -1, -1, 1]) * (D[2, 1] + D[1, 2])
return stencil
if type == 'FE':
raise NotImplementedError('FE not implemented yet')