(vertices, elements, E=1e5, nu=0.3, format=None)
return A.asformat(format), B
P1 elements in 2 or 3 dimensions.

    Parameters
    ----------
    vertices : array_like
        array of vertices of a triangle or tets
    elements : array_like
        array of vertex indices for tri or tet elements
    E : float
        Young's modulus
    nu : float
        Poisson's ratio
    format : string
        'csr', 'csc', 'coo', 'bsr'

    Returns
    -------
    A : csr_matrix
        FE Q1 stiffness matrix

    Notes
    -----
        - works in both 2d and in 3d

    Examples
    --------
    >>> from pyamg.gallery import linear_elasticity_p1
    >>> import numpy as np
    >>> E = np.array([[0, 1, 2],[1, 3, 2]])
    >>> V = np.array([[0.0, 0.0],[1.0, 0.0],[0.0, 1.0],[1.0, 1.0]])
    >>> A, B = linear_elasticity_p1(V, E)

    References
    ----------
    .. [1] J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI
       "Matlab implementation of the finite element method in elasticity"
       Computing, Volume 69,  Issue 3  (November 2002) Pages: 239 - 263
       http://www.math.hu-berlin.de/~cc/
P1 elements in 2 or 3 dimensions.
P1 elements in 2 or 3 dimensions .
def linear_elasticity_p1(vertices, elements, E=1e5, nu=0.3, format=None):
    """P1 elements in 2 or 3 dimensions.

    Parameters
    ----------
    vertices : array_like
        array of vertices of a triangle or tets
    elements : array_like
        array of vertex indices for tri or tet elements
    E : float
        Young's modulus
    nu : float
        Poisson's ratio
    format : string
        'csr', 'csc', 'coo', 'bsr'

    Returns
    -------
    A : csr_matrix
        FE Q1 stiffness matrix

    Notes
    -----
        - works in both 2d and in 3d

    Examples
    --------
    >>> from pyamg.gallery import linear_elasticity_p1
    >>> import numpy as np
    >>> E = np.array([[0, 1, 2],[1, 3, 2]])
    >>> V = np.array([[0.0, 0.0],[1.0, 0.0],[0.0, 1.0],[1.0, 1.0]])
    >>> A, B = linear_elasticity_p1(V, E)

    References
    ----------
    .. [1] J. Alberty, C. Carstensen, S. A. Funken, and R. KloseDOI
       "Matlab implementation of the finite element method in elasticity"
       Computing, Volume 69,  Issue 3  (November 2002) Pages: 239 - 263
       http://www.math.hu-berlin.de/~cc/

    """
    # compute local stiffness matrix
    lame = E * nu / ((1 + nu) * (1 - 2*nu))  # Lame's first parameter
    mu = E / (2 + 2*nu)                   # shear modulus

    vertices = np.asarray(vertices)
    elements = np.asarray(elements)

    D = vertices.shape[1]    # spatial dimension
    DoF = D*vertices.shape[0]  # number of degrees of freedom
    NE = elements.shape[0]    # number of elements

    if elements.shape[1] != D + 1:
        raise ValueError('dimension mismatch')

    if D not in (2, 3):
        raise ValueError('only dimension 2 and 3 are supported')

    if D == 2:
        local_K = p12d_local
    elif D == 3:
        local_K = p13d_local

    row = elements.repeat(D).reshape(-1, D)
    row *= D
    row += np.arange(D)
    row = row.reshape(-1, D*(D+1)).repeat(D*(D+1), axis=0)
    row = row.reshape(-1, D*(D+1), D*(D+1))
    col = row.swapaxes(1, 2)

    data = np.empty((NE, D*(D+1), D*(D+1)), dtype=float)

    for i in range(NE):
        element_indices = elements[i, :]
        element_vertices = vertices[element_indices, :]

        data[i] = local_K(element_vertices, lame, mu)

    row = row.ravel()
    col = col.ravel()
    data = data.ravel()

    # sum duplicates
    A = sparse.coo_matrix((data, (row, col)), shape=(DoF, DoF)).tocsr()
    A = A.tobsr(blocksize=(D, D))

    # compute rigid body modes
    if D == 2:
        B = np.zeros((DoF, 3))
        B[0::2, 0] = 1              # vector field in x direction
        B[1::2, 1] = 1              # vector field in y direction

        B[0::2, 2] = -vertices[:, 1]  # rotation vector field (-y, x)
        B[1::2, 2] = vertices[:, 0]
    else:
        B = np.zeros((DoF, 6))
        B[0::3, 0] = 1              # vector field in x direction
        B[1::3, 1] = 1              # vector field in y direction
        B[2::3, 2] = 1              # vector field in z direction

        B[0::3, 3] = -vertices[:, 1]  # rotation vector field (-y, x, 0)
        B[1::3, 3] = vertices[:, 0]
        B[0::3, 4] = -vertices[:, 2]  # rotation vector field (-z, 0, x)
        B[2::3, 4] = vertices[:, 0]
        B[1::3, 5] = -vertices[:, 2]  # rotation vector field (0,-z, y)
        B[2::3, 5] = vertices[:, 1]

    return A.asformat(format), B
