Source code for pyamg.gallery.diffusion
"""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).
-div Q A Q^T grad u
Q = [cos(theta) -sin(theta)]
[sin(theta) cos(theta)]
A = [1 0 ]
[0 eps ]
"""
import numpy
__docformat__ = "restructuredtext en"
__all__ = ['diffusion_stencil_2d']
[docs]def diffusion_stencil_2d(epsilon=1.0, theta=0.0, type='FE'):
"""
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)
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)
"""
C = numpy.cos(theta)
S = numpy.sin(theta)
CS = C*S
CC = C**2
SS = S**2
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 = numpy.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 ]
"""
C = numpy.cos(theta)
S = numpy.sin(theta)
CC = C**2
SS = S**2
CS = C*S
a = 0.5*(eps - 1)*CS
b = -(eps*SS + CC)
c = -a
d = -(eps*CC + SS)
e = 2.0*(eps + 1)
stencil = numpy.array([[a,b,c],
[d,e,d],
[c,b,a]])
return stencil