Source code for compas.numerical.operators


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())