Skip to content

massmatrix_intrinsic

massmatrix_intrinsic(l_sq, F, n=None, type='voronoi')

FEM intrinsic mass matrix

Builds the finite elements mass matrix for a triangle mesh using a piecewise linear hat function basis, using only intrinsic information (squared halfedge edge lengths).

Parameters:

Name Type Description Default
l_sq (m,3) numpy double array

Vector of squared halfedge lengths as computed by halfedge_lengths_squared

 required
F (m,3) numpy int array

face index list of a triangle mesh (into a V assumed to exist)

 required
n int, optional (default: None)

Integer denoting the number of vertices in the mesh

 None
type str, optional (default: 'voronoi')

Type of mass matrix computation: 'voronoi' (default), 'full' or 'barycentric'

 'voronoi'

Returns:

Name Type Description
M (n,n) scipy sparse.csr_matrix

Intrinsicly computed mass matrix
See Also

massmatrix.

Notes

This implementation is lifted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/massmatrix_intrinsic.m

Examples:

from gpytoolbox import massmatrix_intrinsic
l_sq = np.array([[1.,1.,1.]]) # edge lengths
f = np.array([[0,1,2]],dtype=int) # simple triangle
M_b = massmatrix_intrinsic(l_sq,f, type='barycentric')
# M_b will be a diagonal matrix with 0.25/sqrt(3) on the diagonal
Source code in src/gpytoolbox/massmatrix_intrinsic.py 
  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
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def massmatrix_intrinsic(l_sq,F,n=None,type='voronoi'):
    """FEM intrinsic mass matrix

    Builds the finite elements mass matrix for a triangle mesh using a piecewise
    linear hat function basis, using only intrinsic information (squared
    halfedge edge lengths).

    Parameters
    ----------
    l_sq : (m,3) numpy double array
        Vector of squared halfedge lengths as computed by halfedge_lengths_squared
    F : (m,3) numpy int array
        face index list of a triangle mesh (into a V assumed to exist)
    n : int, optional (default: None)
        Integer denoting the number of vertices in the mesh
    type : str, optional (default: 'voronoi')
        Type of mass matrix computation: 'voronoi' (default), 'full' or 'barycentric'

    Returns
    -------
    M : (n,n) scipy sparse.csr_matrix
        Intrinsicly computed mass matrix

    See Also
    --------
    massmatrix.

    Notes
    -----
    This implementation is lifted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/massmatrix_intrinsic.m

    Examples
    --------
    ```python
    from gpytoolbox import massmatrix_intrinsic
    l_sq = np.array([[1.,1.,1.]]) # edge lengths
    f = np.array([[0,1,2]],dtype=int) # simple triangle
    M_b = massmatrix_intrinsic(l_sq,f, type='barycentric')
    # M_b will be a diagonal matrix with 0.25/sqrt(3) on the diagonal
    ```
    """

    assert F.shape == l_sq.shape
    assert F.shape[1]==3
    assert np.all(l_sq >= 0)

    dictionary ={
    'voronoi' : 0,
    'barycentric' : 1,
    'full' : 2
    }
    massmatrix_type = dictionary.get(type,-1)

    if n==None:
        n = np.max(F)+1


    dblA = doublearea_intrinsic(l_sq,F)
    if massmatrix_type==0:
        #Voronoi

        l = np.sqrt(l_sq)
        cosines = np.stack(( 
        ((l_sq[:,2]+l_sq[:,1]-l_sq[:,0])/(2.*l[:,2]*l[:,1])),
        ((l_sq[:,0]+l_sq[:,2]-l_sq[:,1])/(2.*l[:,0]*l[:,2])),
        ((l_sq[:,1]+l_sq[:,0]-l_sq[:,2])/(2.*l[:,1]*l[:,0]))
        ), axis=-1)

        # cosines = [ ...
        # (l(:,3).^2+l(:,2).^2-l(:,1).^2)./(2*l(:,2).*l(:,3)), ...
        #     (l(:,1).^2+l(:,3).^2-l(:,2).^2)./(2*l(:,1).*l(:,3)), ...
        #     (l(:,1).^2+l(:,2).^2-l(:,3).^2)./(2*l(:,1).*l(:,2))];

        barycentric = cosines*l
        normalized_barycentric = barycentric/np.hstack(( np.sum(barycentric,axis=1)[:,None], np.sum(barycentric,axis=1)[:,None], np.sum(barycentric,axis=1)[:,None] ))

        # barycentric = cosines.*l;
        # normalized_barycentric = barycentric./ ...
        #     [sum(barycentric')' sum(barycentric')' sum(barycentric')'];

        partial_triangle_areas = normalized_barycentric * 0.5 * np.stack((dblA,dblA,dblA), axis=-1)

        # partial_triangle_areas = normalized_barycentric.*[areas areas areas];

        quads = np.stack(( 
        ((partial_triangle_areas[:,1]+ partial_triangle_areas[:,2])*0.5),
        ((partial_triangle_areas[:,0]+ partial_triangle_areas[:,2])*0.5),
        ((partial_triangle_areas[:,0]+ partial_triangle_areas[:,1])*0.5)
        ), axis=-1)

        # quads = [ (partial_triangle_areas(:,2)+ partial_triangle_areas(:,3))*0.5 ...
        #     (partial_triangle_areas(:,1)+ partial_triangle_areas(:,3))*0.5 ...
        #     (partial_triangle_areas(:,1)+ partial_triangle_areas(:,2))*0.5];

        c0s = cosines[:,0]<0
        quads[c0s,:] = np.stack((
            0.25*dblA[c0s], 0.125*dblA[c0s], 0.125*dblA[c0s]
            ), axis=-1)
        c1s = cosines[:,1]<0
        quads[c1s,:] = np.stack((
            0.125*dblA[c1s], 0.25*dblA[c1s], 0.125*dblA[c1s]
            ), axis=-1)
        c2s = cosines[:,2]<0
        quads[c2s,:] = np.stack((
            0.125*dblA[c2s], 0.125*dblA[c2s], 0.25*dblA[c2s]
            ), axis=-1)

        # quads(cosines(:,1)<0,:) = [areas(cosines(:,1)<0,:)*0.5, ...
        #     areas(cosines(:,1)<0,:)*0.25, areas(cosines(:,1)<0,:)*0.25];
        # quads(cosines(:,2)<0,:) = [areas(cosines(:,2)<0,:)*0.25, ...
        #     areas(cosines(:,2)<0,:)*0.5, areas(cosines(:,2)<0,:)*0.25];
        # quads(cosines(:,3)<0,:) = [areas(cosines(:,3)<0,:)*0.25, ...
        #     areas(cosines(:,3)<0,:)*0.25, areas(cosines(:,3)<0,:)*0.5];

        I = np.concatenate((F[:,0],F[:,1],F[:,2]))
        J = I
        vals = np.reshape(quads,(-1,1),order='F').squeeze()

        # i = [i1 i2 i3];
        # j = [i1 i2 i3];
        # v = reshape(quads,size(quads,1)*3,1);

    elif massmatrix_type==1:
        #Barycentric

        I = np.concatenate((F[:,0],F[:,1],F[:,2]))
        J = I
        vals = np.concatenate((dblA,dblA,dblA))/6.

    elif massmatrix_type==2:
        #Full
        I = np.concatenate((F[:,0], F[:,1], F[:,1], F[:,2], F[:,2], F[:,0],
            F[:,0], F[:,1], F[:,2]))
        J = np.concatenate((F[:,1], F[:,0], F[:,2], F[:,1], F[:,0], F[:,2],
            F[:,0], F[:,1], F[:,2]))
        offd = dblA / 24.
        diag = dblA / 12.
        vals = np.concatenate((offd, offd, offd, offd, offd, offd,
            diag, diag, diag))

        # i = [i1 i2 i2 i3 i3 i1  i1 i2 i3];
        # j = [i2 i1 i3 i2 i1 i3  i1 i2 i3];
        # offd_v = dblA/24.;
        # diag_v = dblA/12.;
        # v = [offd_v,offd_v, offd_v,offd_v, offd_v,offd_v, diag_v,diag_v,diag_v];  

    else:
        assert False, "invalid massmatrix type"

    M = csr_matrix((vals,(I,J)),shape=(n,n))

    return M