DAE constraint restructure to full mass-matrix row replacement

include/micm/constraint/constraint.hpp

@@ -2,61 +2,83 @@
 // SPDX-License-Identifier: Apache-2.0
 #pragma once
 
+#include <micm/constraint/equilibrium_constraint.hpp>
+#include <micm/constraint/linear_constraint.hpp>
+
 #include <cstddef>
 #include <string>
+#include <utility>
+#include <variant>
 #include <vector>
+#include <concepts>
+#include <type_traits>
 
 namespace micm
 {
 
-  /// @brief Abstract base class for algebraic constraints G(y) = 0
-  ///
-  /// Constraints define algebraic relations that must be satisfied by the species
-  /// concentrations. They are used in DAE (Differential-Algebraic Equation) solvers
-  /// to enforce conditions like chemical equilibrium or mass conservation.
-  ///
-  /// Each constraint provides:
-  /// - A residual function G(y) that should equal zero when the constraint is satisfied
-  /// - Jacobian entries dG/dy for each species the constraint depends on
+  /// @brief This class uses std::variant to hold different constraint types.
+  ///        Each constraint provides:
+  ///        - A residual function G(y) that should equal zero when the constraint is satisfied
+  ///        - Jacobian entries dG/dy for each species the constraint depends on
   class Constraint
   {
    public:
-    /// @brief Name of the constraint (for identification)
-    std::string name_;
+    using ConstraintVariant = std::variant<EquilibriumConstraint, LinearConstraint>;
 
-    /// @brief Names of species this constraint depends on
-    std::vector<std::string> species_dependencies_;
+    ConstraintVariant constraint_;
 
-    /// @brief Default constructor
-    Constraint() = default;
+    template<typename T>
+      requires std::same_as<std::decay_t<T>, EquilibriumConstraint> || 
+               std::same_as<std::decay_t<T>, LinearConstraint>
+    Constraint(T&& constraint)
+        : constraint_(std::forward<T>(constraint))
+    {
+    }
 
-    /// @brief Constructor with name
-    /// @param name Constraint identifier
-    explicit Constraint(const std::string& name)
-        : name_(name)
+    /// @brief Get the constraint name
+    /// @return Constraint name
+    std::string GetName() const
     {
+      return std::visit([](const auto& c) { return c.name_; }, constraint_);
     }
 
-    /// @brief Virtual destructor
-    virtual ~Constraint() = default;
+    /// @brief Get species dependencies
+    /// @return Vector of species names this constraint depends on
+    const std::vector<std::string>& GetSpeciesDependencies() const
+    {
+      return std::visit([](const auto& c) -> const std::vector<std::string>& { return c.species_dependencies_; }, constraint_);
+    }
+
+    /// @brief Get the number of species this constraint depends on
+    /// @return Number of dependent species
+    std::size_t NumberOfDependencies() const
+    {
+      return std::visit([](const auto& c) { return c.species_dependencies_.size(); }, constraint_);
+    }
 
     /// @brief Evaluate the constraint residual G(y)
     /// @param concentrations Pointer to species concentrations (row of state matrix)
     /// @param indices Pointer to indices mapping species_dependencies_ to positions in concentrations
     /// @return Residual value (should be 0 when constraint is satisfied)
-    virtual double Residual(const double* concentrations, const std::size_t* indices) const = 0;
+    double Residual(const double* concentrations, const std::size_t* indices) const
+    {
+      return std::visit([&](const auto& c) { return c.Residual(concentrations, indices); }, constraint_);
+    }
 
     /// @brief Compute partial derivatives dG/d[species] for each dependent species
     /// @param concentrations Pointer to species concentrations (row of state matrix)
     /// @param indices Pointer to indices mapping species_dependencies_ to positions in concentrations
     /// @param jacobian Output buffer for partial derivatives dG/d[species] (same order as species_dependencies_)
-    virtual void Jacobian(const double* concentrations, const std::size_t* indices, double* jacobian) const = 0;
+    void Jacobian(const double* concentrations, const std::size_t* indices, double* jacobian) const
+    {
+      std::visit([&](const auto& c) { c.Jacobian(concentrations, indices, jacobian); }, constraint_);
+    }
 
-    /// @brief Get the number of species this constraint depends on
-    /// @return Number of dependent species
-    std::size_t NumberOfDependencies() const
+    /// @brief Returns the species whose state row should be replaced by this algebraic constraint
+    /// @return Algebraic species name
+    const std::string& GetAlgebraicSpecies() const
     {
-      return species_dependencies_.size();
+      return std::visit([](const auto& c) -> const std::string& { return c.AlgebraicSpecies(); }, constraint_);
     }
   };