diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml
index 32bdd4316..da5c4f04e 100755
--- a/conda.recipe/meta.yaml
+++ b/conda.recipe/meta.yaml
@@ -52,3 +52,4 @@ test:
imports:
- simcoon
- simcoon._core
+ - simcoon.modular
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/damage_mechanism.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/damage_mechanism.hpp
new file mode 100644
index 000000000..8fdf14042
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/damage_mechanism.hpp
@@ -0,0 +1,194 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file damage_mechanism.hpp
+ * @brief Scalar damage mechanism for continuum damage mechanics (CDM)
+ *
+ * This mechanism implements isotropic scalar damage following the
+ * effective stress concept:
+ * sigma_eff = sigma / (1 - D)
+ *
+ * Damage evolution is driven by an energy-based criterion:
+ * Y = (1/2) * sigma : S : sigma (damage driving force)
+ *
+ * Damage evolution follows:
+ * D = f(Y_max)
+ *
+ * where Y_max is the maximum damage driving force reached in history.
+ *
+ * Supported damage evolution laws:
+ * - Linear: D = (Y - Y_0) / (Y_c - Y_0) for Y > Y_0
+ * Props: Y_0, Y_c
+ * - Exponential: D = 1 - exp(-A * (Y - Y_0) / (Y_c - Y_0)) for Y > Y_0
+ * Props: Y_0, Y_c, A
+ * - Power law: D = ((Y - Y_0) / (Y_c - Y_0))^n for Y > Y_0
+ * Props: Y_0, Y_c, n
+ * - Weibull: D = 1 - exp(-((Y - Y_0) / A)^n) for Y > Y_0
+ * Props: Y_0, Y_c, A (scale), n (shape)
+ *
+ * @see StrainMechanism, ModularUMAT
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include
+#include
+
+namespace simcoon {
+
+/**
+ * @brief Type of damage evolution law
+ */
+enum class DamageType {
+ LINEAR = 0, ///< Linear damage evolution
+ EXPONENTIAL = 1, ///< Exponential damage evolution
+ POWER_LAW = 2, ///< Power-law damage evolution
+ WEIBULL = 3 ///< Weibull-based damage
+};
+
+/**
+ * @brief Scalar damage mechanism
+ *
+ * Implements isotropic damage following continuum damage mechanics.
+ * The damage variable D reduces the effective stiffness:
+ * L_damaged = (1 - D) * L
+ */
+class DamageMechanism final : public StrainMechanism {
+private:
+ DamageType damage_type_; ///< Type of damage evolution law
+
+ // Damage parameters
+ double Y_0_; ///< Damage threshold (no damage below this)
+ double Y_c_; ///< Critical damage driving force
+ double D_c_; ///< Critical damage value (default: 0.99)
+ double A_; ///< Damage evolution parameter
+ double n_; ///< Power law exponent
+
+ // Cached values
+ mutable double Y_current_; ///< Current damage driving force
+ mutable double D_current_; ///< Current damage
+ mutable double dD_dY_; ///< Derivative of damage w.r.t. driving force
+ mutable arma::mat M_cached_; ///< Cached compliance (set from L on first use)
+ mutable bool M_cached_valid_; ///< Whether M_cached_ is valid
+
+public:
+ /**
+ * @brief Constructor
+ * @param type Type of damage evolution law (default: LINEAR)
+ */
+ explicit DamageMechanism(DamageType type = DamageType::LINEAR);
+
+ // StrainMechanism interface
+ void configure(const arma::vec& props, int& offset) override;
+ void register_variables(InternalVariableCollection& ivc) override;
+
+ [[nodiscard]] int num_constraints() const override { return 1; }
+ [[nodiscard]] MechanismType type() const override { return MechanismType::DAMAGE; }
+
+ void compute_constraints(
+ const arma::vec& sigma,
+ const arma::mat& L,
+ double DTime,
+ const InternalVariableCollection& ivc,
+ arma::vec& Phi,
+ arma::vec& Y_crit
+ ) const override;
+
+ void compute_flow_directions(
+ const arma::vec& sigma,
+ const InternalVariableCollection& ivc,
+ std::map& Lambda_map
+ ) const override;
+
+ void compute_jacobian_contribution(
+ const arma::vec& sigma,
+ const arma::mat& L,
+ const InternalVariableCollection& ivc,
+ arma::mat& B,
+ int row_offset
+ ) const override;
+
+ arma::vec inelastic_strain(const InternalVariableCollection& ivc) const override;
+
+ void update(
+ const arma::vec& ds,
+ int offset,
+ InternalVariableCollection& ivc
+ ) override;
+
+ void tangent_contribution(
+ const arma::vec& sigma,
+ const arma::mat& L,
+ const arma::vec& Ds,
+ int offset,
+ const InternalVariableCollection& ivc,
+ arma::mat& Lt
+ ) const override;
+
+ void compute_work(
+ const arma::vec& sigma_start,
+ const arma::vec& sigma,
+ const InternalVariableCollection& ivc,
+ double& Wm_r,
+ double& Wm_ir,
+ double& Wm_d
+ ) const override;
+
+ // Damage-specific methods
+
+ /**
+ * @brief Compute damage driving force from stress
+ * @param sigma Stress (6-component Voigt)
+ * @param S Compliance tensor (6x6)
+ * @return Damage driving force Y
+ */
+ double compute_driving_force(const arma::vec& sigma, const arma::mat& S) const;
+
+ /**
+ * @brief Compute damage from driving force
+ * @param Y Damage driving force
+ * @param Y_max Maximum driving force in history
+ * @return Damage value D
+ */
+ double compute_damage(double Y, double Y_max) const;
+
+ /**
+ * @brief Get the current damage value
+ * @param ivc Internal variable collection
+ * @return Current damage D
+ */
+ double get_damage(const InternalVariableCollection& ivc) const;
+
+ /**
+ * @brief Get damaged stiffness tensor
+ * @param L Undamaged stiffness
+ * @param D Damage value
+ * @return Damaged stiffness (1-D)*L
+ */
+ static arma::mat damaged_stiffness(const arma::mat& L, double D);
+
+ // Accessors
+ [[nodiscard]] DamageType damage_type() const noexcept { return damage_type_; }
+ [[nodiscard]] double Y_0() const noexcept { return Y_0_; }
+ [[nodiscard]] double Y_c() const noexcept { return Y_c_; }
+ [[nodiscard]] double D_c() const noexcept { return D_c_; }
+};
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/elasticity_module.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/elasticity_module.hpp
new file mode 100644
index 000000000..c9e25762a
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/elasticity_module.hpp
@@ -0,0 +1,212 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file elasticity_module.hpp
+ * @brief Elasticity module for modular UMAT.
+ *
+ * This module wraps existing elasticity functions (L_iso, L_ortho, L_isotrans)
+ * to provide a configurable elasticity component.
+ *
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+
+namespace simcoon {
+
+/**
+ * @brief Types of linear elasticity
+ */
+enum class ElasticityType {
+ ISOTROPIC = 0, ///< Isotropic: E, nu, alpha
+ CUBIC = 1, ///< Cubic: E, nu, G, alpha (3 independent elastic constants)
+ TRANSVERSE_ISOTROPIC = 2, ///< Transverse isotropic: EL, ET, nuTL, nuTT, GLT, alpha_L, alpha_T
+ ORTHOTROPIC = 3 ///< Orthotropic: E1, E2, E3, nu12, nu13, nu23, G12, G13, G23, alpha1, alpha2, alpha3
+};
+
+/**
+ * @brief Elasticity module for modular UMAT
+ *
+ * This class encapsulates the computation of the elastic stiffness tensor
+ * and coefficient of thermal expansion for different types of linear elasticity.
+ */
+class ElasticityModule {
+private:
+ ElasticityType type_;
+ arma::mat L_; ///< 6x6 stiffness tensor
+ arma::mat M_; ///< 6x6 compliance tensor
+ arma::vec alpha_; ///< 6-component CTE (Voigt notation)
+ bool configured_;
+
+public:
+ /**
+ * @brief Default constructor
+ */
+ ElasticityModule();
+
+ // Default copy/move/destructor
+ ElasticityModule(const ElasticityModule&) = default;
+ ElasticityModule(ElasticityModule&&) = default;
+ ElasticityModule& operator=(const ElasticityModule&) = default;
+ ElasticityModule& operator=(ElasticityModule&&) = default;
+ ~ElasticityModule() = default;
+
+ // ========== Configuration ==========
+
+ /**
+ * @brief Configure as isotropic elasticity
+ * @param E Young's modulus
+ * @param nu Poisson's ratio
+ * @param alpha Coefficient of thermal expansion (scalar, isotropic)
+ */
+ void configure_isotropic(double E, double nu, double alpha);
+
+ /**
+ * @brief Configure as cubic elasticity
+ * @param E Young's modulus
+ * @param nu Poisson's ratio
+ * @param G Shear modulus (independent from E and nu for cubic symmetry)
+ * @param alpha Coefficient of thermal expansion (scalar, isotropic CTE)
+ *
+ * Cubic symmetry has 3 independent elastic constants. Unlike isotropic
+ * materials where G = E/(2(1+nu)), G is an independent parameter here.
+ * The Zener anisotropy ratio A = 2*G*(1+nu)/E quantifies the deviation
+ * from isotropy (A=1 for isotropic).
+ *
+ * Alternative: use configure_cubic_Cii() with C11, C12, C44 directly.
+ */
+ void configure_cubic(double E, double nu, double G, double alpha);
+
+ /**
+ * @brief Configure as cubic elasticity using Cii convention
+ * @param C11 Stiffness component C11
+ * @param C12 Stiffness component C12
+ * @param C44 Stiffness component C44
+ * @param alpha Coefficient of thermal expansion (scalar, isotropic CTE)
+ */
+ void configure_cubic_Cii(double C11, double C12, double C44, double alpha);
+
+ /**
+ * @brief Configure as transversely isotropic elasticity
+ * @param EL Longitudinal Young's modulus
+ * @param ET Transverse Young's modulus
+ * @param nuTL Poisson's ratio (transverse-longitudinal)
+ * @param nuTT Poisson's ratio (transverse-transverse)
+ * @param GLT Shear modulus
+ * @param alpha_L Longitudinal CTE
+ * @param alpha_T Transverse CTE
+ * @param axis Axis of symmetry (1=x, 2=y, 3=z)
+ */
+ void configure_transverse_isotropic(double EL, double ET, double nuTL, double nuTT,
+ double GLT, double alpha_L, double alpha_T, int axis = 3);
+
+ /**
+ * @brief Configure as orthotropic elasticity
+ * @param E1 Young's modulus in direction 1
+ * @param E2 Young's modulus in direction 2
+ * @param E3 Young's modulus in direction 3
+ * @param nu12 Poisson's ratio 1-2
+ * @param nu13 Poisson's ratio 1-3
+ * @param nu23 Poisson's ratio 2-3
+ * @param G12 Shear modulus 1-2
+ * @param G13 Shear modulus 1-3
+ * @param G23 Shear modulus 2-3
+ * @param alpha1 CTE in direction 1
+ * @param alpha2 CTE in direction 2
+ * @param alpha3 CTE in direction 3
+ */
+ void configure_orthotropic(double E1, double E2, double E3,
+ double nu12, double nu13, double nu23,
+ double G12, double G13, double G23,
+ double alpha1, double alpha2, double alpha3);
+
+ /**
+ * @brief Configure from props array
+ * @param type Elasticity type
+ * @param props Material properties vector
+ * @param offset Current offset in props (will be updated)
+ */
+ void configure(ElasticityType type, const arma::vec& props, int& offset);
+
+ // ========== Accessors ==========
+
+ /**
+ * @brief Get the elasticity type
+ * @return The configured type
+ */
+ [[nodiscard]] ElasticityType type() const noexcept { return type_; }
+
+ /**
+ * @brief Check if the module is configured
+ * @return True if configure has been called
+ */
+ [[nodiscard]] bool is_configured() const noexcept { return configured_; }
+
+ /**
+ * @brief Get the 6x6 stiffness tensor
+ * @return Const reference to L
+ */
+ [[nodiscard]] const arma::mat& L() const noexcept { return L_; }
+
+ /**
+ * @brief Get the 6x6 compliance tensor
+ * @return Const reference to M
+ */
+ [[nodiscard]] const arma::mat& M() const noexcept { return M_; }
+
+ /**
+ * @brief Get the CTE tensor
+ * @return Const reference to alpha (6 components)
+ */
+ [[nodiscard]] const arma::vec& alpha() const noexcept { return alpha_; }
+
+ // ========== Derived Quantities ==========
+
+ /**
+ * @brief Get damaged stiffness tensor
+ * @param d Damage variable (0 = undamaged, 1 = fully damaged)
+ * @return Degraded stiffness (1-d)*L
+ */
+ arma::mat damaged_L(double d) const;
+
+ /**
+ * @brief Get thermal strain
+ * @param DT Temperature increment
+ * @return Thermal strain vector (alpha * DT)
+ */
+ arma::vec thermal_strain(double DT) const;
+
+ /**
+ * @brief Get number of props consumed by this elasticity type
+ * @param type The elasticity type
+ * @return Number of properties required
+ */
+ [[nodiscard]] static constexpr int props_count(ElasticityType type) {
+ switch (type) {
+ case ElasticityType::ISOTROPIC: return 3;
+ case ElasticityType::CUBIC: return 4;
+ case ElasticityType::TRANSVERSE_ISOTROPIC: return 8;
+ case ElasticityType::ORTHOTROPIC: return 12;
+ }
+ return 0;
+ }
+};
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/hardening.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/hardening.hpp
new file mode 100644
index 000000000..360c4640c
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/hardening.hpp
@@ -0,0 +1,370 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file hardening.hpp
+ * @brief Hardening modules for modular UMAT.
+ *
+ * Provides isotropic and kinematic hardening models:
+ * - Isotropic: Linear, Power-law, Voce, Combined Voce
+ * - Kinematic: Prager (linear), Armstrong-Frederick, Chaboche (multi-AF)
+ *
+ * @see PlasticityMechanism, ModularUMAT
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include
+#include
+#include
+#include
+
+namespace simcoon {
+
+// ============================================================================
+// ISOTROPIC HARDENING
+// ============================================================================
+
+/**
+ * @brief Types of isotropic hardening
+ */
+enum class IsoHardType {
+ NONE = 0, ///< No isotropic hardening
+ LINEAR = 1, ///< Linear: R = H * p
+ POWER_LAW = 2, ///< Power law: R = k * p^m
+ VOCE = 3, ///< Voce saturation: R = Q * (1 - exp(-b*p))
+ COMBINED_VOCE = 4 ///< Sum of N Voce terms
+};
+
+/**
+ * @brief Base class for isotropic hardening
+ */
+class IsotropicHardening {
+public:
+ virtual ~IsotropicHardening() = default;
+
+ /**
+ * @brief Configure from props array
+ * @param props Material properties vector
+ * @param offset Current offset in props (will be updated)
+ */
+ virtual void configure(const arma::vec& props, int& offset) = 0;
+
+ /**
+ * @brief Compute hardening stress R(p)
+ * @param p Accumulated plastic strain
+ * @return Hardening stress
+ */
+ [[nodiscard]] virtual double R(double p) const = 0;
+
+ /**
+ * @brief Compute hardening modulus dR/dp
+ * @param p Accumulated plastic strain
+ * @return Derivative of hardening stress
+ */
+ [[nodiscard]] virtual double dR_dp(double p) const = 0;
+
+ /**
+ * @brief Get the hardening type
+ * @return Hardening type enum
+ */
+ [[nodiscard]] virtual IsoHardType type() const = 0;
+
+ /**
+ * @brief Create an isotropic hardening instance
+ * @param type Hardening type
+ * @param N Number of terms (for COMBINED_VOCE)
+ * @return Unique pointer to hardening instance
+ */
+ static std::unique_ptr create(IsoHardType type, int N = 1);
+
+ /**
+ * @brief Get number of props consumed by this hardening type
+ * @param type The hardening type
+ * @param N Number of terms (for COMBINED_VOCE)
+ * @return Number of properties required
+ */
+ [[nodiscard]] static constexpr int props_count(IsoHardType type, int N = 1) {
+ switch (type) {
+ case IsoHardType::NONE: return 0;
+ case IsoHardType::LINEAR: return 1;
+ case IsoHardType::POWER_LAW: return 2;
+ case IsoHardType::VOCE: return 2;
+ case IsoHardType::COMBINED_VOCE: return 2 * N;
+ }
+ return 0;
+ }
+};
+
+/**
+ * @brief No isotropic hardening
+ */
+class NoIsotropicHardening final : public IsotropicHardening {
+public:
+ void configure(const arma::vec& props, int& offset) override {}
+ double R(double p) const override { return 0.0; }
+ double dR_dp(double p) const override { return 0.0; }
+ IsoHardType type() const override { return IsoHardType::NONE; }
+};
+
+/**
+ * @brief Linear isotropic hardening: R = H * p
+ */
+class LinearHardening final : public IsotropicHardening {
+private:
+ double H_; ///< Hardening modulus
+public:
+ LinearHardening() : H_(0.0) {}
+ void configure(const arma::vec& props, int& offset) override;
+ double R(double p) const override { return H_ * p; }
+ double dR_dp(double p) const override { return H_; }
+ IsoHardType type() const override { return IsoHardType::LINEAR; }
+};
+
+/**
+ * @brief Power-law isotropic hardening: R = k * p^m
+ */
+class PowerLawHardening final : public IsotropicHardening {
+private:
+ double k_; ///< Hardening coefficient
+ double m_; ///< Hardening exponent
+public:
+ PowerLawHardening() : k_(0.0), m_(1.0) {}
+ void configure(const arma::vec& props, int& offset) override;
+ double R(double p) const override;
+ double dR_dp(double p) const override;
+ IsoHardType type() const override { return IsoHardType::POWER_LAW; }
+};
+
+/**
+ * @brief Voce saturation hardening: R = Q * (1 - exp(-b*p))
+ */
+class VoceHardening final : public IsotropicHardening {
+private:
+ double Q_; ///< Saturation stress
+ double b_; ///< Hardening rate
+public:
+ VoceHardening() : Q_(0.0), b_(0.0) {}
+ void configure(const arma::vec& props, int& offset) override;
+ double R(double p) const override;
+ double dR_dp(double p) const override;
+ IsoHardType type() const override { return IsoHardType::VOCE; }
+};
+
+/**
+ * @brief Combined Voce hardening: R = sum_i Q_i * (1 - exp(-b_i*p))
+ */
+class CombinedVoceHardening final : public IsotropicHardening {
+private:
+ int N_; ///< Number of Voce terms
+ arma::vec Q_; ///< Saturation stresses
+ arma::vec b_; ///< Hardening rates
+public:
+ explicit CombinedVoceHardening(int N = 1) : N_(N), Q_(arma::zeros(N)), b_(arma::zeros(N)) {}
+ void configure(const arma::vec& props, int& offset) override;
+ double R(double p) const override;
+ double dR_dp(double p) const override;
+ IsoHardType type() const override { return IsoHardType::COMBINED_VOCE; }
+ int num_terms() const { return N_; }
+};
+
+// ============================================================================
+// KINEMATIC HARDENING
+// ============================================================================
+
+/**
+ * @brief Types of kinematic hardening
+ */
+enum class KinHardType {
+ NONE = 0, ///< No kinematic hardening
+ PRAGER = 1, ///< Linear Prager: X = (2/3)*C*alpha
+ ARMSTRONG_FREDERICK = 2, ///< Single AF: dX = (2/3)*C*dep - D*X*dp
+ CHABOCHE = 3 ///< Multiple AF terms
+};
+
+/**
+ * @brief Base class for kinematic hardening
+ */
+class KinematicHardening {
+public:
+ virtual ~KinematicHardening() = default;
+
+ /**
+ * @brief Configure from props array
+ * @param props Material properties vector
+ * @param offset Current offset in props (will be updated)
+ */
+ virtual void configure(const arma::vec& props, int& offset) = 0;
+
+ /**
+ * @brief Register internal variables for this hardening model
+ * @param ivc Internal variable collection
+ */
+ virtual void register_variables(InternalVariableCollection& ivc) = 0;
+
+ /**
+ * @brief Get total backstress X = sum(X_i)
+ * @param ivc Internal variable collection
+ * @return Total backstress tensor (6 components)
+ */
+ virtual arma::vec total_backstress(const InternalVariableCollection& ivc) const = 0;
+
+ /**
+ * @brief Compute flow direction for kinematic variable alpha_i
+ * @param i Index of the backstress (0 for single term)
+ * @param n Plastic flow direction
+ * @param ivc Internal variable collection
+ * @return Flow direction for alpha_i
+ */
+ virtual arma::vec alpha_flow(int i, const arma::vec& n, const InternalVariableCollection& ivc) const = 0;
+
+ /**
+ * @brief Compute total kinematic hardening modulus contribution
+ * @param n Plastic flow direction
+ * @param ivc Internal variable collection
+ * @return Hardening modulus contribution to K
+ */
+ virtual double hardening_modulus(const arma::vec& n, const InternalVariableCollection& ivc) const = 0;
+
+ /**
+ * @brief Update internal variables given plastic multiplier increment
+ * @param dp Plastic multiplier increment
+ * @param n Plastic flow direction
+ * @param ivc Internal variable collection
+ */
+ virtual void update(double dp, const arma::vec& n, InternalVariableCollection& ivc) = 0;
+
+ /**
+ * @brief Get the hardening type
+ * @return Hardening type enum
+ */
+ [[nodiscard]] virtual KinHardType type() const = 0;
+
+ /**
+ * @brief Get number of backstress terms
+ * @return Number of backstress tensors
+ */
+ [[nodiscard]] virtual int num_backstresses() const = 0;
+
+ /**
+ * @brief Create a kinematic hardening instance
+ * @param type Hardening type
+ * @param N Number of backstress terms (for CHABOCHE)
+ * @return Unique pointer to hardening instance
+ */
+ static std::unique_ptr create(KinHardType type, int N = 1);
+
+ /**
+ * @brief Get number of props consumed by this hardening type
+ * @param type The hardening type
+ * @param N Number of backstress terms (for CHABOCHE)
+ * @return Number of properties required
+ */
+ [[nodiscard]] static constexpr int props_count(KinHardType type, int N = 1) {
+ switch (type) {
+ case KinHardType::NONE: return 0;
+ case KinHardType::PRAGER: return 1;
+ case KinHardType::ARMSTRONG_FREDERICK: return 2;
+ case KinHardType::CHABOCHE: return 2 * N;
+ }
+ return 0;
+ }
+};
+
+/**
+ * @brief No kinematic hardening
+ */
+class NoKinematicHardening final : public KinematicHardening {
+public:
+ void configure(const arma::vec& props, int& offset) override {}
+ void register_variables(InternalVariableCollection& ivc) override {}
+ arma::vec total_backstress(const InternalVariableCollection& ivc) const override {
+ return arma::zeros(6);
+ }
+ arma::vec alpha_flow(int i, const arma::vec& n, const InternalVariableCollection& ivc) const override {
+ return arma::zeros(6);
+ }
+ double hardening_modulus(const arma::vec& n, const InternalVariableCollection& ivc) const override {
+ return 0.0;
+ }
+ void update(double dp, const arma::vec& n, InternalVariableCollection& ivc) override {}
+ KinHardType type() const override { return KinHardType::NONE; }
+ int num_backstresses() const override { return 0; }
+};
+
+/**
+ * @brief Linear Prager kinematic hardening: X = (2/3)*C*alpha
+ */
+class PragerHardening final : public KinematicHardening {
+private:
+ double C_; ///< Kinematic hardening modulus
+public:
+ PragerHardening() : C_(0.0) {}
+ void configure(const arma::vec& props, int& offset) override;
+ void register_variables(InternalVariableCollection& ivc) override;
+ arma::vec total_backstress(const InternalVariableCollection& ivc) const override;
+ arma::vec alpha_flow(int i, const arma::vec& n, const InternalVariableCollection& ivc) const override;
+ double hardening_modulus(const arma::vec& n, const InternalVariableCollection& ivc) const override;
+ void update(double dp, const arma::vec& n, InternalVariableCollection& ivc) override;
+ KinHardType type() const override { return KinHardType::PRAGER; }
+ int num_backstresses() const override { return 1; }
+};
+
+/**
+ * @brief Armstrong-Frederick kinematic hardening: dX = (2/3)*C*dep - D*X*dp
+ */
+class ArmstrongFrederickHardening final : public KinematicHardening {
+private:
+ double C_; ///< Hardening parameter
+ double D_; ///< Dynamic recovery parameter
+public:
+ ArmstrongFrederickHardening() : C_(0.0), D_(0.0) {}
+ void configure(const arma::vec& props, int& offset) override;
+ void register_variables(InternalVariableCollection& ivc) override;
+ arma::vec total_backstress(const InternalVariableCollection& ivc) const override;
+ arma::vec alpha_flow(int i, const arma::vec& n, const InternalVariableCollection& ivc) const override;
+ double hardening_modulus(const arma::vec& n, const InternalVariableCollection& ivc) const override;
+ void update(double dp, const arma::vec& n, InternalVariableCollection& ivc) override;
+ KinHardType type() const override { return KinHardType::ARMSTRONG_FREDERICK; }
+ int num_backstresses() const override { return 1; }
+};
+
+/**
+ * @brief Chaboche kinematic hardening: dX_i = (2/3)*C_i*dep - D_i*X_i*dp
+ *
+ * Multiple Armstrong-Frederick terms for capturing different hardening time scales.
+ */
+class ChabocheHardening final : public KinematicHardening {
+private:
+ int N_; ///< Number of backstress terms
+ arma::vec C_; ///< Hardening parameters
+ arma::vec D_; ///< Dynamic recovery parameters
+public:
+ explicit ChabocheHardening(int N = 1) : N_(N), C_(arma::zeros(N)), D_(arma::zeros(N)) {}
+ void configure(const arma::vec& props, int& offset) override;
+ void register_variables(InternalVariableCollection& ivc) override;
+ arma::vec total_backstress(const InternalVariableCollection& ivc) const override;
+ arma::vec alpha_flow(int i, const arma::vec& n, const InternalVariableCollection& ivc) const override;
+ double hardening_modulus(const arma::vec& n, const InternalVariableCollection& ivc) const override;
+ void update(double dp, const arma::vec& n, InternalVariableCollection& ivc) override;
+ KinHardType type() const override { return KinHardType::CHABOCHE; }
+ int num_backstresses() const override { return N_; }
+};
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable.hpp
new file mode 100644
index 000000000..790b9b6f1
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable.hpp
@@ -0,0 +1,280 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file internal_variable.hpp
+ * @brief Type-safe internal variable for constitutive models.
+ *
+ * This class provides storage for internal variables used in constitutive laws
+ * (plasticity, viscoelasticity, damage, etc.) with support for:
+ * - Scalar variables (accumulated plastic strain, damage, etc.)
+ * - 6-component vectors in Voigt notation (plastic strain, backstress)
+ * - 6x6 matrices (for complex tensorial internal variables)
+ *
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include
+#include
+
+namespace simcoon {
+
+/**
+ * @brief Type enumeration for internal variables
+ */
+enum class IVarType {
+ SCALAR, ///< Single scalar value (double)
+ VECTOR_6, ///< 6-component vector (Voigt notation for symmetric tensors)
+ MATRIX_6x6 ///< 6x6 matrix (for fourth-order tensors in Voigt notation)
+};
+
+/**
+ * @brief Type-safe internal variable with automatic serialization
+ *
+ * This class encapsulates internal variables used in constitutive models,
+ * providing type safety, automatic pack/unpack to statev arrays, and
+ * support for objectivity through rotation operations.
+ */
+class InternalVariable {
+private:
+ std::string name_; ///< Identifier for debugging and access
+ IVarType type_; ///< Storage type
+
+ // Storage (only one active based on type_)
+ double scalar_value_; ///< Current scalar value
+ double scalar_start_; ///< Scalar value at start of increment
+ arma::vec vec_value_; ///< Current vector value (6 components)
+ arma::vec vec_start_; ///< Vector value at start of increment
+ arma::mat mat_value_; ///< Current matrix value (6x6)
+ arma::mat mat_start_; ///< Matrix value at start of increment
+
+ bool requires_rotation_; ///< Whether to apply rotation for objectivity
+ unsigned int statev_offset_; ///< Position in flat statev array
+
+public:
+ /**
+ * @brief Construct a scalar internal variable
+ * @param name Identifier for the variable
+ * @param init Initial value (default: 0.0)
+ * @param rotate Whether rotation should be applied (default: false for scalars)
+ */
+ InternalVariable(const std::string& name, double init = 0.0, bool rotate = false);
+
+ /**
+ * @brief Construct a vector internal variable (6 Voigt components)
+ * @param name Identifier for the variable
+ * @param init Initial value (default: zeros)
+ * @param rotate Whether rotation should be applied (default: true for tensors)
+ */
+ InternalVariable(const std::string& name, const arma::vec& init, bool rotate = true);
+
+ /**
+ * @brief Construct a matrix internal variable (6x6)
+ * @param name Identifier for the variable
+ * @param init Initial value (default: zeros)
+ * @param rotate Whether rotation should be applied (default: true for tensors)
+ */
+ InternalVariable(const std::string& name, const arma::mat& init, bool rotate = true);
+
+ // Default copy/move/destructor
+ InternalVariable(const InternalVariable&) = default;
+ InternalVariable(InternalVariable&&) = default;
+ InternalVariable& operator=(const InternalVariable&) = default;
+ InternalVariable& operator=(InternalVariable&&) = default;
+ ~InternalVariable() = default;
+
+ // ========== Type Information ==========
+
+ /**
+ * @brief Get the variable type
+ * @return The IVarType of this variable
+ */
+ [[nodiscard]] IVarType type() const noexcept { return type_; }
+
+ /**
+ * @brief Get the variable name
+ * @return Reference to the name string
+ */
+ [[nodiscard]] const std::string& name() const noexcept { return name_; }
+
+ /**
+ * @brief Get the number of scalar values stored
+ * @return 1 for scalar, 6 for vector, 36 for matrix
+ */
+ unsigned int size() const;
+
+ /**
+ * @brief Check if this variable requires rotation for objectivity
+ * @return True if rotation should be applied
+ */
+ [[nodiscard]] bool requires_rotation() const noexcept { return requires_rotation_; }
+
+ // ========== Value Accessors (throw if wrong type) ==========
+
+ /**
+ * @brief Access scalar value (mutable)
+ * @return Reference to scalar value
+ * @throws std::runtime_error if type is not SCALAR
+ */
+ double& scalar();
+
+ /**
+ * @brief Access scalar value (const)
+ * @return Const reference to scalar value
+ * @throws std::runtime_error if type is not SCALAR
+ */
+ const double& scalar() const;
+
+ /**
+ * @brief Access vector value (mutable)
+ * @return Reference to vector value
+ * @throws std::runtime_error if type is not VECTOR_6
+ */
+ arma::vec& vec();
+
+ /**
+ * @brief Access vector value (const)
+ * @return Const reference to vector value
+ * @throws std::runtime_error if type is not VECTOR_6
+ */
+ const arma::vec& vec() const;
+
+ /**
+ * @brief Access matrix value (mutable)
+ * @return Reference to matrix value
+ * @throws std::runtime_error if type is not MATRIX_6x6
+ */
+ arma::mat& mat();
+
+ /**
+ * @brief Access matrix value (const)
+ * @return Const reference to matrix value
+ * @throws std::runtime_error if type is not MATRIX_6x6
+ */
+ const arma::mat& mat() const;
+
+ // ========== Start Value Accessors ==========
+
+ /**
+ * @brief Access scalar start value
+ * @return Reference to scalar start value
+ */
+ double& scalar_start();
+ const double& scalar_start() const;
+
+ /**
+ * @brief Access vector start value
+ * @return Reference to vector start value
+ */
+ arma::vec& vec_start();
+ const arma::vec& vec_start() const;
+
+ /**
+ * @brief Access matrix start value
+ * @return Reference to matrix start value
+ */
+ arma::mat& mat_start();
+ const arma::mat& mat_start() const;
+
+ // ========== Increment Computation ==========
+
+ /**
+ * @brief Compute scalar increment (current - start)
+ * @return The increment value
+ */
+ double delta_scalar() const;
+
+ /**
+ * @brief Compute vector increment (current - start)
+ * @return The increment vector
+ */
+ arma::vec delta_vec() const;
+
+ /**
+ * @brief Compute matrix increment (current - start)
+ * @return The increment matrix
+ */
+ arma::mat delta_mat() const;
+
+ // ========== State Management ==========
+
+ /**
+ * @brief Copy current value to start value
+ *
+ * Call at the beginning of an increment to save the state.
+ */
+ void to_start();
+
+ /**
+ * @brief Copy start value to current value
+ *
+ * Call to restore state to beginning of increment.
+ */
+ void set_start();
+
+ /**
+ * @brief Apply rotation for objectivity
+ * @param DR Rotation increment matrix (3x3)
+ *
+ * Applies rotation to maintain objectivity during large deformations.
+ * - Scalar variables: not rotated (scalars are frame-invariant)
+ * - Vector variables (6 Voigt): rotated via rotate_strain(value, DR)
+ * - Matrix variables (6x6): not yet supported (placeholder)
+ *
+ * Only applied if requires_rotation_ is true. By default, scalars
+ * have requires_rotation_ = false and tensorial quantities have
+ * requires_rotation_ = true.
+ *
+ * @see rotate_strain() in contimech.hpp
+ */
+ void rotate(const arma::mat& DR);
+
+ // ========== Serialization ==========
+
+ /**
+ * @brief Set the offset in the statev array
+ * @param offset Starting index in statev
+ */
+ void set_offset(unsigned int offset) { statev_offset_ = offset; }
+
+ /**
+ * @brief Get the offset in the statev array
+ * @return The offset value
+ */
+ [[nodiscard]] unsigned int offset() const noexcept { return statev_offset_; }
+
+ /**
+ * @brief Pack current value into statev array
+ * @param statev The state variable vector to pack into
+ *
+ * Writes size() values starting at statev_offset_.
+ */
+ void pack(arma::vec& statev) const;
+
+ /**
+ * @brief Unpack value from statev array
+ * @param statev The state variable vector to unpack from
+ *
+ * Reads size() values starting at statev_offset_.
+ */
+ void unpack(const arma::vec& statev);
+};
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable_collection.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable_collection.hpp
new file mode 100644
index 000000000..f275ddab9
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/internal_variable_collection.hpp
@@ -0,0 +1,246 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file internal_variable_collection.hpp
+ * @brief Collection of internal variables for modular constitutive models.
+ *
+ * This class manages multiple InternalVariable instances, providing:
+ * - Named registration and access
+ * - Automatic offset computation for statev serialization
+ * - Bulk operations (pack/unpack, rotate, state management)
+ *
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace simcoon {
+
+/**
+ * @brief Manages a collection of internal variables for constitutive models
+ *
+ * This class provides a container for InternalVariable instances with
+ * named access, automatic offset computation for statev serialization,
+ * and bulk operations for common tasks.
+ */
+class InternalVariableCollection {
+private:
+ std::vector> variables_;
+ std::unordered_map name_to_index_;
+ unsigned int total_size_;
+ bool offsets_computed_;
+
+public:
+ /**
+ * @brief Default constructor
+ */
+ InternalVariableCollection();
+
+ // Disable copy (due to unique_ptr)
+ InternalVariableCollection(const InternalVariableCollection&) = delete;
+ InternalVariableCollection& operator=(const InternalVariableCollection&) = delete;
+
+ // Allow move
+ InternalVariableCollection(InternalVariableCollection&&) = default;
+ InternalVariableCollection& operator=(InternalVariableCollection&&) = default;
+
+ ~InternalVariableCollection() = default;
+
+ // ========== Registration ==========
+
+ /**
+ * @brief Add a scalar internal variable
+ * @param name Unique identifier for the variable
+ * @param init Initial value (default: 0.0)
+ * @param rotate Whether rotation should be applied (default: false)
+ * @return Reference to the created variable
+ * @throws std::runtime_error if name already exists
+ */
+ InternalVariable& add_scalar(const std::string& name, double init = 0.0, bool rotate = false);
+
+ /**
+ * @brief Add a 6-component vector internal variable
+ * @param name Unique identifier for the variable
+ * @param init Initial value (default: zeros)
+ * @param rotate Whether rotation should be applied (default: true)
+ * @return Reference to the created variable
+ * @throws std::runtime_error if name already exists
+ */
+ InternalVariable& add_vec(const std::string& name, const arma::vec& init = arma::zeros(6), bool rotate = true);
+
+ /**
+ * @brief Add a 6x6 matrix internal variable
+ * @param name Unique identifier for the variable
+ * @param init Initial value (default: zeros)
+ * @param rotate Whether rotation should be applied (default: true)
+ * @return Reference to the created variable
+ * @throws std::runtime_error if name already exists
+ */
+ InternalVariable& add_mat(const std::string& name, const arma::mat& init = arma::zeros(6, 6), bool rotate = true);
+
+ // ========== Access ==========
+
+ /**
+ * @brief Get a variable by name
+ * @param name The variable name
+ * @return Reference to the variable
+ * @throws std::runtime_error if name not found
+ */
+ InternalVariable& get(const std::string& name);
+
+ /**
+ * @brief Get a variable by name (const)
+ * @param name The variable name
+ * @return Const reference to the variable
+ * @throws std::runtime_error if name not found
+ */
+ const InternalVariable& get(const std::string& name) const;
+
+ /**
+ * @brief Check if a variable exists
+ * @param name The variable name
+ * @return True if the variable exists
+ */
+ bool has(const std::string& name) const;
+
+ /**
+ * @brief Get variable by index
+ * @param i Index in the collection
+ * @return Reference to the variable
+ */
+ InternalVariable& operator[](size_t i);
+
+ /**
+ * @brief Get variable by index (const)
+ * @param i Index in the collection
+ * @return Const reference to the variable
+ */
+ const InternalVariable& operator[](size_t i) const;
+
+ /**
+ * @brief Get the number of variables
+ * @return Number of registered variables
+ */
+ [[nodiscard]] size_t size() const noexcept { return variables_.size(); }
+
+ /**
+ * @brief Check if collection is empty
+ * @return True if no variables are registered
+ */
+ [[nodiscard]] bool empty() const noexcept { return variables_.empty(); }
+
+ // ========== Offset Management ==========
+
+ /**
+ * @brief Get the total size needed in statev
+ * @return Total number of scalar values across all variables
+ */
+ [[nodiscard]] unsigned int total_statev_size() const noexcept { return total_size_; }
+
+ /**
+ * @brief Compute offsets for all variables
+ * @param base_offset Starting offset in statev (default: 0)
+ *
+ * This assigns sequential offsets to each variable based on their sizes.
+ * Must be called before pack_all/unpack_all.
+ */
+ void compute_offsets(unsigned int base_offset = 0);
+
+ /**
+ * @brief Check if offsets have been computed
+ * @return True if compute_offsets has been called
+ */
+ [[nodiscard]] bool offsets_computed() const noexcept { return offsets_computed_; }
+
+ // ========== Bulk Operations ==========
+
+ /**
+ * @brief Pack all variables into statev array
+ * @param statev The state variable vector to pack into
+ *
+ * Requires compute_offsets() to have been called first.
+ */
+ void pack_all(arma::vec& statev) const;
+
+ /**
+ * @brief Unpack all variables from statev array
+ * @param statev The state variable vector to unpack from
+ *
+ * Requires compute_offsets() to have been called first.
+ * Also sets start values to current values.
+ */
+ void unpack_all(const arma::vec& statev);
+
+ /**
+ * @brief Apply rotation to all variables that require it
+ * @param DR Rotation increment matrix (3x3)
+ */
+ void rotate_all(const arma::mat& DR);
+
+ /**
+ * @brief Copy current values to start values for all variables
+ */
+ void to_start_all();
+
+ /**
+ * @brief Copy start values to current values for all variables
+ */
+ void set_start_all();
+
+ /**
+ * @brief Reset all variables to zero
+ */
+ void reset_all();
+
+ // ========== Iteration ==========
+
+ /**
+ * @brief Iterator type for range-based for loops
+ */
+ using iterator = std::vector>::iterator;
+ using const_iterator = std::vector>::const_iterator;
+
+ iterator begin() { return variables_.begin(); }
+ iterator end() { return variables_.end(); }
+ const_iterator begin() const { return variables_.begin(); }
+ const_iterator end() const { return variables_.end(); }
+ const_iterator cbegin() const { return variables_.cbegin(); }
+ const_iterator cend() const { return variables_.cend(); }
+
+ // ========== Utility ==========
+
+ /**
+ * @brief Get a list of all variable names
+ * @return Vector of variable names in registration order
+ */
+ std::vector names() const;
+
+ /**
+ * @brief Clear all variables
+ */
+ void clear();
+};
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/modular_umat.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/modular_umat.hpp
new file mode 100644
index 000000000..59e2496d9
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/modular_umat.hpp
@@ -0,0 +1,388 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file modular_umat.hpp
+ * @brief Modular UMAT orchestrator for composing constitutive models.
+ *
+ * This class orchestrates the composition of:
+ * - Elasticity module (isotropic, cubic, transverse isotropic, orthotropic)
+ * - Multiple strain mechanisms (plasticity, viscoelasticity, damage)
+ *
+ * It provides a unified return mapping algorithm that handles all
+ * coupled constraints from the different mechanisms.
+ *
+ * Usage example (isotropic elasticity + von Mises plasticity with Voce hardening):
+ * @code
+ * ModularUMAT mumat;
+ * int offset = 0;
+ * arma::vec props = {210000, 0.3, 1.2e-5, 300, 100, 10};
+ * mumat.set_elasticity(ElasticityType::ISOTROPIC, props, offset);
+ * mumat.add_plasticity(YieldType::VON_MISES, IsoHardType::VOCE,
+ * KinHardType::NONE, 1, 1, props, offset);
+ * mumat.initialize(nstatev, statev);
+ * mumat.run(...);
+ * @endcode
+ *
+ * @see ElasticityModule, StrainMechanism, PlasticityMechanism,
+ * ViscoelasticMechanism, DamageMechanism
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace simcoon {
+
+// Forward declarations
+class PlasticityMechanism;
+class ViscoelasticMechanism;
+class DamageMechanism;
+enum class YieldType;
+enum class IsoHardType;
+enum class KinHardType;
+enum class DamageType;
+
+/**
+ * @brief Modular UMAT orchestrator
+ *
+ * This class provides a composable constitutive model framework where
+ * different elasticity types and strain mechanisms can be combined.
+ */
+class ModularUMAT {
+private:
+ // Modules
+ ElasticityModule elasticity_;
+ std::vector> mechanisms_;
+ InternalVariableCollection ivc_;
+
+ // State
+ double T_init_; ///< Initial temperature
+ arma::vec sigma_start_; ///< Stress at start of increment
+ bool initialized_; ///< Whether initialize() has been called
+
+ // Solver parameters
+ int maxiter_; ///< Maximum iterations for return mapping
+ double precision_; ///< Convergence tolerance
+
+public:
+ /**
+ * @brief Default constructor
+ */
+ ModularUMAT();
+
+ // Disable copy, allow move
+ ModularUMAT(const ModularUMAT&) = delete;
+ ModularUMAT& operator=(const ModularUMAT&) = delete;
+ ModularUMAT(ModularUMAT&&) = default;
+ ModularUMAT& operator=(ModularUMAT&&) = default;
+ ~ModularUMAT() = default;
+
+ // ========== Configuration ==========
+
+ /**
+ * @brief Set elasticity module
+ * @param type Elasticity type
+ * @param props Material properties
+ * @param offset Current offset in props (will be updated)
+ */
+ void set_elasticity(ElasticityType type, const arma::vec& props, int& offset);
+
+ /**
+ * @brief Add a plasticity mechanism
+ * @param yield_type Type of yield criterion
+ * @param iso_type Type of isotropic hardening
+ * @param kin_type Type of kinematic hardening
+ * @param N_iso Number of isotropic hardening terms
+ * @param N_kin Number of kinematic hardening terms
+ * @param props Material properties
+ * @param offset Current offset in props (will be updated)
+ * @return Reference to the added mechanism
+ */
+ PlasticityMechanism& add_plasticity(
+ YieldType yield_type,
+ IsoHardType iso_type,
+ KinHardType kin_type,
+ int N_iso,
+ int N_kin,
+ const arma::vec& props,
+ int& offset
+ );
+
+ /**
+ * @brief Add a viscoelastic mechanism (Prony series)
+ * @param N_prony Number of Prony terms
+ * @param props Material properties (g_i, tau_i pairs for each term)
+ * @param offset Current offset in props (will be updated)
+ * @return Reference to the added mechanism
+ */
+ ViscoelasticMechanism& add_viscoelasticity(
+ int N_prony,
+ const arma::vec& props,
+ int& offset
+ );
+
+ /**
+ * @brief Add a damage mechanism
+ * @param damage_type Type of damage evolution law
+ * @param props Material properties
+ * @param offset Current offset in props (will be updated)
+ * @return Reference to the added mechanism
+ */
+ DamageMechanism& add_damage(
+ DamageType damage_type,
+ const arma::vec& props,
+ int& offset
+ );
+
+ /**
+ * @brief Configure from props array
+ *
+ * Props format:
+ * - props[0]: elasticity_type (0=iso, 1=cubic, 2=trans_iso, 3=ortho)
+ * - props[1..N_el]: elasticity parameters
+ * - props[N_el+1]: num_mechanisms
+ * - For each mechanism:
+ * - props[i]: mechanism_type (0=plasticity, 1=viscoelasticity, 2=damage)
+ * - For plasticity (type 0):
+ * - yield_type, iso_type, kin_type, N_iso, N_kin, sigma_Y, [yield params], [iso params], [kin params]
+ * - For viscoelasticity (type 1):
+ * - N_prony, then N_prony pairs of (g_i, tau_i)
+ * - For damage (type 2):
+ * - damage_type (0=linear, 1=exp, 2=power, 3=weibull), Y_0, Y_c, [type-specific params]
+ *
+ * @param props Material properties vector
+ * @param offset Starting offset (default: 0)
+ */
+ void configure_from_props(const arma::vec& props, int offset = 0);
+
+ /**
+ * @brief Initialize internal variables
+ * @param nstatev Number of state variables
+ * @param statev State variable vector
+ */
+ void initialize(int nstatev, arma::vec& statev);
+
+ // ========== Accessors ==========
+
+ /**
+ * @brief Get the elasticity module
+ * @return Const reference to elasticity module
+ */
+ [[nodiscard]] const ElasticityModule& elasticity() const noexcept { return elasticity_; }
+
+ /**
+ * @brief Get number of mechanisms
+ * @return Number of strain mechanisms
+ */
+ [[nodiscard]] size_t num_mechanisms() const noexcept { return mechanisms_.size(); }
+
+ /**
+ * @brief Get mechanism by index
+ * @param i Index
+ * @return Reference to mechanism
+ */
+ StrainMechanism& mechanism(size_t i) { return *mechanisms_[i]; }
+ const StrainMechanism& mechanism(size_t i) const { return *mechanisms_[i]; }
+
+ /**
+ * @brief Get internal variable collection
+ * @return Reference to internal variable collection
+ */
+ InternalVariableCollection& internal_variables() { return ivc_; }
+ const InternalVariableCollection& internal_variables() const { return ivc_; }
+
+ /**
+ * @brief Get total number of state variables required
+ * @return nstatev
+ */
+ int required_nstatev() const;
+
+ /**
+ * @brief Check if initialized
+ * @return True if initialize() has been called
+ */
+ [[nodiscard]] bool is_initialized() const noexcept { return initialized_; }
+
+ // ========== Solver Parameters ==========
+
+ /**
+ * @brief Set maximum iterations
+ * @param maxiter Maximum number of return mapping iterations
+ */
+ void set_max_iterations(int maxiter) { maxiter_ = maxiter; }
+
+ /**
+ * @brief Set convergence precision
+ * @param precision Relative tolerance for convergence
+ */
+ void set_precision(double precision) { precision_ = precision; }
+
+ // ========== Main UMAT Entry Point ==========
+
+ /**
+ * @brief Run the constitutive model update
+ *
+ * This is the main entry point that performs:
+ * 1. Unpack state variables
+ * 2. Apply rotation for objectivity
+ * 3. Elastic prediction
+ * 4. Return mapping (if inelastic)
+ * 5. Compute consistent tangent
+ * 6. Compute work quantities
+ * 7. Pack state variables
+ *
+ * @param umat_name UMAT identifier
+ * @param Etot Total strain at end of increment
+ * @param DEtot Strain increment
+ * @param sigma Output: stress at end of increment
+ * @param Lt Output: consistent tangent modulus
+ * @param L Output: elastic stiffness
+ * @param DR Rotation increment matrix
+ * @param nprops Number of properties
+ * @param props Material properties
+ * @param nstatev Number of state variables
+ * @param statev State variables
+ * @param T Temperature at end of increment
+ * @param DT Temperature increment
+ * @param Time Current time
+ * @param DTime Time increment
+ * @param Wm Total mechanical work
+ * @param Wm_r Recoverable work
+ * @param Wm_ir Irrecoverable work
+ * @param Wm_d Dissipated work
+ * @param ndi Number of direct stress components
+ * @param nshr Number of shear stress components
+ * @param start True if first increment
+ * @param tnew_dt Suggested new time step ratio
+ */
+ void run(
+ const std::string& umat_name,
+ const arma::vec& Etot,
+ const arma::vec& DEtot,
+ arma::vec& sigma,
+ arma::mat& Lt,
+ arma::mat& L,
+ const arma::mat& DR,
+ int nprops,
+ const arma::vec& props,
+ int nstatev,
+ arma::vec& statev,
+ double T,
+ double DT,
+ double Time,
+ double DTime,
+ double& Wm,
+ double& Wm_r,
+ double& Wm_ir,
+ double& Wm_d,
+ int ndi,
+ int nshr,
+ bool start,
+ double& tnew_dt
+ );
+
+private:
+ /**
+ * @brief Perform return mapping algorithm
+ *
+ * Uses Newton iteration with Fischer-Burmeister complementarity
+ * to solve the coupled constraint equations from all mechanisms.
+ * The algorithm:
+ * 1. Computes elastic prediction (trial stress)
+ * 2. Evaluates constraint functions (Phi) from all mechanisms
+ * 3. Solves for multiplier increments via Fischer-Burmeister
+ * 4. Updates internal variables and recomputes stress
+ * 5. Repeats until convergence (error < precision_)
+ *
+ * @param Etot Total strain at start of increment
+ * @param DEtot Strain increment
+ * @param sigma Output: stress at end of increment
+ * @param T_init Reference temperature
+ * @param T Current temperature
+ * @param DT Temperature increment
+ * @param DTime Time increment
+ * @param ndi Number of direct stress components
+ * @param Ds_total Output: converged multiplier increments
+ *
+ * @see Fischer_Burmeister_m() in num_solve.hpp
+ */
+ void return_mapping(
+ const arma::vec& Etot,
+ const arma::vec& DEtot,
+ arma::vec& sigma,
+ double T_init,
+ double T,
+ double DT,
+ double DTime,
+ int ndi,
+ arma::vec& Ds_total
+ );
+
+ /**
+ * @brief Compute consistent tangent modulus
+ * @param sigma Current stress
+ * @param Ds_total Total multiplier increments
+ * @param Lt Output: tangent modulus
+ */
+ void compute_tangent(
+ const arma::vec& sigma,
+ const arma::vec& Ds_total,
+ arma::mat& Lt
+ );
+};
+
+/**
+ * @brief UMAT function for modular constitutive model
+ *
+ * Standard UMAT interface for the modular constitutive model.
+ * Can be registered in umat_smart.cpp for dispatch.
+ */
+void umat_modular(
+ const std::string& umat_name,
+ const arma::vec& Etot,
+ const arma::vec& DEtot,
+ arma::vec& sigma,
+ arma::mat& Lt,
+ arma::mat& L,
+ const arma::mat& DR,
+ const int& nprops,
+ const arma::vec& props,
+ const int& nstatev,
+ arma::vec& statev,
+ const double& T,
+ const double& DT,
+ const double& Time,
+ const double& DTime,
+ double& Wm,
+ double& Wm_r,
+ double& Wm_ir,
+ double& Wm_d,
+ const int& ndi,
+ const int& nshr,
+ const bool& start,
+ double& tnew_dt
+);
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/plasticity_mechanism.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/plasticity_mechanism.hpp
new file mode 100644
index 000000000..09cb9ff12
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/plasticity_mechanism.hpp
@@ -0,0 +1,175 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file plasticity_mechanism.hpp
+ * @brief Plasticity strain mechanism for modular UMAT.
+ *
+ * Combines yield criterion, isotropic hardening, and kinematic hardening
+ * into a complete plasticity model.
+ *
+ * @see StrainMechanism, YieldCriterion, IsotropicHardening, KinematicHardening, ModularUMAT
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include
+#include
+#include
+
+namespace simcoon {
+
+/**
+ * @brief Plasticity strain mechanism
+ *
+ * This class implements rate-independent plasticity by combining:
+ * - A yield criterion (von Mises, Hill, etc.)
+ * - Isotropic hardening (power-law, Voce, etc.)
+ * - Kinematic hardening (Prager, Chaboche, etc.)
+ *
+ * Internal variables:
+ * - "p": accumulated plastic strain (scalar)
+ * - "EP": plastic strain tensor (6 Voigt)
+ * - "X_i": backstress tensors from kinematic hardening
+ */
+class PlasticityMechanism final : public StrainMechanism {
+private:
+ YieldType yield_type_; ///< Yield criterion type (for deferred configuration)
+ std::unique_ptr yield_;
+ std::unique_ptr iso_hard_;
+ std::unique_ptr kin_hard_;
+ double sigma_Y_; ///< Initial yield stress
+
+ // Cached quantities from last constraint computation
+ mutable arma::vec flow_dir_; ///< Flow direction
+ mutable arma::vec kappa_; ///< L * flow_direction
+ mutable double dPhi_dp_; ///< dPhi/dp
+ mutable double H_total_; ///< Total hardening modulus
+
+public:
+ /**
+ * @brief Constructor
+ * @param yield_type Type of yield criterion
+ * @param iso_type Type of isotropic hardening
+ * @param kin_type Type of kinematic hardening
+ * @param N_iso Number of isotropic hardening terms (for COMBINED_VOCE)
+ * @param N_kin Number of kinematic hardening terms (for CHABOCHE)
+ */
+ PlasticityMechanism(
+ YieldType yield_type = YieldType::VON_MISES,
+ IsoHardType iso_type = IsoHardType::POWER_LAW,
+ KinHardType kin_type = KinHardType::NONE,
+ int N_iso = 1,
+ int N_kin = 1
+ );
+
+ // Default move, no copy (due to unique_ptr)
+ PlasticityMechanism(const PlasticityMechanism&) = delete;
+ PlasticityMechanism& operator=(const PlasticityMechanism&) = delete;
+ PlasticityMechanism(PlasticityMechanism&&) = default;
+ PlasticityMechanism& operator=(PlasticityMechanism&&) = default;
+ ~PlasticityMechanism() override = default;
+
+ // ========== Configuration ==========
+
+ void configure(const arma::vec& props, int& offset) override;
+ void register_variables(InternalVariableCollection& ivc) override;
+
+ // ========== Mechanism Properties ==========
+
+ [[nodiscard]] MechanismType type() const override { return MechanismType::PLASTICITY; }
+ [[nodiscard]] int num_constraints() const override { return 1; } // Single yield surface
+
+ /**
+ * @brief Get the yield criterion
+ * @return Reference to yield criterion
+ */
+ [[nodiscard]] const YieldCriterion& yield_criterion() const noexcept { return *yield_; }
+
+ /**
+ * @brief Get the isotropic hardening
+ * @return Reference to isotropic hardening
+ */
+ [[nodiscard]] const IsotropicHardening& isotropic_hardening() const noexcept { return *iso_hard_; }
+
+ /**
+ * @brief Get the kinematic hardening
+ * @return Reference to kinematic hardening
+ */
+ [[nodiscard]] const KinematicHardening& kinematic_hardening() const noexcept { return *kin_hard_; }
+
+ /**
+ * @brief Get initial yield stress
+ * @return sigma_Y
+ */
+ [[nodiscard]] double initial_yield_stress() const noexcept { return sigma_Y_; }
+
+ // ========== Constitutive Computations ==========
+
+ void compute_constraints(
+ const arma::vec& sigma,
+ const arma::mat& L,
+ double DTime,
+ const InternalVariableCollection& ivc,
+ arma::vec& Phi,
+ arma::vec& Y_crit
+ ) const override;
+
+ void compute_flow_directions(
+ const arma::vec& sigma,
+ const InternalVariableCollection& ivc,
+ std::map& Lambda_map
+ ) const override;
+
+ void compute_jacobian_contribution(
+ const arma::vec& sigma,
+ const arma::mat& L,
+ const InternalVariableCollection& ivc,
+ arma::mat& B,
+ int row_offset
+ ) const override;
+
+ arma::vec inelastic_strain(const InternalVariableCollection& ivc) const override;
+
+ void update(
+ const arma::vec& ds,
+ int offset,
+ InternalVariableCollection& ivc
+ ) override;
+
+ void tangent_contribution(
+ const arma::vec& sigma,
+ const arma::mat& L,
+ const arma::vec& Ds,
+ int offset,
+ const InternalVariableCollection& ivc,
+ arma::mat& Lt
+ ) const override;
+
+ void compute_work(
+ const arma::vec& sigma_start,
+ const arma::vec& sigma,
+ const InternalVariableCollection& ivc,
+ double& Wm_r,
+ double& Wm_ir,
+ double& Wm_d
+ ) const override;
+};
+
+} // namespace simcoon
diff --git a/include/simcoon/Continuum_mechanics/Umat/Modular/strain_mechanism.hpp b/include/simcoon/Continuum_mechanics/Umat/Modular/strain_mechanism.hpp
new file mode 100644
index 000000000..2838d60c4
--- /dev/null
+++ b/include/simcoon/Continuum_mechanics/Umat/Modular/strain_mechanism.hpp
@@ -0,0 +1,228 @@
+/* This file is part of simcoon.
+
+simcoon is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+simcoon is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with simcoon. If not, see .
+
+*/
+
+/**
+ * @file strain_mechanism.hpp
+ * @brief Base class for strain mechanisms in modular UMAT.
+ *
+ * Strain mechanisms are pluggable components that contribute to the
+ * total strain decomposition: Etot = Eel + E_plastic + E_viscous + ...
+ *
+ * Each mechanism:
+ * - Registers its internal variables
+ * - Contributes constraint functions (yield/evolution equations)
+ * - Contributes to the Jacobian matrix
+ * - Updates its internal variables during return mapping
+ *
+ * @see PlasticityMechanism, ViscoelasticMechanism, DamageMechanism, ModularUMAT
+ * @version 1.0
+ */
+
+#pragma once
+
+#include
+#include