adding DEC operators

src/gpytoolbox/__init__.py

@@ -117,4 +117,18 @@
 from .adjacency_matrix import adjacency_matrix
 from .non_manifold_edges import non_manifold_edges
 from .connected_components import connected_components
-from .rotation_matrix_from_vectors import rotation_matrix_from_vectors
+from .rotation_matrix_from_vectors import rotation_matrix_from_vectors
+from .dec_d0 import dec_d0
+from .dec_d1 import dec_d1
+from .dec_h0 import dec_h0
+from .dec_h0inv import dec_h0inv
+from .dec_h1 import dec_h1
+from .dec_h1inv import dec_h1inv
+from .dec_h2 import dec_h2
+from .dec_h2inv import dec_h2inv
+from .dec_h0_intrinsic import dec_h0_intrinsic
+from .dec_h0inv_intrinsic import dec_h0inv_intrinsic
+from .dec_h1_intrinsic import dec_h1_intrinsic
+from .dec_h1inv_intrinsic import dec_h1inv_intrinsic
+from .dec_h2_intrinsic import dec_h2_intrinsic
+from .dec_h2inv_intrinsic import dec_h2inv_intrinsic

src/gpytoolbox/cotangent_laplacian_intrinsic.py

@@ -25,7 +25,7 @@ def cotangent_laplacian_intrinsic(l_sq,F,n=None):
     ```python
     # Mesh in V,F
     from gpytoolbox import halfedge_lengths_squared, cotangent_laplacian_intrinsic
-    l = halfedge_lengths_squared(V,F)
+    l_sq = halfedge_lengths_squared(V,F)
     L = cotangent_laplacian_intrinsic(l_sq,F,n=V.shape[0])
     ```
 

src/gpytoolbox/dec_d0.py

@@ -0,0 +1,53 @@
+import numpy as np
+import scipy as sp
+from .halfedge_edge_map import halfedge_edge_map
+
+def dec_d0(F,E=None,n=None):
+    """Builds the DEC d0 operator as described, for example, in Crane et al.
+    2013. "Digital Geometry Processing with Discrete Exterior Calculus".
+
+    The edge labeling in E follows the convention from Gpytoolbox's
+    `halfedge_edge_map`.
+
+    The input mesh _must_ be a manifold mesh.
+
+    Parameters
+    ----------
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    E : (e,2) numpy int array, optional (default None)
+        edge index list of a triangle mesh.
+        If absent, will be computed using `halfedge_edge_map`
+    n : int, optional (default None)
+        number of vertices in the mesh.
+        If absent, will try to infer from F.
+
+    Returns
+    -------
+    d0 : (e,n) scipy csr_matrix
+        DEC operator d0
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    d0 = gpy.dec_d0(F)
+    ```
+
+    """
+
+    assert F.shape[1] == 3
+
+    if n is None:
+        n = np.max(F)+1
+
+    if E is None:
+        _,E,_,_ = halfedge_edge_map(F, assume_manifold=True)
+
+    i = np.concatenate((np.arange(E.shape[0]), np.arange(E.shape[0])), axis=0)
+    j = np.concatenate((E[:,0], E[:,1]), axis=0)
+    k = np.concatenate((np.ones(E.shape[0], dtype=float),
+        -np.ones(E.shape[0], dtype=float)))
+    d0 = sp.sparse.csr_matrix((k, (i,j)), shape=(E.shape[0],n))
+
+    return d0

src/gpytoolbox/dec_d1.py

