Skip to content

linear_elasticity_stiffness

linear_elasticity_stiffness(V, F, K=1.75, mu=0.0115, volumes=None, mass=None)

Differential operators needed for linear elasticity calculations

Returns the linear elastic stiffness and strain matrices for a given shape and material parameters

Parameters:

Name Type Description Default
V numpy double array

Matrix of vertex coordinates

 required
F numpy int array

Matrix of triangle indices

 required
K double (optional, default 1.75)

Bulk modulus

 1.75
mu double (optional, default 0.0115)

Material shear modulus

 0.0115
volumes numpy double array (optional, default None)

Vector with the volumes (in 2D) or areas (in 3D) of each mesh element (if None, will be computed)

 None
mass scipy sparse_csr (optional, default None)

The mesh's sparse mass matrix (if None, will be computed)

 None

Returns:

Name Type Description
K scipy sparse.csr_matrix

Stiffness matrix
C scipy sparse.csr_matrix

Constituitive model matrix

strain scipy sparse.csr_matrix

Strain matrix
A scipy csr_matrix

Diagonal element area matrix
M scipy sparse.csr_matrix

Mass matrix (if input mass is not None, this returns the input)

See Also

linear_elasticity.

Notes

This implementation only works for 2D triangle meshes. Tetrahedral meshes will be supported soon.

Examples:

from gpytoolbox import regular_square_mesh, linear_elasticity_stiffness
V, F = regular_square_mesh(3) # Make regular mesh
V = (V + 1.)/2. # Normalize mesh
# Compute linear elasticity operators
K, C, strain, A, M = linear_elasticity_stiffness(V,F)
Source code in src/gpytoolbox/linear_elasticity_stiffness.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def linear_elasticity_stiffness(V,F,K=1.75,mu=0.0115,volumes=None,mass=None):
    """Differential operators needed for linear elasticity calculations

    Returns the linear elastic stiffness and strain matrices for a given shape and material parameters

    Parameters
    ----------
    V : numpy double array
        Matrix of vertex coordinates
    F : numpy int array
        Matrix of triangle indices
    K : double (optional, default 1.75)
        Bulk modulus
    mu : double (optional, default 0.0115)
        Material shear modulus
    volumes : numpy double array (optional, default None)
        Vector with the volumes (in 2D) or areas (in 3D) of each mesh element (if None, will be computed)
    mass : scipy sparse_csr (optional, default None)
        The mesh's sparse mass matrix (if None, will be computed)


    Returns
    -------
    K  : scipy sparse.csr_matrix 
        Stiffness matrix
    C : scipy sparse.csr_matrix 
        Constituitive model matrix 
    strain : scipy sparse.csr_matrix 
        Strain matrix
    A : scipy csr_matrix 
        Diagonal element area matrix
    M : scipy sparse.csr_matrix 
        Mass matrix (if input mass is not None, this returns the input)

    See Also
    --------
    linear_elasticity.

    Notes
    -----
    This implementation only works for 2D triangle meshes. Tetrahedral meshes will be supported soon.

    Examples
    --------
    ```python
    from gpytoolbox import regular_square_mesh, linear_elasticity_stiffness
    V, F = regular_square_mesh(3) # Make regular mesh
    V = (V + 1.)/2. # Normalize mesh
    # Compute linear elasticity operators
    K, C, strain, A, M = linear_elasticity_stiffness(V,F)
    ```
    """

    l = K - (2/3)*mu
    Z = csr_matrix((F.shape[0],V.shape[0]))
    dim = V.shape[1]
    m = F.shape[0]
    I = identity(m)

    if dim==2:
        G = grad(np.concatenate((V, np.zeros((V.shape[0],1))), axis=1),F)
        G = G[0:(2*F.shape[0]),0:(2*F.shape[0])]
        G0 = G[0:F.shape[0],:]
        G1 = G[F.shape[0]:(2*F.shape[0])]
        strain0 = hstack((G0, Z))
        strain1 = hstack((Z, G1))
        strain2 = hstack((G1, G0))
        strain = vstack((strain0,strain1,strain2))
        C = vstack(( hstack(( (l+2*mu)*I,l*I,0*I )),hstack(( l*I,(l+2*mu)*I,0*I )),hstack(( 0*I, 0*I, mu*I )) ) )
        if (volumes is None):
            A = diags(doublearea(V,F)*0.5)
        else:
            A = diags(volumes)
        A = block_diag((A,A,A))
        if (mass is None):
            M = massmatrix(V,F,'voronoi')
        else:
            M = mass

        M = block_diag((M,M))
    if dim==3:
        print("DIMENSION 3 NOT SUPPORTED YET")

    Z = csr_matrix((V.shape[0],F.shape[0]))
    D = strain.transpose()
    K = D * A * C * strain
    K = csr_matrix(K)

    return K, C, strain, A, M