456 improve semantics of functions to clarify what simplex they operate on

src/wmtk/attribute/AttributeHandle.hpp

@@ -74,6 +74,7 @@ class MeshAttributeHandle
                m_primitive_type == o.m_primitive_type;
     }
     bool is_valid() const { return m_base_handle.is_valid(); }
+    PrimitiveType primitive_type() const { return m_primitive_type; }
 };
 } // namespace attribute
 using AttributeHandle = attribute::AttributeHandle;

src/wmtk/function/AMIPS.cpp

@@ -3,7 +3,7 @@
 
 namespace wmtk::function {
 AMIPS::AMIPS(const TriMesh& mesh, const MeshAttributeHandle<double>& vertex_attribute_handle)
-    : AutodiffFunction(mesh, vertex_attribute_handle)
+    : AutodiffFunction(mesh, PrimitiveType::Face, vertex_attribute_handle)
 {}
 
 } // namespace wmtk::function

src/wmtk/function/AMIPS2D.cpp

@@ -8,23 +8,24 @@ namespace wmtk::function {
 AMIPS2D::AMIPS2D(const TriMesh& mesh, const MeshAttributeHandle<double>& vertex_attribute_handle)
     : AMIPS(mesh, vertex_attribute_handle)
 {
-    assert(get_vertex_attribute_handle().is_valid());
+    assert(get_variable_attribute_handle().is_valid());
     // check the dimension of the position
     assert(embedded_dimension() == 2);
 }
 
 
-auto AMIPS2D::get_value_autodiff(const Tuple& tuple) const -> DScalar
+auto AMIPS2D::get_value_autodiff(const Simplex& simplex) const -> DScalar
 {
-    return function_eval<DScalar>(tuple);
+    assert(simplex.primitive_type() == PrimitiveType::Vertex);
+    return function_eval<DScalar>(simplex.tuple());
 }
 
 template <typename T>
 T AMIPS2D::function_eval(const Tuple& tuple) const
 {
     // get_autodiff_value sets the autodiff size if necessary
     // get the uv coordinates of the triangle
-    ConstAccessor<double> pos = mesh().create_const_accessor(get_vertex_attribute_handle());
+    ConstAccessor<double> pos = mesh().create_const_accessor(get_variable_attribute_handle());
 
     auto tuple_value = pos.const_vector_attribute(tuple);
     Vector2<T> uv0;

src/wmtk/function/AMIPS2D.hpp

@@ -7,7 +7,7 @@ class AMIPS2D : public AMIPS
     AMIPS2D(const TriMesh& mesh, const MeshAttributeHandle<double>& vertex_attribute_handle);
 
 protected:
-    DScalar get_value_autodiff(const Tuple& tuple) const override;
+    DScalar get_value_autodiff(const Simplex& simplex) const override;
 
 
     template <typename T>

src/wmtk/function/AMIPS3D.cpp

@@ -8,23 +8,24 @@ namespace wmtk::function {
 AMIPS3D::AMIPS3D(const TriMesh& mesh, const MeshAttributeHandle<double>& vertex_attribute_handle)
     : AMIPS(mesh, vertex_attribute_handle)
 {
-    assert(get_vertex_attribute_handle().is_valid());
+    assert(get_variable_attribute_handle().is_valid());
     // check the dimension of the position
     assert(embedded_dimension() == 3);
 }
 
 
-auto AMIPS3D::get_value_autodiff(const Tuple& tuple) const -> DScalar
+auto AMIPS3D::get_value_autodiff(const Simplex& simplex) const -> DScalar
 {
-    return function_eval<DScalar>(tuple);
+    assert(simplex.primitive_type() == PrimitiveType::Vertex);
+    return function_eval<DScalar>(simplex.tuple());
 }
 
 template <typename T>
 T AMIPS3D::function_eval(const Tuple& tuple) const
 {
     // get_autodiff_value sets the autodiff size if necessary
     // get the pos coordinates of the triangle
-    ConstAccessor<double> pos = mesh().create_const_accessor(get_vertex_attribute_handle());
+    ConstAccessor<double> pos = mesh().create_const_accessor(get_variable_attribute_handle());
 
     auto tuple_value = pos.const_vector_attribute(tuple);
     Vector3<T> pos0;

src/wmtk/function/AMIPS3D.hpp

@@ -8,7 +8,7 @@ class AMIPS3D : public AMIPS
     AMIPS3D(const TriMesh& mesh, const MeshAttributeHandle<double>& vertex_attribute_handle);
 
 protected:
-    DScalar get_value_autodiff(const Tuple& tuple) const override;
+    DScalar get_value_autodiff(const Simplex& simplex) const override;
 
     template <typename T>
     T function_eval(const Tuple& tuple) const;

src/wmtk/function/AutodiffFunction.cpp

@@ -5,25 +5,26 @@ namespace wmtk::function {
 
 AutodiffFunction::AutodiffFunction(
     const Mesh& mesh,
-    const MeshAttributeHandle<double>& vertex_attribute_handle)
-    : DifferentiableFunction(mesh, vertex_attribute_handle)
+    const PrimitiveType& simplex_type,
+    const MeshAttributeHandle<double>& variable_attribute_handle)
+    : DifferentiablePerSimplexFunction(mesh, simplex_type, variable_attribute_handle)
 {}
 
 AutodiffFunction::~AutodiffFunction() = default;
-double AutodiffFunction::get_value(const Tuple& tuple) const
+double AutodiffFunction::get_value(const Simplex& simplex) const
 {
     auto scope = utils::AutoDiffRAII(embedded_dimension());
-    return get_value_autodiff(tuple).getValue();
+    return get_value_autodiff(simplex).getValue();
 }
-Eigen::VectorXd AutodiffFunction::get_gradient(const Tuple& tuple) const
+Eigen::VectorXd AutodiffFunction::get_gradient(const Simplex& simplex) const
 {
     auto scope = utils::AutoDiffRAII(embedded_dimension());
-    auto v = get_value_autodiff(tuple);
+    auto v = get_value_autodiff(simplex);
     return v.getGradient();
 }
-Eigen::MatrixXd AutodiffFunction::get_hessian(const Tuple& tuple) const
+Eigen::MatrixXd AutodiffFunction::get_hessian(const Simplex& simplex) const
 {
     auto scope = utils::AutoDiffRAII(embedded_dimension());
-    return get_value_autodiff(tuple).getHessian();
+    return get_value_autodiff(simplex).getHessian();
 }
 } // namespace wmtk::function

src/wmtk/function/AutodiffFunction.hpp

@@ -1,25 +1,28 @@
 #pragma once
 #include <wmtk/function/utils/autodiff.h>
-#include "DifferentiableFunction.hpp"
+#include "DifferentiablePerSimplexFunction.hpp"
 namespace wmtk::function {
 
-class AutodiffFunction : public DifferentiableFunction
+class AutodiffFunction : public DifferentiablePerSimplexFunction
 {
 public:
     using DScalar = DScalar2<double, Eigen::Matrix<double, -1, 1>, Eigen::Matrix<double, -1, -1>>;
     using Scalar = typename DScalar::Scalar;
     static_assert(
         std::is_same_v<Scalar, double>); // MTAO: i'm leaving scalar here but is it ever not double?
-    AutodiffFunction(const Mesh& mesh, const MeshAttributeHandle<double>& vertex_attribute_handle);
+    AutodiffFunction(
+        const Mesh& mesh,
+        const PrimitiveType& simplex_type,
+        const attribute::MeshAttributeHandle<double>& variable_attribute_handle);
 
     virtual ~AutodiffFunction();
 
 public:
-    double get_value(const Tuple& tuple) const override;
-    Eigen::VectorXd get_gradient(const Tuple& tuple) const override;
-    Eigen::MatrixXd get_hessian(const Tuple& tuple) const override;
+    double get_value(const Simplex& simplex) const override;
+    Eigen::VectorXd get_gradient(const Simplex& simplex) const override;
+    Eigen::MatrixXd get_hessian(const Simplex& simplex) const override;
 
 protected:
-    virtual DScalar get_value_autodiff(const Tuple& tuple) const = 0;
+    virtual DScalar get_value_autodiff(const Simplex& simplex) const = 0;
 };
 } // namespace wmtk::function

src/wmtk/function/CMakeLists.txt

@@ -2,8 +2,14 @@
 set(SRC_FILES
     Function.hpp    
     Function.cpp
-    DifferentiableFunction.cpp
-    DifferentiableFunction.hpp
+
+    LocallyDifferentiableFunction.hpp
+    LocallyDifferentiableFunction.cpp
+    PerSimplexFunction.hpp
+    PerSimplexFunction.cpp
+    DifferentiablePerSimplexFunction.hpp
+    DifferentiablePerSimplexFunction.cpp
+
     AMIPS.hpp
     AMIPS.cpp
     AMIPS2D.hpp
@@ -12,9 +18,11 @@ set(SRC_FILES
     AMIPS3D.cpp
     PositionMapAMIPS2D.hpp
     PositionMapAMIPS2D.cpp
-
+    ValenceEnergyPerEdge.hpp
+    ValenceEnergyPerEdge.cpp
     TriMeshValenceFunction.hpp
     TriMeshValenceFunction.cpp
+
     AutodiffFunction.hpp
     AutodiffFunction.cpp
 )

src/wmtk/function/DifferentiablePerSimplexFunction.cpp

@@ -0,0 +1,40 @@
+#include "DifferentiablePerSimplexFunction.hpp"
+
+namespace wmtk {
+namespace function {
+DifferentiablePerSimplexFunction::DifferentiablePerSimplexFunction(
+    const Mesh& mesh,
+    const PrimitiveType& simplex_type,
+    const attribute::MeshAttributeHandle<double>& variable_attribute_handle)
+    : PerSimplexFunction(mesh, simplex_type)
+    , m_attribute_handle(variable_attribute_handle)
+{}
+
+DifferentiablePerSimplexFunction::~DifferentiablePerSimplexFunction() = default;
+
+const MeshAttributeHandle<double>& DifferentiablePerSimplexFunction::get_variable_attribute_handle()
+    const
+{
+    return m_attribute_handle;
+}
+
+long DifferentiablePerSimplexFunction::embedded_dimension() const
+{
+    return mesh().get_attribute_dimension(get_variable_attribute_handle());
+}
+
+
+void DifferentiablePerSimplexFunction::assert_function_type(const PrimitiveType& s_type) const
+{
+    if (get_variable_attribute_handle().primitive_type() != s_type) {
+        throw std::runtime_error("Differentiation of the DifferentiableFunction must be taken wrt "
+                                 "the attribute of the simplex type ");
+    }
+    if (get_simplex_type() < s_type) {
+        throw std::runtime_error(
+            "The DifferentiableFunction must be defined on the cofaces of higher dimension to the "
+            "simplex type ");
+    }
+}
+} // namespace function
+} // namespace wmtk

src/wmtk/function/DifferentiablePerSimplexFunction.hpp

@@ -0,0 +1,45 @@
+#pragma once
+#include <wmtk/Primitive.hpp>
+#include <wmtk/attribute/AttributeHandle.hpp>
+#include <wmtk/attribute/MeshAttributes.hpp>
+#include "PerSimplexFunction.hpp"
+namespace wmtk {
+namespace function {
+class DifferentiablePerSimplexFunction : public PerSimplexFunction
+{
+public:
+    /**
+     * @brief Construct a new DifferentiablePerSimplexFunction object where the function is defined
+     * over simplices of simplex_type. And the differentiation is taken wrt the
+     * attribute_handle.primitive_type()
+     *
+     * @param mesh
+     * @param simplex_type
+     */
+    DifferentiablePerSimplexFunction(
+        const Mesh& mesh,
+        const PrimitiveType& simplex_type,
+        const attribute::MeshAttributeHandle<double>& variable_attribute_handle);
+    virtual ~DifferentiablePerSimplexFunction();
+
+public:
+    virtual Eigen::VectorXd get_gradient(const Simplex& s) const = 0;
+    virtual Eigen::MatrixXd get_hessian(const Simplex& s) const = 0;
+
+    const attribute::MeshAttributeHandle<double>& get_variable_attribute_handle() const;
+
+    /**
+     * @brief the function should be defined on the simplex of type A and the differntiation is
+     * taken wrt simplex type B. The definition mandate (1) A should be coface with B, (2) the type
+     * of the attribute_handle should be of the same type as B.
+     *
+     * @param s_type simplex Type B as defined above
+     */
+    void assert_function_type(const PrimitiveType& s_type) const;
+    long embedded_dimension() const;
+
+private:
+    const attribute::MeshAttributeHandle<double> m_attribute_handle;
+};
+} // namespace function
+} // namespace wmtk

src/wmtk/function/Function.cpp

@@ -1,30 +1,14 @@
 #include "Function.hpp"
 #include <wmtk/simplex/top_level_cofaces.hpp>
-namespace wmtk::function {
+namespace wmtk {
+namespace function {
 
-Function::Function(const Mesh& mesh)
-    : m_mesh(mesh)
+Function::Function(std::unique_ptr<PerSimplexFunction>&& function)
+    : m_function(std::move(function))
 {}
 
 Function::~Function() = default;
 
-double Function::get_one_ring_value(const Tuple& vertex) const
-{
-    auto simplices =
-        simplex::top_level_cofaces_tuples(mesh(), Simplex(PrimitiveType::Vertex, vertex));
-    return get_value_sum(simplices);
-}
-double Function::get_value_sum(const std::vector<Tuple>& top_level_simplices) const
-{
-    double v = 0;
-    for (const Tuple& cell : top_level_simplices) {
-        v += get_value(cell);
-    }
-    return v;
-}
 
-const Mesh& Function::mesh() const
-{
-    return m_mesh;
-}
-} // namespace wmtk::function
+} // namespace function
+} // namespace wmtk

src/wmtk/function/Function.hpp

@@ -1,25 +1,21 @@
 #pragma once
-#include <wmtk/Mesh.hpp>
-#include <wmtk/Tuple.hpp>
+#include <memory>
+#include <wmtk/Primitive.hpp>
+#include "PerSimplexFunction.hpp"
 namespace wmtk {
 namespace function {
 class Function
 {
-private:
-    const Mesh& m_mesh;
-
-
 public:
-    Function(const Mesh& mesh);
+    Function(std::unique_ptr<PerSimplexFunction>&& function);
     virtual ~Function();
 
-    const Mesh& mesh() const;
-    double get_one_ring_value(const Tuple& vertex) const;
-    double get_value_sum(const std::vector<Tuple>& top_level_simplices) const;
-
 public:
     // evaluate the function on the top level simplex of the tuple
-    virtual double get_value(const Tuple& top_level_simplex) const = 0;
+    virtual double get_value(const Simplex& simplex) const = 0;
+
+private:
+    std::unique_ptr<PerSimplexFunction> m_function;
 };
 } // namespace function
 } // namespace wmtk