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 | def grad(V,F=None):
"""Finite element gradient matrix
Given a triangle mesh or a polyline, computes the finite element gradient matrix assuming piecewise linear hat function basis.
Parameters
----------
V : (n,d) numpy array
vertex list of a polyline or triangle mesh
F : numpy int array, optional (default: None)
if None, interpret input as ordered polyline;
if (m,3) numpy int array, interpred as face index list of a triangle
mesh
Returns
-------
G : (d*m,n) scipy sparse.csr_matrix
Sparse FEM gradient matrix
See Also
--------
cotangent_laplacian.
Notes
-----
Adapted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad.m
Examples
--------
```python
# Mesh in V,F
from gpytoolbox import grad
G = grad(V,F)
```
"""
if (F is None):
F = edge_indices(V.shape[0])
dim = V.shape[1]
simplex_size = F.shape[1]
# Option 1: simplex size is two
if simplex_size==2:
# Then this is just finite difference with varying edge lengths
edge_lengths = np.linalg.norm(V[F[:,1],:] - V[F[:,0],:],axis=1)
I = np.linspace(0,F.shape[0]-1,F.shape[0],dtype=int)
I = np.concatenate((I,I))
J = np.concatenate((F[:,0],F[:,1]))
vals = np.ones(F.shape[0])/edge_lengths
vals = np.concatenate((-vals,vals))
G = csr_matrix((vals,(I,J)),shape=(F.shape[0],V.shape[0]))
if simplex_size==3:
# There are two options: dimension two or three. If it's two, we'll add a third zero dimension for convenience
if dim==2:
V = np.hstack((V,np.zeros((V.shape[0],1))))
# Gradient of a scalar function defined on piecewise linear elements (mesh)
# is constant on each triangle i,j,k:
#
# renaming indices of vertices of triangles for convenience
i0 = F[:,0]
i1 = F[:,1]
i2 = F[:,2]
# F x 3 matrices of triangle edge vectors, named after opposite vertices
v21 = V[i2,:] - V[i1,:]
v02 = V[i0,:] - V[i2,:]
v10 = V[i1,:] - V[i0,:]
# area of parallelogram is twice area of triangle
n = np.cross(v21,v02,axis=1)
# This does correct l2 norm of rows, so that it contains #F list of twice
# triangle areas
dblA = np.linalg.norm(n,axis=1)
u = n/np.tile(dblA[:,None],(1,3))
eperp10 = np.cross(u,v10,axis=1)*np.tile(np.linalg.norm(v10,axis=1)[:,None]/(dblA[:,None]*np.linalg.norm(np.cross(u,v10,axis=1),axis=1)[:,None]),(1,3))
eperp02 = np.cross(u,v02,axis=1)*np.tile(np.linalg.norm(v02,axis=1)[:,None]/(dblA[:,None]*np.linalg.norm(np.cross(u,v02,axis=1),axis=1)[:,None]) ,(1,3))
Find = np.linspace(0,F.shape[0]-1,F.shape[0],dtype=int)
I = np.concatenate((Find,Find,Find,Find,
F.shape[0]+Find,F.shape[0]+Find,F.shape[0]+Find,F.shape[0]+Find,
2*F.shape[0]+Find,2*F.shape[0]+Find,2*F.shape[0]+Find,2*F.shape[0]+Find))
J = np.concatenate((F[:,1],F[:,0],F[:,2],F[:,0]))
J = np.concatenate((J,J,J))
vals = np.concatenate((eperp02[:,0],-eperp02[:,0],eperp10[:,0],-eperp10[:,0],
eperp02[:,1],-eperp02[:,1],eperp10[:,1],-eperp10[:,1],
eperp02[:,2],-eperp02[:,2],eperp10[:,2],-eperp10[:,2]))
G = csr_matrix((vals,(I,J)),shape=(3*F.shape[0],V.shape[0]))
if dim == 2:
G = G[0:(2*F.shape[0]),:]
return G
|