from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from numpy import arange
from numpy import divide
from numpy import hstack
from numpy import tile
from scipy import cross
from scipy.sparse import coo_matrix
from compas.numerical.linalg import normrow
from compas.numerical.linalg import normalizerow
from compas.numerical.linalg import rot90
__all__ = [
'grad',
'div',
'curl'
]
[docs]def grad(V, F, rtype='array'):
"""Construct the gradient operator of a trianglular mesh.
Parameters
----------
V : array
Vertex coordinates of the mesh.
F : array
Face vertex indices of the mesh.
rtype : {'array', 'csc', 'csr', 'coo', 'list'}
Format of the result.
Returns
-------
array-like
Depending on rtype return type.
Notes
-----
The gradient operator is fully determined by the connectivity of the mesh
and the coordinate difference vectors associated with the edges
"""
v = V.shape[0]
f = F.shape[0]
f0 = F[:, 0] # Index of first vertex of each face
f1 = F[:, 1] # Index of second vertex of each face
f2 = F[:, 2] # Index of last vertex of each face
v01 = V[f1, :] - V[f0, :] # Vector from vertex 0 to 1 for each face
v12 = V[f2, :] - V[f1, :] # Vector from vertex 1 to 2 for each face
v20 = V[f0, :] - V[f2, :] # Vector from vertex 2 to 0 for each face
n = cross(v12, v20) # Normal vector to each face
A2 = normrow(n) # Length of normal vector is twice the area of the face
A2 = tile(A2, (1, 3))
u = normalizerow(n) # Unit normals for each face
v01_ = divide(rot90(v01, u), A2) # Vector perpendicular to v01, normalized by A2
v20_ = divide(rot90(v20, u), A2) # Vector perpendicular to v20, normalized by A2
i = hstack(( # Nonzero rows
0 * f + tile(arange(f), (1, 4)),
1 * f + tile(arange(f), (1, 4)),
2 * f + tile(arange(f), (1, 4))
)).flatten()
j = tile(hstack((f1, f0, f2, f0)), (1, 3)).flatten() # Nonzero columns
data = hstack((
hstack((v20_[:, 0], - v20_[:, 0], v01_[:, 0], - v01_[:, 0])),
hstack((v20_[:, 1], - v20_[:, 1], v01_[:, 1], - v01_[:, 1])),
hstack((v20_[:, 2], - v20_[:, 2], v01_[:, 2], - v01_[:, 2])),
)).flatten()
G = coo_matrix((data, (i, j)), shape=(3 * f, v))
if rtype == 'array':
return G.toarray()
elif rtype == 'csr':
return G.tocsr()
elif rtype == 'csc':
return G.tocsc()
elif rtype == 'coo':
return G
else:
return G
def div():
raise NotImplementedError
def curl():
raise NotImplementedError
# ==============================================================================
# Main
# ==============================================================================
if __name__ == "__main__":
import doctest
doctest.testmod(globs=globals())