Moved coboundary closure function.

src/dctkit/dec/cochain.py

@@ -261,46 +261,6 @@ def coboundary(c: Cochain) -> Cochain:
     return dc
 
 
-def coboundary_closure(c: CochainP) -> CochainD:
-    """Implements the operator that complements the coboundary on the boundary
-    of dual (n-1)-simplices, where n is the dimension of the complex.
-
-    Args:
-        c: a primal (n-1)-cochain
-    Returns:
-        the coboundary closure of c, resulting in a dual n-cochain with non-zero
-        coefficients in the "uncompleted" cells.
-    """
-    n = c.complex.dim
-    num_tets = c.complex.S[n].shape[0]
-    num_dual_faces = c.complex.S[n-1].shape[0]
-
-    # to extract only the boundary components with the right orientation
-    # we construct a dual n-2 cochain and we take the (true) coboundary.
-    # In this way the obtain a cochain such that an entry is 0 if it's in
-    # the interior of the complex and ±1 if it's in the boundary
-    ones = CochainD(dim=n-2, complex=c.complex, coeffs=jnp.ones(num_tets))
-    diagonal_elems = coboundary(ones).coeffs
-    diagonal_matrix_rows = jnp.arange(num_dual_faces)
-    diagonal_matrix_cols = diagonal_matrix_rows
-    diagonal_matrix_COO = [diagonal_matrix_rows, diagonal_matrix_cols, diagonal_elems]
-
-    # build the absolute value of the (n-1)-coboundary
-    abs_dual_coboundary_faces = c.complex.boundary[n-1].copy()
-    # same of doing abs(dual_coboundary_faces)
-    abs_dual_coboundary_faces[2] = abs_dual_coboundary_faces[2]**2
-    # with this product, we extract with the right orientation the boundary pieces
-    diagonal_times_c = spmv.spmm(diagonal_matrix_COO, c.coeffs,
-                                 transpose=False,
-                                 shape=c.complex.S[n-1].shape[0])
-    # here we sum their contribution taking into account the orientation
-    d_closure_coeffs = spmv.spmm(abs_dual_coboundary_faces, diagonal_times_c,
-                                 transpose=False,
-                                 shape=c.complex.num_nodes)
-    d_closure = CochainD(dim=n, complex=c.complex, coeffs=0.5*d_closure_coeffs)
-    return d_closure
-
-
 def star(c: Cochain) -> Cochain:
     """Implements the diagonal Hodge star operator (see Grinspun et al.).
 
@@ -388,21 +348,6 @@ def laplacian(c: Cochain) -> Cochain:
     return laplacian
 
 
-def deformation_gradient(c: Cochain) -> Cochain:
-    """Compute the deformation gradient of a primal vector-valued 0-cochain
-    representing the node coordinates.
-
-    Args:
-        c: a primal vector-valued 0 cochain.
-
-    Returns:
-        the deformation gradient for each dual nodes, i.e. a dual tensor-valued
-            0-cochain.
-    """
-    F = c.complex.get_deformation_gradient(c.coeffs)
-    return Cochain(0, not c.is_primal, c.complex, F)
-
-
 def transpose(c: Cochain) -> Cochain:
     """Compute the transpose of a tensor-valued cochain.
 

src/dctkit/dec/flat.py

@@ -1,5 +1,4 @@
 import jax.numpy as jnp
-import dctkit as dt
 from dctkit.dec import cochain as C
 
 
@@ -63,28 +62,3 @@ def flat_DPP(c: C.CochainD0V | C.CochainD0T) -> C.CochainP1:
 
 def flat_PDP(c: C.CochainP0) -> C.CochainP1:
     return C.CochainP1(c.complex, c.complex.flat_PDP_weights @ c.coeffs)
-
-
-def flat_PDD(c: C.CochainD0, scheme: str) -> C.CochainD1:
-    # NOTE: we use periodic boundary conditions
-    # NOTE: only implemented for dim = 1, where dim is the dimension
-    # of the complex
-    dual_volumes = c.complex.dual_volumes[0]
-    # FIXME: rewrite this!
-    if c.coeffs.ndim == 1:
-        flat_c_coeffs = jnp.zeros(c.complex.num_nodes, dtype=dt.float_dtype)
-    elif c.coeffs.ndim == 2:
-        flat_c_coeffs = jnp.zeros(
-            (c.complex.num_nodes, c.coeffs.shape[1]), dtype=dt.float_dtype)
-    if scheme == "upwind":
-        # periodic bc
-        flat_c_coeffs = flat_c_coeffs.at[0].set(dual_volumes[0]*c.coeffs[-1])
-        # upwind implementation
-        flat_c_coeffs = flat_c_coeffs.at[1:].set(dual_volumes[1:]*c.coeffs)
-    elif scheme == "parabolic":
-        # periodic bc
-        flat_c_coeffs = flat_c_coeffs.at[0].set(dual_volumes[0]*c.coeffs[-1])
-        flat_c_coeffs = flat_c_coeffs.at[-1].set(dual_volumes[-1]*c.coeffs[0])
-        flat_c_coeffs = flat_c_coeffs.at[1:-1].set(0.5 * (dual_volumes[1:-1] *
-                                                   (c.coeffs[:-1] + c.coeffs[1:]).T).T)
-    return C.CochainD1(c.complex, flat_c_coeffs)