def linear_elasticity_p1 ( vertices , elements , E = 1e5 , nu = 0.3 , format = None ) : # compute local stiffness matrix lame = E * nu / ( ( 1 + nu ) * ( 1 - 2 * nu ) ) # Lame's first parameter mu = E / ( 2 + 2 * nu ) # shear modulus vertices = np . asarray ( vertices ) elements = np . asarray ( elements ) D = vertices . shape [ 1 ] # spatial dimension DoF = D * vertices . shape [ 0 ] # number of degrees of freedom NE = elements . shape [ 0 ] # number of elements if elements . shape [ 1 ] != D + 1 : raise ValueError ( 'dimension mismatch' ) if D not in ( 2 , 3 ) : raise ValueError ( 'only dimension 2 and 3 are supported' ) if D == 2 : local_K = p12d_local elif D == 3 : local_K = p13d_local row = elements . repeat ( D ) . reshape ( - 1 , D ) row *= D row += np . arange ( D ) row = row . reshape ( - 1 , D * ( D + 1 ) ) . repeat ( D * ( D + 1 ) , axis = 0 ) row = row . reshape ( - 1 , D * ( D + 1 ) , D * ( D + 1 ) ) col = row . swapaxes ( 1 , 2 ) data = np . empty ( ( NE , D * ( D + 1 ) , D * ( D + 1 ) ) , dtype = float ) for i in range ( NE ) : element_indices = elements [ i , : ] element_vertices = vertices [ element_indices , : ] data [ i ] = local_K ( element_vertices , lame , mu ) row = row . ravel ( ) col = col . ravel ( ) data = data . ravel ( ) # sum duplicates A = sparse . coo_matrix ( ( data , ( row , col ) ) , shape = ( DoF , DoF ) ) . tocsr ( ) A = A . tobsr ( blocksize = ( D , D ) ) # compute rigid body modes if D == 2 : B = np . zeros ( ( DoF , 3 ) ) B [ 0 : : 2 , 0 ] = 1 # vector field in x direction B [ 1 : : 2 , 1 ] = 1 # vector field in y direction B [ 0 : : 2 , 2 ] = - vertices [ : , 1 ] # rotation vector field (-y, x) B [ 1 : : 2 , 2 ] = vertices [ : , 0 ] else : B = np . zeros ( ( DoF , 6 ) ) B [ 0 : : 3 , 0 ] = 1 # vector field in x direction B [ 1 : : 3 , 1 ] = 1 # vector field in y direction B [ 2 : : 3 , 2 ] = 1 # vector field in z direction B [ 0 : : 3 , 3 ] = - vertices [ : , 1 ] # rotation vector field (-y, x, 0) B [ 1 : : 3 , 3 ] = vertices [ : , 0 ] B [ 0 : : 3 , 4 ] = - vertices [ : , 2 ] # rotation vector field (-z, 0, x) B [ 2 : : 3 , 4 ] = vertices [ : , 0 ] B [ 1 : : 3 , 5 ] = - vertices [ : , 2 ] # rotation vector field (0,-z, y) B [ 2 : : 3 , 5 ] = vertices [ : , 1 ] return A . asformat ( format ) , B