(epsilony=1.0, epsilonz=1.0, theta=0.0, phi=0.0,
                         psi=0.0, type='FD')
return stencil
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 as sp
    >>> from pyamg.gallery.diffusion import diffusion_stencil_2d
    >>> sten = diffusion_stencil_2d(epsilon=0.0001,theta=sp.pi/6,type='FD')
    >>> print(sten)
    [[-0.2164847 -0.750025   0.2164847]
     [-0.250075   2.0002    -0.250075 ]
     [ 0.2164847 -0.750025  -0.2164847]]
Rotated Anisotropic diffusion in 3d of the form.
Rotated Anisotropic diffusion in 3d of the form .
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 as sp
    >>> from pyamg.gallery.diffusion import diffusion_stencil_2d
    >>> sten = diffusion_stencil_2d(epsilon=0.0001,theta=sp.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 == 'FE':
        raise NotImplementedError('FE not implemented yet')

    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
def diffusion_stencil_3d ( epsilony = 1.0 , epsilonz = 1.0 , theta = 0.0 , phi = 0.0 , psi = 0.0 , type = 'FD' ) : 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 == 'FE' : raise NotImplementedError ( 'FE not implemented yet' ) 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