@@ -0,0 +1,53 @@
+import numpy as np
+import scipy as sp
+from .halfedge_edge_map import halfedge_edge_map
+
+def dec_d1(F,E_to_he=None):
+    """Builds the DEC d1 operator as described, for example, in Crane et al.
+    2013. "Digital Geometry Processing with Discrete Exterior Calculus".
+
+    The edge labeling in E_to_he follows the convention from Gpytoolbox's
+    `halfedge_edge_map`.
+
+    The input mesh _must_ be a manifold mesh.
+
+    Parameters
+    ----------
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    E_to_he : (e,2,2) numpy int array, optional (default None)
+        index map from e to corresponding row and col in the list of
+        all halfedges `he` as computed by `halfedge_edge_map` for two
+        halfedges (or -1 if only one halfedge exists)
+        If absent, will be computed using `halfedge_edge_map`
+
+    Returns
+    -------
+    d1 : (m,e) scipy csr_matrix
+        DEC operator d1
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    d1 = gpy.dec_d1(F)
+    ```
+
+    """
+
+    assert F.shape[1] == 3
+
+    if E_to_he is None:
+        _,_,_,E_to_he = halfedge_edge_map(F, assume_manifold=True)
+
+    # A second halfedge exists for these
+    se = E_to_he[:,1,0] >= 0
+
+    i = np.concatenate((E_to_he[:,0,0], E_to_he[se,1,0]), axis=0)
+    j = np.concatenate((np.arange(E_to_he.shape[0]),
+        np.arange(E_to_he.shape[0])[se]), axis=0)
+    k = np.concatenate((np.ones(E_to_he.shape[0], dtype=float),
+        -np.ones(np.sum(se), dtype=float)), axis=0)
+    d1 = sp.sparse.csr_matrix((k, (i,j)), shape=(F.shape[0],E_to_he.shape[0]))
+
+    return d1

src/gpytoolbox/dec_h0.py

@@ -0,0 +1,33 @@
+import numpy as np
+import scipy as sp
+from .halfedge_lengths_squared import halfedge_lengths_squared
+from .dec_h0_intrinsic import dec_h0_intrinsic
+
+def dec_h0(V,F):
+    """Builds the DEC 0-Hodge-star operator as described, for example, in Crane
+    et al. 2013. "Digital Geometry Processing with Discrete Exterior Calculus".
+
+    Parameters
+    ----------
+    V : (n,d) numpy array
+        vertex list of a triangle mesh
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+
+    Returns
+    -------
+    h0 : (n,n) scipy csr_matrix
+        DEC operator h0
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    h0 = gpy.dec_h0(V,F)
+    ```
+
+    """
+
+    l_sq = halfedge_lengths_squared(V,F)
+    return dec_h0_intrinsic(l_sq,F,n=V.shape[0])
+

src/gpytoolbox/dec_h0_intrinsic.py

@@ -0,0 +1,46 @@
+import numpy as np
+import scipy as sp
+from .doublearea_intrinsic import doublearea_intrinsic
+
+def dec_h0_intrinsic(l_sq,F,n=None):
+    """Builds the DEC 0-Hodge-star operator as described, for example, in Crane
+    et al. 2013. "Digital Geometry Processing with Discrete Exterior Calculus".
+
+    Parameters
+    ----------
+    l_sq : (m,3) numpy array
+        squared halfedge lengths as computed by halfedge_lengths_squared
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    n : int, optional (default None)
+        number of vertices in the mesh.
+        If absent, will try to infer from F.
+
+    Returns
+    -------
+    h0 : (n,n) scipy csr_matrix
+        DEC operator h0
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    l_sq = gpy.halfedge_lengths_squared(V,F)
+    h0 = gpy.dec_h0_intrinsic(l_sq,F)
+    ```
+
+    """
+
+    assert F.shape[1] == 3
+
+    if n is None:
+        n = np.max(F)+1
+
+    A3 = 0.5/3. * doublearea_intrinsic(l_sq,F)
+    i = np.concatenate((F[:,0],F[:,1],F[:,2]), axis=0)
+    j = np.concatenate((F[:,0],F[:,1],F[:,2]), axis=0)
+    k = np.concatenate((A3,A3,A3), axis=0)
+    h0 = sp.sparse.csr_matrix((k, (i,j)), shape=(n,n))
+
+    return h0
+

src/gpytoolbox/dec_h0inv.py

