(mesh, kappa=None, f=None, degree=1)
return A, b
Finite element discretization of a Poisson problem.

    - div . kappa(x,y) grad u = f(x,y)

    Parameters
    ----------
    V : ndarray
        nv x 2 list of coordinates

    E : ndarray
        ne x 3 or 6 list of vertices

    kappa : function
        diffusion coefficient, kappa(x,y) with vector input

    f : function
        right hand side, f(x,y) with vector input

    degree : 1 or 2
        polynomial degree of the bases (assumed to be Lagrange locally)

    Returns
    -------
    A : sparse matrix
        finite element matrix where A_ij = <kappa grad phi_i, grad phi_j>

    b : array
        finite element rhs where b_ij = <f, phi_j>

    Notes
    -----
        - modepy is used to generate the quadrature points
          q = modepy.XiaoGimbutasSimplexQuadrature(4,2)

    Example
    -------
    >>> import numpy as np
    >>> from pyamg.gallery import fem
    >>> import scipy.sparse.linalg as sla
    >>> V = np.array(
    ... [[  0,  0],
    ...  [  1,  0],
    ...  [2*1,  0],
    ...  [  0,  1],
    ...  [  1,  1],
    ...  [2*1,  1],
    ...  [  0,2*1],
    ...  [  1,2*1],
    ...  [2*1,2*1]])
    >>> E = np.array(
    ... [[0, 1, 3],
    ...  [1, 2, 4],
    ...  [1, 4, 3],
    ...  [2, 5, 4],
    ...  [3, 4, 6],
    ...  [4, 5, 7],
    ...  [4, 7, 6],
    ...  [5, 8, 7]])
    >>> mesh = Mesh(V, E)
    >>> A, b = fem.gradgradform(mesh)
    >>> print(A.toarray())
    [[ 1.  -0.5  0.  -0.5  0.   0.   0.   0.   0. ]
     [-0.5  2.  -0.5  0.  -1.   0.   0.   0.   0. ]
     [ 0.  -0.5  1.   0.   0.  -0.5  0.   0.   0. ]
     [-0.5  0.   0.   2.  -1.   0.  -0.5  0.   0. ]
     [ 0.  -1.   0.  -1.   4.  -1.   0.  -1.   0. ]
     [ 0.   0.  -0.5  0.  -1.   2.   0.   0.  -0.5]
     [ 0.   0.   0.  -0.5  0.   0.   1.  -0.5  0. ]
     [ 0.   0.   0.   0.  -1.   0.  -0.5  2.  -0.5]
     [ 0.   0.   0.   0.   0.  -0.5  0.  -0.5  1. ]]
    >>> print(b)
    [0. 0. 0. 0. 0. 0. 0. 0. 0.]
    >>> f = lambda x, y : 0*x + 1.0
    >>> g = lambda x, y : 0*x + 0.0
    >>> g1 = lambda x, y : 0*x + 1.0
    >>> tol = 1e-12
    >>> X, Y = V[:,0], V[:,1]
    >>> id1 = np.where(abs(Y) < tol)[0]    # north
    >>> id2 = np.where(abs(Y-2) < tol)[0]  # south
    >>> bc = [{'id': id1, 'g': g},
    ...       {'id': id2, 'g': g}]
    >>> A, b = fem.gradgradform(mesh, f=f)
    >>> A, b = fem.applybc(A, b, mesh, bc)
    >>> A = A.tocsr()
    >>> u = sla.spsolve(A, b)
    >>> print(A.toarray())
    [[ 1.  0.  0.  0.  0.  0.  0.  0.  0.]
     [ 0.  1.  0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  1.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  2. -1.  0.  0.  0.  0.]
     [ 0.  0.  0. -1.  4. -1.  0.  0.  0.]
     [ 0.  0.  0.  0. -1.  2.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  1.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.  1.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.  0.  1.]]
    >>> print(b)
    [0.  0.  0.  0.5 1.  0.5 0.  0.  0. ]
    >>> print(u)
    [0.  0.  0.  0.5 0.5 0.5 0.  0.  0. ]