src/dctkit/experimental/coboundary_closure.py

@@ -0,0 +1,43 @@
+import dctkit.dec.cochain as C
+from dctkit.math import spmv
+import jax.numpy as jnp
+
+
+def coboundary_closure(c: C.CochainP) -> C.CochainD:
+    """Implements the operator that complements the coboundary on the boundary
+    of dual (n-1)-simplices, where n is the dimension of the complex.
+
+    Args:
+        c: a primal (n-1)-cochain
+    Returns:
+        the coboundary closure of c, resulting in a dual n-cochain with non-zero
+        coefficients in the "uncompleted" cells.
+    """
+    n = c.complex.dim
+    num_tets = c.complex.S[n].shape[0]
+    num_dual_faces = c.complex.S[n-1].shape[0]
+
+    # to extract only the boundary components with the right orientation
+    # we construct a dual n-2 cochain and we take the (true) coboundary.
+    # In this way the obtain a cochain such that an entry is 0 if it's in
+    # the interior of the complex and ±1 if it's in the boundary
+    ones = C.CochainD(dim=n-2, complex=c.complex, coeffs=jnp.ones(num_tets))
+    diagonal_elems = C.coboundary(ones).coeffs
+    diagonal_matrix_rows = jnp.arange(num_dual_faces)
+    diagonal_matrix_cols = diagonal_matrix_rows
+    diagonal_matrix_COO = [diagonal_matrix_rows, diagonal_matrix_cols, diagonal_elems]
+
+    # build the absolute value of the (n-1)-coboundary
+    abs_dual_coboundary_faces = c.complex.boundary[n-1].copy()
+    # same of doing abs(dual_coboundary_faces)
+    abs_dual_coboundary_faces[2] = abs_dual_coboundary_faces[2]**2
+    # with this product, we extract with the right orientation the boundary pieces
+    diagonal_times_c = spmv.spmm(diagonal_matrix_COO, c.coeffs,
+                                 transpose=False,
+                                 shape=c.complex.S[n-1].shape[0])
+    # here we sum their contribution taking into account the orientation
+    d_closure_coeffs = spmv.spmm(abs_dual_coboundary_faces, diagonal_times_c,
+                                 transpose=False,
+                                 shape=c.complex.num_nodes)
+    d_closure = C.CochainD(dim=n, complex=c.complex, coeffs=0.5*d_closure_coeffs)
+    return d_closure

src/dctkit/physics/elasticity.py

@@ -1,5 +1,6 @@
 import numpy.typing as npt
 from dctkit.mesh.simplex import SimplicialComplex
+from dctkit.experimental.coboundary_closure import coboundary_closure
 import dctkit.dec.cochain as C
 import dctkit.dec.flat as V
 from jax import Array
@@ -163,7 +164,7 @@ def get_dual_balance(self, node_coords: C.CochainP0,
         # set tractions on given sub-portions of the boundary
         forces_closure_update = self.set_boundary_tractions(
             forces_closure, boundary_tractions)
-        balance_forces_closure = C.coboundary_closure(forces_closure_update)
+        balance_forces_closure = coboundary_closure(forces_closure_update)
         balance = C.add(C.coboundary(forces), balance_forces_closure)
         # balance = C.coboundary(forces)
         return balance

tests/test_cochain.py

@@ -2,6 +2,7 @@
 import dctkit
 from dctkit.mesh import util
 from dctkit.dec import cochain as C
+from dctkit.experimental.coboundary_closure import coboundary_closure
 
 
 def test_coboundary(setup_test):
@@ -607,6 +608,6 @@ def test_coboundary_closure(setup_test):
     S_2.get_hodge_star()
 
     c = C.CochainP1(complex=S_2, coeffs=np.arange(1, 9, dtype=dctkit.float_dtype))
-    cob_clos_c = C.coboundary_closure(c)
+    cob_clos_c = coboundary_closure(c)
     cob_clos_c_true = np.array([-0.5,  2.5,  5.,  2.,  0.], dtype=dctkit.float_dtype)
     assert np.allclose(cob_clos_c.coeffs, cob_clos_c_true)