@@ -0,0 +1,33 @@
+import numpy as np
+import scipy as sp
+from .halfedge_lengths_squared import halfedge_lengths_squared
+from .dec_h0inv_intrinsic import dec_h0inv_intrinsic
+
+def dec_h0inv(V,F):
+    """Builds the inverse DEC 0-Hodge-star operator as described, for example,
+    in Crane et al. 2013. "Digital Geometry Processing with Discrete Exterior
+    Calculus".
+
+    Parameters
+    ----------
+    V : (n,d) numpy array
+        vertex list of a triangle mesh
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+
+    Returns
+    -------
+    h0inv : (n,n) scipy csr_matrix
+        inverse of DEC operator h0
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    h0inv = gpy.dec_h0inv(V,F)
+    ```
+
+    """
+
+    l_sq = halfedge_lengths_squared(V,F)
+    return dec_h0inv_intrinsic(l_sq,F,n=V.shape[0])

src/gpytoolbox/dec_h0inv_intrinsic.py

@@ -0,0 +1,46 @@
+import numpy as np
+import scipy as sp
+from .doublearea_intrinsic import doublearea_intrinsic
+
+def dec_h0inv_intrinsic(l_sq,F,n=None):
+    """Builds the inverse DEC 0-Hodge-star operator as described, for example,
+    in Crane et al. 2013. "Digital Geometry Processing with Discrete Exterior
+    Calculus".
+
+    Parameters
+    ----------
+    l_sq : (m,3) numpy array
+        squared halfedge lengths as computed by halfedge_lengths_squared
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    n : int, optional (default None)
+        number of vertices in the mesh.
+        If absent, will try to infer from F.
+
+    Returns
+    -------
+    h0inv : (n,n) scipy csr_matrix
+        inverse of DEC operator h0
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    l_sq = gpy.halfedge_lengths_squared(V,F)
+    h0inv = gpy.dec_h0inv_intrinsic(l_sq,F)
+    ```
+
+    """
+
+    assert F.shape[1] == 3
+
+    if n is None:
+        n = np.max(F)+1
+
+    A3 = 0.5/3. * doublearea_intrinsic(l_sq,F)
+    i = np.concatenate((F[:,0],F[:,1],F[:,2]), axis=0)
+    j = np.concatenate((F[:,0],F[:,1],F[:,2]), axis=0)
+    k = np.nan_to_num(1. / np.concatenate((A3,A3,A3), axis=0))
+    h0inv = sp.sparse.csr_matrix((k, (i,j)), shape=(n,n))
+
+    return h0inv

src/gpytoolbox/dec_h1.py

@@ -0,0 +1,42 @@
+import numpy as np
+import scipy as sp
+from .halfedge_lengths_squared import halfedge_lengths_squared
+from .dec_h1_intrinsic import dec_h1_intrinsic
+
+def dec_h1(V,F,E_to_he=None):
+    """Builds the DEC 1-Hodge-star operator as described, for example, in Crane
+    et al. 2013. "Digital Geometry Processing with Discrete Exterior Calculus".
+
+    The edge labeling in E_to_he follows the convention from Gpytoolbox's
+    `halfedge_edge_map`.
+
+    The input mesh _must_ be a manifold mesh.
+
+    Parameters
+    ----------
+    V : (n,d) numpy array
+        vertex list of a triangle mesh
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    E_to_he : (e,2,2) numpy int array, optional (default None)
+        index map from e to corresponding row and col in the list of
+        all halfedges `he` as computed by `halfedge_edge_map` for two
+        halfedges (or -1 if only one halfedge exists)
+        If absent, will be computed using `halfedge_edge_map`
+
+    Returns
+    -------
+    h1 : (e,e) scipy csr_matrix
+        DEC operator h1
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    h1 = gpy.dec_h1(V,F)
+    ```
+
+    """
+
+    l_sq = halfedge_lengths_squared(V,F)
+    return dec_h1_intrinsic(l_sq,F,E_to_he=E_to_he)