Finite element discretization of a Poisson problem.
Finite element discretization of a Poisson problem .
def gradgradform(mesh, kappa=None, f=None, degree=1):
    """Finite element discretization of a Poisson problem.

    - div . kappa(x,y) grad u = f(x,y)

    Parameters
    ----------
    V : ndarray
        nv x 2 list of coordinates

    E : ndarray
        ne x 3 or 6 list of vertices

    kappa : function
        diffusion coefficient, kappa(x,y) with vector input

    f : function
        right hand side, f(x,y) with vector input

    degree : 1 or 2
        polynomial degree of the bases (assumed to be Lagrange locally)

    Returns
    -------
    A : sparse matrix
        finite element matrix where A_ij = <kappa grad phi_i, grad phi_j>

    b : array
        finite element rhs where b_ij = <f, phi_j>

    Notes
    -----
        - modepy is used to generate the quadrature points
          q = modepy.XiaoGimbutasSimplexQuadrature(4,2)

    Example
    -------
    >>> import numpy as np
    >>> from pyamg.gallery import fem
    >>> import scipy.sparse.linalg as sla
    >>> V = np.array(
    ... [[  0,  0],
    ...  [  1,  0],
    ...  [2*1,  0],
    ...  [  0,  1],
    ...  [  1,  1],
    ...  [2*1,  1],
    ...  [  0,2*1],
    ...  [  1,2*1],
    ...  [2*1,2*1]])
    >>> E = np.array(
    ... [[0, 1, 3],
    ...  [1, 2, 4],
    ...  [1, 4, 3],
    ...  [2, 5, 4],
    ...  [3, 4, 6],
    ...  [4, 5, 7],
    ...  [4, 7, 6],
    ...  [5, 8, 7]])
    >>> mesh = Mesh(V, E)
    >>> A, b = fem.gradgradform(mesh)
    >>> print(A.toarray())
    [[ 1.  -0.5  0.  -0.5  0.   0.   0.   0.   0. ]
     [-0.5  2.  -0.5  0.  -1.   0.   0.   0.   0. ]
     [ 0.  -0.5  1.   0.   0.  -0.5  0.   0.   0. ]
     [-0.5  0.   0.   2.  -1.   0.  -0.5  0.   0. ]
     [ 0.  -1.   0.  -1.   4.  -1.   0.  -1.   0. ]
     [ 0.   0.  -0.5  0.  -1.   2.   0.   0.  -0.5]
     [ 0.   0.   0.  -0.5  0.   0.   1.  -0.5  0. ]
     [ 0.   0.   0.   0.  -1.   0.  -0.5  2.  -0.5]
     [ 0.   0.   0.   0.   0.  -0.5  0.  -0.5  1. ]]
    >>> print(b)
    [0. 0. 0. 0. 0. 0. 0. 0. 0.]
    >>> f = lambda x, y : 0*x + 1.0
    >>> g = lambda x, y : 0*x + 0.0
    >>> g1 = lambda x, y : 0*x + 1.0
    >>> tol = 1e-12
    >>> X, Y = V[:,0], V[:,1]
    >>> id1 = np.where(abs(Y) < tol)[0]    # north
    >>> id2 = np.where(abs(Y-2) < tol)[0]  # south
    >>> bc = [{'id': id1, 'g': g},
    ...       {'id': id2, 'g': g}]
    >>> A, b = fem.gradgradform(mesh, f=f)
    >>> A, b = fem.applybc(A, b, mesh, bc)
    >>> A = A.tocsr()
    >>> u = sla.spsolve(A, b)
    >>> print(A.toarray())
    [[ 1.  0.  0.  0.  0.  0.  0.  0.  0.]
     [ 0.  1.  0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  1.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  2. -1.  0.  0.  0.  0.]
     [ 0.  0.  0. -1.  4. -1.  0.  0.  0.]
     [ 0.  0.  0.  0. -1.  2.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  1.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.  1.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.  0.  1.]]
    >>> print(b)
    [0.  0.  0.  0.5 1.  0.5 0.  0.  0. ]
    >>> print(u)
    [0.  0.  0.  0.5 0.5 0.5 0.  0.  0. ]
    """
    if degree not in [1, 2]:
        raise ValueError('degree = 1 or 2 supported')

    if f is None:
        def f(_x, _y):
            return 0.0

    if kappa is None:
        def kappa(_x, _y):
            return 1.0

    if not callable(f) or not callable(kappa):
        raise ValueError('f, kappa must be callable functions')

    ne = mesh.ne

    if degree == 1:
        E = mesh.E
        X = mesh.X
        Y = mesh.Y

    if degree == 2:
        E = mesh.E2
        X = mesh.X2
        Y = mesh.Y2

    # allocate sparse matrix arrays
    m = 3 if degree == 1 else 6
    AA = np.zeros((ne, m**2))
    IA = np.zeros((ne, m**2), dtype=int)
    JA = np.zeros((ne, m**2), dtype=int)
    bb = np.zeros((ne, m))
    ib = np.zeros((ne, m), dtype=int)
    jb = np.zeros((ne, m), dtype=int)

    # Assemble A and b
    for ei in range(0, ne):
        # Step 1: set the vertices and indices
        K = E[ei, :]
        x0, y0 = X[K[0]], Y[K[0]]
        x1, y1 = X[K[1]], Y[K[1]]
        x2, y2 = X[K[2]], Y[K[2]]

        # Step 2: compute the Jacobian, inv, and det
        J = np.array([[x1 - x0, x2 - x0],
                      [y1 - y0, y2 - y0]])
        invJ = np.linalg.inv(J.T)
        detJ = np.linalg.det(J)

        if degree == 1:
            # Step 3, define the gradient of the basis
            dbasis = np.array([[-1, 1, 0],
                               [-1, 0, 1]])

            # Step 4
            dphi = invJ.dot(dbasis)

            # Step 5, 1-point gauss quadrature
            Aelem = kappa(X[K].mean(), Y[K].mean()) * (detJ / 2.0) * (dphi.T).dot(dphi)

            # Step 6, 1-point gauss quadrature
            belem = f(X[K].mean(), Y[K].mean()) * (detJ / 6.0) * np.ones((3,))

        if degree == 2:
            ww = np.array([0.44676317935602256, 0.44676317935602256, 0.44676317935602256,
                           0.21990348731064327, 0.21990348731064327, 0.21990348731064327])
            xy = np.array([[-0.10810301816807008, -0.78379396366385990],
                           [-0.10810301816806966, -0.10810301816807061],
                           [-0.78379396366386020, -0.10810301816806944],
                           [-0.81684757298045740, -0.81684757298045920],
                           [0.63369514596091700, -0.81684757298045810],
                           [-0.81684757298045870, 0.63369514596091750]])
            xx, yy = (xy[:, 0]+1)/2, (xy[:, 1]+1)/2
            ww *= 0.5

            Aelem = np.zeros((m, m))
            belem = np.zeros((m,))

            for w, x, y in zip(ww, xx, yy):
                # Step 3
                basis = np.array([(1-x-y)*(1-2*x-2*y),
                                  x*(2*x-1),
                                  y*(2*y-1),
                                  4*x*(1-x-y),
                                  4*x*y,
                                  4*y*(1-x-y)])

                dbasis = np.array([
                    [4*x + 4*y - 3, 4*x-1,     0, -8*x - 4*y + 4, 4*y,           -4*y],
                    [4*x + 4*y - 3,     0, 4*y-1,           -4*x, 4*x, -4*x - 8*y + 4]
                ])

                # Step 4
                dphi = invJ.dot(dbasis)

                # Step 5
                xt, yt = J.dot(np.array([x, y])) + np.array([x0, y0])
                Aelem += (detJ / 2) * w * kappa(xt, yt) * dphi.T.dot(dphi)

                # Step 6
                belem += (detJ / 2) * w * f(xt, yt) * basis

        # Step 7
        AA[ei, :] = Aelem.ravel()
        IA[ei, :] = np.repeat(K[np.arange(m)], m)
        JA[ei, :] = np.tile(K[np.arange(m)], m)
        bb[ei, :] = belem.ravel()
        ib[ei, :] = K[np.arange(m)]
        jb[ei, :] = 0

    # convert matrices
    A = sparse.coo_matrix((AA.ravel(), (IA.ravel(), JA.ravel())))
    A.sum_duplicates()
    b = sparse.coo_matrix((bb.ravel(), (ib.ravel(), jb.ravel()))).toarray().ravel()

    # A = A.tocsr()
    return A, b
