(mesh)
return BX, BY
Calculate the (div u , p) form that arises in Stokes assumes P2-P1 elements.

    Parameters
    ----------
    mesh : Mesh
        Mesh object

    Returns
    -------
    BX, BY : ndarray
        div block B = [BX, BY].T in [[A, B], [B.T 0]]
Calculate the (div u , p) form that arises in Stokes assumes P2-P1 elements.
Calculate the ( div u p ) form that arises in Stokes assumes P2 - P1 elements .
def divform(mesh):
    """Calculate the (div u , p) form that arises in Stokes assumes P2-P1 elements.

    Parameters
    ----------
    mesh : Mesh
        Mesh object

    Returns
    -------
    BX, BY : ndarray
        div block B = [BX, BY].T in [[A, B], [B.T 0]]
    """
    if mesh.V2 is None:
        mesh.generate_quadratic()

    X, Y = mesh.X, mesh.Y
    ne = mesh.ne
    E = mesh.E2

    m1 = 6
    m2 = 3
    DX = np.zeros((ne, m1*m2))
    DXI = np.zeros((ne, m1*m2), dtype=int)
    DXJ = np.zeros((ne, m1*m2), dtype=int)
    DY = np.zeros((ne, m1*m2))
    DYI = np.zeros((ne, m1*m2), dtype=int)
    DYJ = np.zeros((ne, m1*m2), dtype=int)

    # Assemble A and b
    for ei in range(0, ne):
        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]]

        J = np.array([[x1 - x0, x2 - x0],
                      [y1 - y0, y2 - y0]])
        invJ = np.linalg.inv(J.T)
        detJ = np.linalg.det(J)

        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],   # noqa: E201
                       [-0.81684757298045870,  0.63369514596091750]])
        xx, yy = (xy[:, 0]+1)/2, (xy[:, 1]+1)/2
        ww *= 0.5

        DXelem = np.zeros((3, 6))
        DYelem = np.zeros((3, 6))

        for w, x, y in zip(ww, xx, yy):
            basis1 = np.array([1-x-y, x, y])

            # basis2 = 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]
            ])

            dphi = invJ.dot(dbasis)

            DXelem += (detJ / 2) * w * (np.outer(basis1, dphi[0, :]))
            DYelem += (detJ / 2) * w * (np.outer(basis1, dphi[1, :]))
            dphi.T.dot(dphi)

        # Step 7
        DX[ei, :] = DXelem.ravel()
        DXI[ei, :] = np.repeat(K[np.arange(m2)], m1)
        DXJ[ei, :] = np.tile(K[np.arange(m1)], m2)
        BX = sparse.coo_matrix((DX.ravel(), (DXI.ravel(), DXJ.ravel())))
        BX.sum_duplicates()

        DY[ei, :] = DYelem.ravel()
        DYI[ei, :] = np.repeat(K[np.arange(m2)], m1)
        DYJ[ei, :] = np.tile(K[np.arange(m1)], m2)
        BY = sparse.coo_matrix((DY.ravel(), (DYI.ravel(), DYJ.ravel())))
        BY.sum_duplicates()
    return BX, BY
def divform ( mesh ) : if mesh . V2 is None : mesh . generate_quadratic ( ) X , Y = mesh . X , mesh . Y ne = mesh . ne E = mesh . E2 m1 = 6 m2 = 3 DX = np . zeros ( ( ne , m1 * m2 ) ) DXI = np . zeros ( ( ne , m1 * m2 ) , dtype = int ) DXJ = np . zeros ( ( ne , m1 * m2 ) , dtype = int ) DY = np . zeros ( ( ne , m1 * m2 ) ) DYI = np . zeros ( ( ne , m1 * m2 ) , dtype = int ) DYJ = np . zeros ( ( ne , m1 * m2 ) , dtype = int ) # Assemble A and b for ei in range ( 0 , ne ) : 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 ] ] J = np . array ( [ [ x1 - x0 , x2 - x0 ] , [ y1 - y0 , y2 - y0 ] ] ) invJ = np . linalg . inv ( J . T ) detJ = np . linalg . det ( J ) 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 ] , # noqa: E201 [ - 0.81684757298045870 , 0.63369514596091750 ] ] ) xx , yy = ( xy [ : , 0 ] + 1 ) / 2 , ( xy [ : , 1 ] + 1 ) / 2 ww *= 0.5 DXelem = np . zeros ( ( 3 , 6 ) ) DYelem = np . zeros ( ( 3 , 6 ) ) for w , x , y in zip ( ww , xx , yy ) : basis1 = np . array ( [ 1 - x - y , x , y ] ) # basis2 = 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 ] ] ) dphi = invJ . dot ( dbasis ) DXelem += ( detJ / 2 ) * w * ( np . outer ( basis1 , dphi [ 0 , : ] ) ) DYelem += ( detJ / 2 ) * w * ( np . outer ( basis1 , dphi [ 1 , : ] ) ) dphi . T . dot ( dphi ) # Step 7 DX [ ei , : ] = DXelem . ravel ( ) DXI [ ei , : ] = np . repeat ( K [ np . arange ( m2 ) ] , m1 ) DXJ [ ei , : ] = np . tile ( K [ np . arange ( m1 ) ] , m2 ) BX = sparse . coo_matrix ( ( DX . ravel ( ) , ( DXI . ravel ( ) , DXJ . ravel ( ) ) ) ) BX . sum_duplicates ( ) DY [ ei , : ] = DYelem . ravel ( ) DYI [ ei , : ] = np . repeat ( K [ np . arange ( m2 ) ] , m1 ) DYJ [ ei , : ] = np . tile ( K [ np . arange ( m1 ) ] , m2 ) BY = sparse . coo_matrix ( ( DY . ravel ( ) , ( DYI . ravel ( ) , DYJ . ravel ( ) ) ) ) BY . sum_duplicates ( ) return BX , BY