def gradgradform ( mesh , kappa = None , f = None , degree = 1 ) : if degree not in [ 1 , 2 ] : raise ValueError ( 'degree = 1 or 2 supported' ) if f is None : def f ( _x , _y ) : return 0.0 if kappa is None : def kappa ( _x , _y ) : return 1.0 if not callable ( f ) or not callable ( kappa ) : raise ValueError ( 'f, kappa must be callable functions' ) ne = mesh . ne if degree == 1 : E = mesh . E X = mesh . X Y = mesh . Y if degree == 2 : E = mesh . E2 X = mesh . X2 Y = mesh . Y2 # allocate sparse matrix arrays m = 3 if degree == 1 else 6 AA = np . zeros ( ( ne , m ** 2 ) ) IA = np . zeros ( ( ne , m ** 2 ) , dtype = int ) JA = np . zeros ( ( ne , m ** 2 ) , dtype = int ) bb = np . zeros ( ( ne , m ) ) ib = np . zeros ( ( ne , m ) , dtype = int ) jb = np . zeros ( ( ne , m ) , dtype = int ) # Assemble A and b for ei in range ( 0 , ne ) : # Step 1: set the vertices and indices K = E [ ei , : ] x0 , y0 = X [ K [ 0 ] ] , Y [ K [ 0 ] ] x1 , y1 = X [ K [ 1 ] ] , Y [ K [ 1 ] ] x2 , y2 = X [ K [ 2 ] ] , Y [ K [ 2 ] ] # Step 2: compute the Jacobian, inv, and det J = np . array ( [ [ x1 - x0 , x2 - x0 ] , [ y1 - y0 , y2 - y0 ] ] ) invJ = np . linalg . inv ( J . T ) detJ = np . linalg . det ( J ) if degree == 1 : # Step 3, define the gradient of the basis dbasis = np . array ( [ [ - 1 , 1 , 0 ] , [ - 1 , 0 , 1 ] ] ) # Step 4 dphi = invJ . dot ( dbasis ) # Step 5, 1-point gauss quadrature Aelem = kappa ( X [ K ] . mean ( ) , Y [ K ] . mean ( ) ) * ( detJ / 2.0 ) * ( dphi . T ) . dot ( dphi ) # Step 6, 1-point gauss quadrature belem = f ( X [ K ] . mean ( ) , Y [ K ] . mean ( ) ) * ( detJ / 6.0 ) * np . ones ( ( 3 , ) ) if degree == 2 : ww = np . array ( [ 0.44676317935602256 , 0.44676317935602256 , 0.44676317935602256 , 0.21990348731064327 , 0.21990348731064327 , 0.21990348731064327 ] ) xy = np . array ( [ [ - 0.10810301816807008 , - 0.78379396366385990 ] , [ - 0.10810301816806966 , - 0.10810301816807061 ] , [ - 0.78379396366386020 , - 0.10810301816806944 ] , [ - 0.81684757298045740 , - 0.81684757298045920 ] , [ 0.63369514596091700 , - 0.81684757298045810 ] , [ - 0.81684757298045870 , 0.63369514596091750 ] ] ) xx , yy = ( xy [ : , 0 ] + 1 ) / 2 , ( xy [ : , 1 ] + 1 ) / 2 ww *= 0.5 Aelem = np . zeros ( ( m , m ) ) belem = np . zeros ( ( m , ) ) for w , x , y in zip ( ww , xx , yy ) : # Step 3 basis = np . array ( [ ( 1 - x - y ) * ( 1 - 2 * x - 2 * y ) , x * ( 2 * x - 1 ) , y * ( 2 * y - 1 ) , 4 * x * ( 1 - x - y ) , 4 * x * y , 4 * y * ( 1 - x - y ) ] ) dbasis = np . array ( [ [ 4 * x + 4 * y - 3 , 4 * x - 1 , 0 , - 8 * x - 4 * y + 4 , 4 * y , - 4 * y ] , [ 4 * x + 4 * y - 3 , 0 , 4 * y - 1 , - 4 * x , 4 * x , - 4 * x - 8 * y + 4 ] ] ) # Step 4 dphi = invJ . dot ( dbasis ) # Step 5 xt , yt = J . dot ( np . array ( [ x , y ] ) ) + np . array ( [ x0 , y0 ] ) Aelem += ( detJ / 2 ) * w * kappa ( xt , yt ) * dphi . T . dot ( dphi ) # Step 6 belem += ( detJ / 2 ) * w * f ( xt , yt ) * basis # Step 7 AA [ ei , : ] = Aelem . ravel ( ) IA [ ei , : ] = np . repeat ( K [ np . arange ( m ) ] , m ) JA [ ei , : ] = np . tile ( K [ np . arange ( m ) ] , m ) bb [ ei , : ] = belem . ravel ( ) ib [ ei , : ] = K [ np . arange ( m ) ] jb [ ei , : ] = 0 # convert matrices A = sparse . coo_matrix ( ( AA . ravel ( ) , ( IA . ravel ( ) , JA . ravel ( ) ) ) ) A . sum_duplicates ( ) b = sparse . coo_matrix ( ( bb . ravel ( ) , ( ib . ravel ( ) , jb . ravel ( ) ) ) ) . toarray ( ) . ravel ( ) # A = A.tocsr() return A , b