From 19da18a35039d041415a642b332b77a14beaffd8 Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 6 Mar 2025 14:50:24 -0600 Subject: [PATCH 01/18] The updated adhesion with the development branch --- src/Model/Solid/Finite/Adhesion.H | 225 ++++++++++++++++++++++++++++++ src/mechanics.cc | 2 + src/test.cc | 4 +- 3 files changed, 230 insertions(+), 1 deletion(-) create mode 100644 src/Model/Solid/Finite/Adhesion.H diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H new file mode 100644 index 0000000000..2a16f7b4d7 --- /dev/null +++ b/src/Model/Solid/Finite/Adhesion.H @@ -0,0 +1,225 @@ + +#ifndef MODEL_SOLID_FINITE_ADHESION_H_ +#define MODEL_SOLID_FINITE_ADHESION_H_ + +#include "IO/ParmParse.H" +#include "Model/Solid/Solid.H" +#include "Set/Set.H" +#include // For Eigen::Matrix3d +#include // For std::array to represent fourth-order tensors +#include // For std::pow and std::exp + +namespace Set { + +// Outer product for two 3x3 matrices producing a 3x3x3x3 tensor +inline std::array, 3>, 3>, 3> +Outer(const Eigen::Matrix3d &A, const Eigen::Matrix3d &B) { +std::array, 3>, 3>, 3> result = {}; +for (int i = 0; i < 3; ++i) +for (int j = 0; j < 3; ++j) +for (int k = 0; k < 3; ++k) +for (int l = 0; l < 3; ++l) +result[i][j][k][l] = A(i, j) * B(k, l); +return result; +} + +// Derivative of inverse transpose: d(F^-T) / dF +inline std::array, 3>, 3>, 3> +Derivative(const Eigen::Matrix3d &FinvT) { +std::array, 3>, 3>, 3> result = {}; +for (int i = 0; i < 3; ++i) +for (int j = 0; j < 3; ++j) +for (int k = 0; k < 3; ++k) +for (int l = 0; l < 3; ++l) +result[i][j][k][l] = -FinvT(i, j) * FinvT(j, i); +return result; +} + +} // namespace Set + +namespace Model { +namespace Solid { +namespace Finite { + +class Adhesion : public Solid { +public: +Adhesion() {} +Adhesion(Solid base) : Solid(base) {} +virtual ~Adhesion() {} + +// Material parameters (declared in order) +Set::Scalar d = NAN; ///< Adhesion strength +Set::Scalar mu = NAN; ///< Shear modulus +Set::Scalar kappa = NAN; ///< Bulk modulus +Set::Scalar zeta = NAN; ///< Adhesion decay parameter +Set::Scalar n = NAN; ///< Decay exponent + +static Adhesion Zero() { +Adhesion adhesion; +adhesion.d = 0; +adhesion.mu = 0; +adhesion.kappa = 0; +adhesion.zeta = 0; +adhesion.n = 0; +return adhesion; +} + +static Adhesion Random() { +Adhesion adhesion; +adhesion.d = Util::Random(); +adhesion.mu = Util::Random(); +adhesion.kappa = Util::Random(); +adhesion.zeta = Util::Random(); +adhesion.n = Util::Random(); +return adhesion; +} + +// Updated Parse function: +// If "E" and "nu" (and optionally "F0") are provided, use these to compute mu and kappa. +// Otherwise, query for "mu" and "kappa" directly. +static void Parse(Adhesion &value, IO::ParmParse &pp) { +if (pp.contains("E") && pp.contains("nu")) { +Set::Scalar E, nu; +pp_query("E", E); // Elastic modulus +pp_query("nu", nu); // Poisson's ratio +// Map E and nu to shear and bulk moduli (using a convention similar to NeoHookean citeturn0file2) +value.mu = E / (2.0 * (1.0 + nu)); +value.kappa = E / (3.0 * (1.0 - 2.0 * nu)); +// Consume F0 if provided, even though Adhesion does not use it. +if (pp.contains("F0")) { +Set::Matrix dummy; +pp_queryarr("F0", dummy); // F0 +} +} else { +pp_query("mu", value.mu); // Shear modulus +pp_query("kappa", value.kappa); // Bulk modulus +} +// Query the remaining Adhesion-specific parameters. +pp_query("d", value.d); // Adhesion strength +pp_query("zeta", value.zeta); // Adhesion decay parameter +pp_query("n", value.n); // Decay exponent +} + +Set::Scalar W(const Set::Matrix &a_F) const override { +#if AMREX_SPACEDIM == 2 +Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); +F(0,0) = a_F(0,0); +F(0,1) = a_F(0,1); +F(1,0) = a_F(1,0); +F(1,1) = a_F(1,1); +#elif AMREX_SPACEDIM == 3 +Eigen::Matrix3d F = a_F; +#endif +Set::Scalar J = F.determinant(); +Set::Scalar J23 = std::pow(fabs(J), 2.0/3.0); +Set::Scalar w = 0.0; +// Shear and bulk contributions (scaled by adhesion factor d) +w += d * (0.5 * mu * (((F.transpose() * F).trace()) - 3) / J23 + +0.5 * kappa * (J - 1.0) * (J - 1.0)); +// Adhesion decay term +w += zeta / std::pow(J, n); +return w; +} + +Set::Matrix DW(const Set::Matrix &a_F) const override { +#if AMREX_SPACEDIM == 2 +Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); +F(0,0) = a_F(0,0); +F(0,1) = a_F(0,1); +F(1,0) = a_F(1,0); +F(1,1) = a_F(1,1); +#elif AMREX_SPACEDIM == 3 +Eigen::Matrix3d F = a_F; +#endif +Set::Scalar J = F.determinant(); +Set::Scalar J23 = std::pow(fabs(J), 2.0/3.0); +Eigen::Matrix3d FinvT = F.inverse().transpose(); + +Eigen::Matrix3d dw = d * mu * (F / J23 - +(1.0/3.0) * (((F.transpose() * F).trace()) - 3) * FinvT / J23) ++ d * kappa * (J - 1.0) * J * FinvT +- zeta * n * std::pow(J, -n) * FinvT; +#if AMREX_SPACEDIM == 2 +Set::Matrix r_dw; +r_dw(0,0) = dw(0,0); +r_dw(0,1) = dw(0,1); +r_dw(1,0) = dw(1,0); +r_dw(1,1) = dw(1,1); +return r_dw; +#elif AMREX_SPACEDIM == 3 +return dw; +#endif +} + +Set::Matrix4 DDW(const Set::Matrix &a_F) const override { +#if AMREX_SPACEDIM == 2 +Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); +F(0,0) = a_F(0,0); +F(0,1) = a_F(0,1); +F(1,0) = a_F(1,0); +F(1,1) = a_F(1,1); +#elif AMREX_SPACEDIM == 3 +Eigen::Matrix3d F = a_F; +#endif +Set::Scalar J = F.determinant(); +Set::Scalar J23 = std::pow(fabs(J), 2.0/3.0); +Eigen::Matrix3d FinvT = F.inverse().transpose(); +Set::Scalar trace_FTF = (F.transpose() * F).trace(); + +auto Deriv = ::Set::Derivative(FinvT); + +Set::Matrix4<3, Set::Sym::Major> ddw; +for (int i = 0; i < 3; ++i) { +for (int j = 0; j < 3; ++j) { +for (int k = 0; k < 3; ++k) { +for (int l = 0; l < 3; ++l) { +ddw(i, j, k, l) = 0.0; +Set::Scalar t1 = 0.0, t2 = 0.0, t3 = 0.0; +// Shear-related term (mu) +if (i == k && j == l) +t1 += 1.0; +t1 -= (2.0/3.0)*F(i,j)*FinvT(k,l); +t1 -= (2.0/3.0)*FinvT(i,j)*F(k,l); +t1 += (2.0/9.0)*(trace_FTF - 3)*FinvT(i,j)*FinvT(k,l); +t1 += (1.0/3.0)*(trace_FTF - 3)*FinvT(i,l)*FinvT(k,j); +// Bulk-related term (kappa) +t2 += (2.0*J*J - J)*FinvT(i,j)*FinvT(k,l) ++ (J - J*J)*FinvT(i,l)*FinvT(k,j); +// Adhesion-related term (zeta) +t3 += (n*n)*std::pow(J, -n)*FinvT(i,j)*FinvT(k,l) ++ n*std::pow(J, -n)*FinvT(i,l)*FinvT(k,j); +ddw(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; +} +} +} +} +#if AMREX_SPACEDIM == 2 +Set::Matrix4 r_ddw; +for (int i = 0; i < 2; i++) +for (int j = 0; j < 2; j++) +for (int k = 0; k < 2; k++) +for (int l = 0; l < 2; l++) +r_ddw(i, j, k, l) = ddw(i, j, k, l); +return r_ddw; +#elif AMREX_SPACEDIM == 3 +return ddw; +#endif +} + +virtual void Print(std::ostream &out) const override { +out << "d = " << d << ", mu = " << mu << ", kappa = " << kappa +<< ", zeta = " << zeta << ", n = " << n; +} + +#define OP_CLASS Adhesion +#define OP_VARS X(d) X(mu) X(kappa) X(zeta) X(n) +#include "Model/Solid/InClassOperators.H" +}; + +#include "Model/Solid/ExtClassOperators.H" + +} // namespace Finite +} // namespace Solid +} // namespace Model + +#endif diff --git a/src/mechanics.cc b/src/mechanics.cc index b302d3a3d6..d1a2dab401 100644 --- a/src/mechanics.cc +++ b/src/mechanics.cc @@ -12,6 +12,7 @@ #include "Model/Solid/Affine/Isotropic.H" #include "Model/Solid/Linear/Cubic.H" #include "Model/Solid/Affine/Cubic.H" +#include "Model/Solid/Finite/Adhesion.H" #include "Model/Solid/Finite/NeoHookean.H" #include "Model/Solid/Finite/NeoHookeanPredeformed.H" #include "Model/Solid/Finite/PseudoLinear/Cubic.H" @@ -47,6 +48,7 @@ int main (int argc, char* argv[]) else if (model == "finite.neohookeanpre") integrator = new Integrator::Mechanics(pp); else if (model == "finite.pseudolinear.cubic") integrator = new Integrator::Mechanics(pp); else if (model == "finite.pseudoaffine.cubic") integrator = new Integrator::Mechanics(pp); + else if (model == "finite.adhesion") integrator = new Integrator::Mechanics(pp); else if (model == "affine.j2") integrator = new Integrator::Mechanics(pp); else if (model == "finite.crystalplastic") integrator = new Integrator::Mechanics(pp); else Util::Abort(INFO,model," is not a valid model"); diff --git a/src/test.cc b/src/test.cc index c5deb8f3df..ca05cf3902 100644 --- a/src/test.cc +++ b/src/test.cc @@ -15,6 +15,7 @@ #include "Model/Solid/Linear/Laplacian.H" #include "Model/Solid/Affine/Isotropic.H" #include "Model/Solid/Affine/Cubic.H" +#include "Model/Solid/Finite/Adhesion.H" #include "Model/Solid/Finite/NeoHookean.H" #include "Model/Solid/Finite/NeoHookeanPredeformed.H" #include "Model/Solid/Finite/PseudoLinear/Cubic.H" @@ -57,7 +58,8 @@ int main (int argc, char* argv[]) MODELTEST(Model::Solid::Finite::PseudoLinear::Cubic); MODELTEST(Model::Solid::Finite::NeoHookeanPredeformed); MODELTEST(Model::Solid::Finite::PseudoAffine::Cubic); - + MODELTEST(Model::Solid::Finite::Adhesion); + Util::Test::Message("Set::Matrix4"); { From 0bbc22b4b794a5100bff1945b0d6f5584bc3994f Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 27 Mar 2025 14:14:05 -0500 Subject: [PATCH 02/18] Merging development into adhesion modifying the style error --- src/mechanics.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mechanics.cc b/src/mechanics.cc index 05fd9cd355..1c83b440f9 100644 --- a/src/mechanics.cc +++ b/src/mechanics.cc @@ -57,7 +57,7 @@ int main (int argc, char* argv[]) pp.select_only>(integrator); else if (model == "finite.neohookeanpre") pp.select_only>(integrator); - else if (model == "finite.adhesion") + else if (model == "finite.adhesion") pp.select_only>(integrator); else if (model == "finite.pseudolinear.cubic") pp.select_only>(integrator); From 4021170fea1d65b6d8ae79ba8bdb304bbaa5999b Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 27 Mar 2025 15:07:54 -0500 Subject: [PATCH 03/18] Merging development into adhesion modifying all errors --- src/mechanics.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mechanics.cc b/src/mechanics.cc index 1c83b440f9..8542b8b262 100644 --- a/src/mechanics.cc +++ b/src/mechanics.cc @@ -57,7 +57,7 @@ int main (int argc, char* argv[]) pp.select_only>(integrator); else if (model == "finite.neohookeanpre") pp.select_only>(integrator); - else if (model == "finite.adhesion") + else if (model == "finite.adhesion") pp.select_only>(integrator); else if (model == "finite.pseudolinear.cubic") pp.select_only>(integrator); From d6718e7556d56d50c978dc627bc2dfb282456647 Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Tue, 1 Apr 2025 10:26:45 -0500 Subject: [PATCH 04/18] Modifying clang and documentation and adding regression tests --- scripts/runtests.py | 2 +- src/Model/Solid/Finite/Adhesion.H | 30 ++++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/scripts/runtests.py b/scripts/runtests.py index d28c12573d..9272c0a4be 100755 --- a/scripts/runtests.py +++ b/scripts/runtests.py @@ -629,4 +629,4 @@ class stats: return_code += stats.timeouts # Return nonzero only if no tests failed or were unexpectedly skipped -exit(return_code) +exit(return_code) \ No newline at end of file diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 2a16f7b4d7..1b81ae5f91 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -1,3 +1,33 @@ +// +// This models a finite-strain, solid with adhesion behavior governed by +// a combination of elastic energy and an adhesion decay term. +// +// The energy and derivatives are: +// +// .. math:: +// :nowrap: +// +// \begin{gather} +// W = d \left[ \frac{\mu}{2J^{2/3}}(F^T F : I - 3) + \frac{\kappa}{2}(J - 1)^2 \right] + \frac{\zeta}{J^n} \\ +// DW = \text{(first Piola-Kirchhoff stress)} \\ +// DDW = \text{(tangent stiffness tensor)} +// \end{gather} +// +// where: +// - :math:`d` is the adhesion strength, +// - :math:`\mu` is the shear modulus, +// - :math:`\kappa` is the bulk modulus, +// - :math:`\zeta` is the adhesion decay coefficient, +// - :math:`n` is the decay exponent, +// - :math:`F` is the deformation gradient, +// - :math:`J = \det(F)` is the volume ratio, +// - :math:`F^T F : I = \text{trace}(F^T F)` is used to capture isochoric deformation. +// +// Notes: +// +// * The model does not include rate-dependence or plasticity directly, +// and is suitable for modeling reversible adhesion and separation under deformation. +// #ifndef MODEL_SOLID_FINITE_ADHESION_H_ #define MODEL_SOLID_FINITE_ADHESION_H_ From 237c76d79cf445541361017eabd609801ae3056b Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Tue, 1 Apr 2025 12:52:49 -0500 Subject: [PATCH 05/18] Modifying in clang format and documentation and adding regression tests and amall modofication in the adhesion --- src/Model/Solid/Finite/Adhesion.H | 380 ++++++++++++++++-------------- 1 file changed, 205 insertions(+), 175 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 1b81ae5f91..1357767fd1 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -1,5 +1,5 @@ // -// This models a finite-strain, solid with adhesion behavior governed by +// This models a finite-strain, solid with adhesion behavior governed by // a combination of elastic energy and an adhesion decay term. // // The energy and derivatives are: @@ -14,10 +14,10 @@ // \end{gather} // // where: -// - :math:`d` is the adhesion strength, +// - :math:`d` The coefficient of degradable part, // - :math:`\mu` is the shear modulus, // - :math:`\kappa` is the bulk modulus, -// - :math:`\zeta` is the adhesion decay coefficient, +// - :math:`\zeta` is the adhesion decay coefficient or the coefficient of non-degradable part, // - :math:`n` is the decay exponent, // - :math:`F` is the deformation gradient, // - :math:`J = \det(F)` is the volume ratio, @@ -25,8 +25,11 @@ // // Notes: // -// * The model does not include rate-dependence or plasticity directly, -// and is suitable for modeling reversible adhesion and separation under deformation. +// Adhesion energy model W(F): +// Consists of two parts: +// 1. Degradable (elastic shear + volumetric terms) +// 2. Non-degradable (adhesion term decaying with volume expansion) +// This model is rate-independent, reversible, and purely elastic. // #ifndef MODEL_SOLID_FINITE_ADHESION_H_ @@ -35,211 +38,238 @@ #include "IO/ParmParse.H" #include "Model/Solid/Solid.H" #include "Set/Set.H" -#include // For Eigen::Matrix3d -#include // For std::array to represent fourth-order tensors -#include // For std::pow and std::exp +#include // For std::array to represent fourth-order tensors +#include // For std::pow and std::exp +#include // For Eigen::Matrix3d -namespace Set { +namespace Set +{ // Outer product for two 3x3 matrices producing a 3x3x3x3 tensor inline std::array, 3>, 3>, 3> -Outer(const Eigen::Matrix3d &A, const Eigen::Matrix3d &B) { -std::array, 3>, 3>, 3> result = {}; -for (int i = 0; i < 3; ++i) -for (int j = 0; j < 3; ++j) -for (int k = 0; k < 3; ++k) -for (int l = 0; l < 3; ++l) -result[i][j][k][l] = A(i, j) * B(k, l); -return result; +Outer(const Eigen::Matrix3d &A, const Eigen::Matrix3d &B) +{ + std::array, 3>, 3>, 3> result = {}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + for (int k = 0; k < 3; ++k) + for (int l = 0; l < 3; ++l) + result[i][j][k][l] = A(i, j) * B(k, l); + return result; } // Derivative of inverse transpose: d(F^-T) / dF inline std::array, 3>, 3>, 3> -Derivative(const Eigen::Matrix3d &FinvT) { -std::array, 3>, 3>, 3> result = {}; -for (int i = 0; i < 3; ++i) -for (int j = 0; j < 3; ++j) -for (int k = 0; k < 3; ++k) -for (int l = 0; l < 3; ++l) -result[i][j][k][l] = -FinvT(i, j) * FinvT(j, i); -return result; +Derivative(const Eigen::Matrix3d &FinvT) +{ + std::array, 3>, 3>, 3> result = {}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + for (int k = 0; k < 3; ++k) + for (int l = 0; l < 3; ++l) + result[i][j][k][l] = -FinvT(i, j) * FinvT(j, i); + return result; } } // namespace Set -namespace Model { -namespace Solid { -namespace Finite { +namespace Model +{ +namespace Solid +{ +namespace Finite +{ -class Adhesion : public Solid { +class Adhesion : public Solid +{ public: -Adhesion() {} -Adhesion(Solid base) : Solid(base) {} -virtual ~Adhesion() {} - -// Material parameters (declared in order) -Set::Scalar d = NAN; ///< Adhesion strength -Set::Scalar mu = NAN; ///< Shear modulus -Set::Scalar kappa = NAN; ///< Bulk modulus -Set::Scalar zeta = NAN; ///< Adhesion decay parameter -Set::Scalar n = NAN; ///< Decay exponent - -static Adhesion Zero() { -Adhesion adhesion; -adhesion.d = 0; -adhesion.mu = 0; -adhesion.kappa = 0; -adhesion.zeta = 0; -adhesion.n = 0; -return adhesion; -} + Adhesion() {} + Adhesion(Solid base) : Solid(base) {} + virtual ~Adhesion() {} -static Adhesion Random() { -Adhesion adhesion; -adhesion.d = Util::Random(); -adhesion.mu = Util::Random(); -adhesion.kappa = Util::Random(); -adhesion.zeta = Util::Random(); -adhesion.n = Util::Random(); -return adhesion; -} + // Material parameters (declared in order) + Set::Scalar d = NAN; ///< Adhesion strength + Set::Scalar mu = NAN; ///< Shear modulus + Set::Scalar kappa = NAN; ///< Bulk modulus + Set::Scalar zeta = NAN; ///< Adhesion decay parameter + Set::Scalar n = NAN; ///< Decay exponent -// Updated Parse function: -// If "E" and "nu" (and optionally "F0") are provided, use these to compute mu and kappa. -// Otherwise, query for "mu" and "kappa" directly. -static void Parse(Adhesion &value, IO::ParmParse &pp) { -if (pp.contains("E") && pp.contains("nu")) { -Set::Scalar E, nu; -pp_query("E", E); // Elastic modulus -pp_query("nu", nu); // Poisson's ratio -// Map E and nu to shear and bulk moduli (using a convention similar to NeoHookean citeturn0file2) -value.mu = E / (2.0 * (1.0 + nu)); -value.kappa = E / (3.0 * (1.0 - 2.0 * nu)); -// Consume F0 if provided, even though Adhesion does not use it. -if (pp.contains("F0")) { -Set::Matrix dummy; -pp_queryarr("F0", dummy); // F0 -} -} else { -pp_query("mu", value.mu); // Shear modulus -pp_query("kappa", value.kappa); // Bulk modulus -} -// Query the remaining Adhesion-specific parameters. -pp_query("d", value.d); // Adhesion strength -pp_query("zeta", value.zeta); // Adhesion decay parameter -pp_query("n", value.n); // Decay exponent -} + static Adhesion + Zero() + { + Adhesion adhesion; + adhesion.d = 0; + adhesion.mu = 0; + adhesion.kappa = 0; + adhesion.zeta = 0; + adhesion.n = 0; + return adhesion; + } -Set::Scalar W(const Set::Matrix &a_F) const override { + static Adhesion + Random() + { + Adhesion adhesion; + adhesion.d = Util::Random(); + adhesion.mu = Util::Random(); + adhesion.kappa = Util::Random(); + adhesion.zeta = Util::Random(); + adhesion.n = Util::Random(); + return adhesion; + } + + // Updated Parse function: + // If "E" and "nu" (and optionally "F0") are provided, use these to compute mu and kappa. + // Otherwise, query for "mu" and "kappa" directly. + static void + Parse(Adhesion &value, IO::ParmParse &pp) + { + if (pp.contains("E") && pp.contains("nu")) + { + Set::Scalar E, nu; + pp_query("E", E); // Elastic modulus + pp_query("nu", nu); // Poisson's ratio + // Map E and nu to shear and bulk moduli (using a convention similar to NeoHookean citeturn0file2) + value.mu = E / (2.0 * (1.0 + nu)); + value.kappa = E / (3.0 * (1.0 - 2.0 * nu)); + // Consume F0 if provided, even though Adhesion does not use it. + if (pp.contains("F0")) + { + Set::Matrix dummy; + pp_queryarr("F0", dummy); // F0 + } + } + else + { + pp_query("mu", value.mu); // Shear modulus + pp_query("kappa", value.kappa); // Bulk modulus + } + // Query the remaining Adhesion-specific parameters. + pp_query("d", value.d); // Adhesion strength + pp_query("zeta", value.zeta); // Adhesion decay parameter + pp_query("n", value.n); // Decay exponent + } + + Set::Scalar + W(const Set::Matrix &a_F) const override + { #if AMREX_SPACEDIM == 2 -Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); -F(0,0) = a_F(0,0); -F(0,1) = a_F(0,1); -F(1,0) = a_F(1,0); -F(1,1) = a_F(1,1); + Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); + F(0, 0) = a_F(0, 0); + F(0, 1) = a_F(0, 1); + F(1, 0) = a_F(1, 0); + F(1, 1) = a_F(1, 1); #elif AMREX_SPACEDIM == 3 -Eigen::Matrix3d F = a_F; + Eigen::Matrix3d F = a_F; #endif -Set::Scalar J = F.determinant(); -Set::Scalar J23 = std::pow(fabs(J), 2.0/3.0); -Set::Scalar w = 0.0; -// Shear and bulk contributions (scaled by adhesion factor d) -w += d * (0.5 * mu * (((F.transpose() * F).trace()) - 3) / J23 + -0.5 * kappa * (J - 1.0) * (J - 1.0)); -// Adhesion decay term -w += zeta / std::pow(J, n); -return w; -} + Set::Scalar J = F.determinant(); + Set::Scalar J23 = std::pow(fabs(J), 2.0 / 3.0); + Set::Scalar w = 0.0; + // Shear and bulk contributions (scaled by adhesion factor d) + w += d * (0.5 * mu * (((F.transpose() * F).trace()) - 3) / J23 + 0.5 * kappa * (J - 1.0) * (J - 1.0)); + // Adhesion decay term + w += zeta / std::pow(J, n); + return w; + } -Set::Matrix DW(const Set::Matrix &a_F) const override { + Set::Matrix + DW(const Set::Matrix &a_F) const override + { #if AMREX_SPACEDIM == 2 -Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); -F(0,0) = a_F(0,0); -F(0,1) = a_F(0,1); -F(1,0) = a_F(1,0); -F(1,1) = a_F(1,1); + Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); + F(0, 0) = a_F(0, 0); + F(0, 1) = a_F(0, 1); + F(1, 0) = a_F(1, 0); + F(1, 1) = a_F(1, 1); #elif AMREX_SPACEDIM == 3 -Eigen::Matrix3d F = a_F; + Eigen::Matrix3d F = a_F; #endif -Set::Scalar J = F.determinant(); -Set::Scalar J23 = std::pow(fabs(J), 2.0/3.0); -Eigen::Matrix3d FinvT = F.inverse().transpose(); - -Eigen::Matrix3d dw = d * mu * (F / J23 - -(1.0/3.0) * (((F.transpose() * F).trace()) - 3) * FinvT / J23) -+ d * kappa * (J - 1.0) * J * FinvT -- zeta * n * std::pow(J, -n) * FinvT; + Set::Scalar J = F.determinant(); + Set::Scalar J23 = std::pow(fabs(J), 2.0 / 3.0); + Eigen::Matrix3d FinvT = F.inverse().transpose(); + + Eigen::Matrix3d dw = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3) * FinvT / J23) + + d * kappa * (J - 1.0) * J * FinvT + - zeta * n * std::pow(J, -n) * FinvT; #if AMREX_SPACEDIM == 2 -Set::Matrix r_dw; -r_dw(0,0) = dw(0,0); -r_dw(0,1) = dw(0,1); -r_dw(1,0) = dw(1,0); -r_dw(1,1) = dw(1,1); -return r_dw; + Set::Matrix r_dw; + r_dw(0, 0) = dw(0, 0); + r_dw(0, 1) = dw(0, 1); + r_dw(1, 0) = dw(1, 0); + r_dw(1, 1) = dw(1, 1); + return r_dw; #elif AMREX_SPACEDIM == 3 -return dw; + return dw; #endif -} + } -Set::Matrix4 DDW(const Set::Matrix &a_F) const override { + Set::Matrix4 + DDW(const Set::Matrix &a_F) const override + { #if AMREX_SPACEDIM == 2 -Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); -F(0,0) = a_F(0,0); -F(0,1) = a_F(0,1); -F(1,0) = a_F(1,0); -F(1,1) = a_F(1,1); + Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); + F(0, 0) = a_F(0, 0); + F(0, 1) = a_F(0, 1); + F(1, 0) = a_F(1, 0); + F(1, 1) = a_F(1, 1); #elif AMREX_SPACEDIM == 3 -Eigen::Matrix3d F = a_F; + Eigen::Matrix3d F = a_F; #endif -Set::Scalar J = F.determinant(); -Set::Scalar J23 = std::pow(fabs(J), 2.0/3.0); -Eigen::Matrix3d FinvT = F.inverse().transpose(); -Set::Scalar trace_FTF = (F.transpose() * F).trace(); - -auto Deriv = ::Set::Derivative(FinvT); - -Set::Matrix4<3, Set::Sym::Major> ddw; -for (int i = 0; i < 3; ++i) { -for (int j = 0; j < 3; ++j) { -for (int k = 0; k < 3; ++k) { -for (int l = 0; l < 3; ++l) { -ddw(i, j, k, l) = 0.0; -Set::Scalar t1 = 0.0, t2 = 0.0, t3 = 0.0; -// Shear-related term (mu) -if (i == k && j == l) -t1 += 1.0; -t1 -= (2.0/3.0)*F(i,j)*FinvT(k,l); -t1 -= (2.0/3.0)*FinvT(i,j)*F(k,l); -t1 += (2.0/9.0)*(trace_FTF - 3)*FinvT(i,j)*FinvT(k,l); -t1 += (1.0/3.0)*(trace_FTF - 3)*FinvT(i,l)*FinvT(k,j); -// Bulk-related term (kappa) -t2 += (2.0*J*J - J)*FinvT(i,j)*FinvT(k,l) -+ (J - J*J)*FinvT(i,l)*FinvT(k,j); -// Adhesion-related term (zeta) -t3 += (n*n)*std::pow(J, -n)*FinvT(i,j)*FinvT(k,l) -+ n*std::pow(J, -n)*FinvT(i,l)*FinvT(k,j); -ddw(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; -} -} -} -} + Set::Scalar J = F.determinant(); + Set::Scalar J23 = std::pow(fabs(J), 2.0 / 3.0); + Eigen::Matrix3d FinvT = F.inverse().transpose(); + Set::Scalar trace_FTF = (F.transpose() * F).trace(); + + auto Deriv = ::Set::Derivative(FinvT); + + Set::Matrix4<3, Set::Sym::Major> ddw; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + for (int k = 0; k < 3; ++k) + { + for (int l = 0; l < 3; ++l) + { + ddw(i, j, k, l) = 0.0; + Set::Scalar t1 = 0.0, t2 = 0.0, t3 = 0.0; + // Shear-related term (mu) + if (i == k && j == l) + t1 += 1.0; + t1 -= (2.0 / 3.0) * F(i, j) * FinvT(k, l); + t1 -= (2.0 / 3.0) * FinvT(i, j) * F(k, l); + t1 += (2.0 / 9.0) * (trace_FTF - 3) * FinvT(i, j) * FinvT(k, l); + t1 += (1.0 / 3.0) * (trace_FTF - 3) * FinvT(i, l) * FinvT(k, j); + // Bulk-related term (kappa) + t2 += (2.0 * J * J - J) * FinvT(i, j) * FinvT(k, l) + + (J - J * J) * FinvT(i, l) * FinvT(k, j); + // Adhesion-related term (zeta) + t3 += (n * n) * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) + + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); + ddw(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; + } + } + } + } #if AMREX_SPACEDIM == 2 -Set::Matrix4 r_ddw; -for (int i = 0; i < 2; i++) -for (int j = 0; j < 2; j++) -for (int k = 0; k < 2; k++) -for (int l = 0; l < 2; l++) -r_ddw(i, j, k, l) = ddw(i, j, k, l); -return r_ddw; + Set::Matrix4 r_ddw; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + for (int k = 0; k < 2; k++) + for (int l = 0; l < 2; l++) + r_ddw(i, j, k, l) = ddw(i, j, k, l); + return r_ddw; #elif AMREX_SPACEDIM == 3 -return ddw; + return ddw; #endif -} + } -virtual void Print(std::ostream &out) const override { -out << "d = " << d << ", mu = " << mu << ", kappa = " << kappa -<< ", zeta = " << zeta << ", n = " << n; -} + virtual void + Print(std::ostream &out) const override + { + out << "d = " << d << ", mu = " << mu << ", kappa = " << kappa + << ", zeta = " << zeta << ", n = " << n; + } #define OP_CLASS Adhesion #define OP_VARS X(d) X(mu) X(kappa) X(zeta) X(n) From e545217916330ebcc82dfb339a9a7d920722f5b8 Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Tue, 1 Apr 2025 13:49:54 -0500 Subject: [PATCH 06/18] Modifying in clang format and documentation and adding regression tests and editor config --- src/Model/Solid/Finite/Adhesion.H | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 1357767fd1..f6dfe64381 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -189,8 +189,8 @@ public: Eigen::Matrix3d FinvT = F.inverse().transpose(); Eigen::Matrix3d dw = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3) * FinvT / J23) - + d * kappa * (J - 1.0) * J * FinvT - - zeta * n * std::pow(J, -n) * FinvT; + + d * kappa * (J - 1.0) * J * FinvT + - zeta * n * std::pow(J, -n) * FinvT; #if AMREX_SPACEDIM == 2 Set::Matrix r_dw; r_dw(0, 0) = dw(0, 0); @@ -242,10 +242,10 @@ public: t1 += (1.0 / 3.0) * (trace_FTF - 3) * FinvT(i, l) * FinvT(k, j); // Bulk-related term (kappa) t2 += (2.0 * J * J - J) * FinvT(i, j) * FinvT(k, l) - + (J - J * J) * FinvT(i, l) * FinvT(k, j); + + (J - J * J) * FinvT(i, l) * FinvT(k, j); // Adhesion-related term (zeta) t3 += (n * n) * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) - + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); + + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); ddw(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; } } From b09fa4456181bea366bfc03904e5daf7653671f6 Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 3 Apr 2025 11:34:06 -0500 Subject: [PATCH 07/18] Documentation and adding regression tests and clang format --- src/Model/Solid/Finite/Adhesion.H | 82 ++++++++------- tests/Solid/input | 159 +++++------------------------- 2 files changed, 72 insertions(+), 169 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index f6dfe64381..1c35a7bc30 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -1,4 +1,3 @@ -// // This models a finite-strain, solid with adhesion behavior governed by // a combination of elastic energy and an adhesion decay term. // @@ -30,7 +29,6 @@ // 1. Degradable (elastic shear + volumetric terms) // 2. Non-degradable (adhesion term decaying with volume expansion) // This model is rate-independent, reversible, and purely elastic. -// #ifndef MODEL_SOLID_FINITE_ADHESION_H_ #define MODEL_SOLID_FINITE_ADHESION_H_ @@ -39,7 +37,7 @@ #include "Model/Solid/Solid.H" #include "Set/Set.H" #include // For std::array to represent fourth-order tensors -#include // For std::pow and std::exp +#include // For std::pow and std::abs #include // For Eigen::Matrix3d namespace Set @@ -88,12 +86,13 @@ public: virtual ~Adhesion() {} // Material parameters (declared in order) - Set::Scalar d = NAN; ///< Adhesion strength + Set::Scalar d = NAN; ///< Adhesion strength (degradable part coefficient) Set::Scalar mu = NAN; ///< Shear modulus Set::Scalar kappa = NAN; ///< Bulk modulus - Set::Scalar zeta = NAN; ///< Adhesion decay parameter + Set::Scalar zeta = NAN; ///< Adhesion decay parameter (non-degradable part) Set::Scalar n = NAN; ///< Decay exponent + // Return a zero-parameter adhesion model. static Adhesion Zero() { @@ -106,6 +105,7 @@ public: return adhesion; } + // Return a random-parameter adhesion model. static Adhesion Random() { @@ -129,7 +129,7 @@ public: Set::Scalar E, nu; pp_query("E", E); // Elastic modulus pp_query("nu", nu); // Poisson's ratio - // Map E and nu to shear and bulk moduli (using a convention similar to NeoHookean citeturn0file2) + // Map E and nu to shear and bulk moduli (using a convention similar to NeoHookean :contentReference[oaicite:0]{index=0}) value.mu = E / (2.0 * (1.0 + nu)); value.kappa = E / (3.0 * (1.0 - 2.0 * nu)); // Consume F0 if provided, even though Adhesion does not use it. @@ -150,6 +150,8 @@ public: pp_query("n", value.n); // Decay exponent } + // Strain energy function: + // W = d * [ 0.5*mu*(trace(F^T F)-3)/J^(2/3) + 0.5*kappa*(J-1)^2 ] + zeta/J^n Set::Scalar W(const Set::Matrix &a_F) const override { @@ -163,15 +165,16 @@ public: Eigen::Matrix3d F = a_F; #endif Set::Scalar J = F.determinant(); - Set::Scalar J23 = std::pow(fabs(J), 2.0 / 3.0); - Set::Scalar w = 0.0; - // Shear and bulk contributions (scaled by adhesion factor d) - w += d * (0.5 * mu * (((F.transpose() * F).trace()) - 3) / J23 + 0.5 * kappa * (J - 1.0) * (J - 1.0)); - // Adhesion decay term - w += zeta / std::pow(J, n); - return w; + Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); + Set::Scalar energy = 0.0; + // Degradable part: shear and volumetric terms (scaled by d) + energy += d * (0.5 * mu * ((F.transpose() * F).trace() - 3.0) / J23 + 0.5 * kappa * (J - 1.0) * (J - 1.0)); + // Non-degradable (adhesion) part: decays with volume expansion + energy += zeta / std::pow(J, n); + return energy; } + // First Piola-Kirchhoff stress (derivative of W) Set::Matrix DW(const Set::Matrix &a_F) const override { @@ -185,24 +188,27 @@ public: Eigen::Matrix3d F = a_F; #endif Set::Scalar J = F.determinant(); - Set::Scalar J23 = std::pow(fabs(J), 2.0 / 3.0); + Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); Eigen::Matrix3d FinvT = F.inverse().transpose(); - Eigen::Matrix3d dw = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3) * FinvT / J23) - + d * kappa * (J - 1.0) * J * FinvT - - zeta * n * std::pow(J, -n) * FinvT; + // Derivative from degradable elastic part and bulk contribution + Eigen::Matrix3d stress = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3.0) * FinvT / J23) + + d * kappa * (J - 1.0) * J * FinvT; + // Adhesion term contribution + stress -= zeta * n * std::pow(J, -n) * FinvT; #if AMREX_SPACEDIM == 2 - Set::Matrix r_dw; - r_dw(0, 0) = dw(0, 0); - r_dw(0, 1) = dw(0, 1); - r_dw(1, 0) = dw(1, 0); - r_dw(1, 1) = dw(1, 1); - return r_dw; + Set::Matrix r_stress; + r_stress(0, 0) = stress(0, 0); + r_stress(0, 1) = stress(0, 1); + r_stress(1, 0) = stress(1, 0); + r_stress(1, 1) = stress(1, 1); + return r_stress; #elif AMREX_SPACEDIM == 3 - return dw; + return stress; #endif } + // Tangent stiffness tensor (second derivative of W) Set::Matrix4 DDW(const Set::Matrix &a_F) const override { @@ -216,13 +222,14 @@ public: Eigen::Matrix3d F = a_F; #endif Set::Scalar J = F.determinant(); - Set::Scalar J23 = std::pow(fabs(J), 2.0 / 3.0); + Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); Eigen::Matrix3d FinvT = F.inverse().transpose(); Set::Scalar trace_FTF = (F.transpose() * F).trace(); + // Precompute derivative helper for FinvT auto Deriv = ::Set::Derivative(FinvT); - Set::Matrix4<3, Set::Sym::Major> ddw; + Set::Matrix4<3, Set::Sym::Major> stiffness; for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) @@ -231,39 +238,40 @@ public: { for (int l = 0; l < 3; ++l) { - ddw(i, j, k, l) = 0.0; + stiffness(i, j, k, l) = 0.0; Set::Scalar t1 = 0.0, t2 = 0.0, t3 = 0.0; // Shear-related term (mu) if (i == k && j == l) t1 += 1.0; t1 -= (2.0 / 3.0) * F(i, j) * FinvT(k, l); t1 -= (2.0 / 3.0) * FinvT(i, j) * F(k, l); - t1 += (2.0 / 9.0) * (trace_FTF - 3) * FinvT(i, j) * FinvT(k, l); - t1 += (1.0 / 3.0) * (trace_FTF - 3) * FinvT(i, l) * FinvT(k, j); + t1 += (2.0 / 9.0) * (trace_FTF - 3.0) * FinvT(i, j) * FinvT(k, l); + t1 += (1.0 / 3.0) * (trace_FTF - 3.0) * FinvT(i, l) * FinvT(k, j); // Bulk-related term (kappa) t2 += (2.0 * J * J - J) * FinvT(i, j) * FinvT(k, l) - + (J - J * J) * FinvT(i, l) * FinvT(k, j); + + (J - J * J) * FinvT(i, l) * FinvT(k, j); // Adhesion-related term (zeta) t3 += (n * n) * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) - + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); - ddw(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; + + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); + stiffness(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; } } } } #if AMREX_SPACEDIM == 2 - Set::Matrix4 r_ddw; + Set::Matrix4 r_stiffness; for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) for (int k = 0; k < 2; k++) for (int l = 0; l < 2; l++) - r_ddw(i, j, k, l) = ddw(i, j, k, l); - return r_ddw; + r_stiffness(i, j, k, l) = stiffness(i, j, k, l); + return r_stiffness; #elif AMREX_SPACEDIM == 3 - return ddw; + return stiffness; #endif } + // Print material parameters virtual void Print(std::ostream &out) const override { @@ -282,4 +290,4 @@ public: } // namespace Solid } // namespace Model -#endif +#endif // MODEL_SOLID_FINITE_ADHESION_H_ diff --git a/tests/Solid/input b/tests/Solid/input index fac0d54c40..7e7a9c4519 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -1,142 +1,37 @@ -#@ -#@ [linear-isotropic] -#@ exe=mechanics -#@ dim = 3 -#@ check = false -#@ -#@ [linear-cubic] -#@ exe=mechanics -#@ dim = 3 -#@ check-file = reference/linear-cubic.dat -#@ args = alamo.program.mechanics.model=linear.cubic -#@ args = model1.C11 = 168.3 -#@ args = model1.C12 = 1.221 -#@ args = model1.C44 = 0.757 -#@ ignore = model1.E model1.nu -#@ coverage = true -#@ -#@ [affine-cubic] -#@ exe=mechanics -#@ dim = 3 -#@ check-file = reference/affine-cubic.dat -#@ args = alamo.program.mechanics.model=affine.cubic -#@ args = model1.C11 = 168.3 -#@ args = model1.C12 = 1.221 -#@ args = model1.C44 = 0.757 -#@ args = model1.F0=0.001 0 0 0 0.001 0 0 0 0.001 -#@ ignore = model1.E model1.nu -#@ coverage = true -#@ -#@ -#@ [affine-hexagonal] -#@ exe=mechanics -#@ dim = 3 -#@ check-file = reference/affine-hexagonal.dat -#@ args = alamo.program.mechanics.model=affine.hexagonal -#@ args = model1.C11 = 0.597 -#@ args = model1.C12 = 0.262 -#@ args = model1.C13 = 0.217 -#@ args = model1.C33 = 0.617 -#@ args = model1.C44 = 0.164 -#@ args = model1.F0=0.001 0 0 0 0.001 0 0 0 0.001 -#@ ignore = model1.E model1.nu -#@ coverage = true -#@ -#@ [j2] -#@ exe=mechanics -#@ dim = 3 -#@ check-file = reference/j2.dat -#@ args = timestep = 0.001 -#@ args = alamo.program.mechanics.model=affine.j2 -#@ args = model1.E=210 -#@ args = model1.nu=0.3 -#@ args = model1.sigma0=0.2 -#@ args = bc.tensiontest.disp = (0,0.25,0.75,1.0:0,0.002,-0.002,0) -#@ benchmark-beaker = 10.42 -#@ benchmark-statler = 7.02 -#@ benchmark-github = 11.09 -#@ -#@ [affine-isotropic] -#@ exe=mechanics -#@ dim = 3 -#@ check = false -#@ args = alamo.program.mechanics.model=affine.isotropic -#@ args = model1.F0=0.001 0 0 0 0.001 0 0 0 0.001 -#@ args = model1.E=210 -#@ args = model1.nu=0.3 -#@ coverage=true -#@ -#@ [neo-hookean] -#@ exe=mechanics -#@ dim=3 -#@ check-file = reference/neo-hookean.dat -#@ args = timestep=0.01 -#@ args = alamo.program.mechanics.model=finite.neohookean -#@ args = model1.mu=3.0 -#@ args = model1.kappa=6.5 -#@ args = bc.tensiontest.disp=(0,1:0,1) -#@ args = solver.nriters=10 -#@ ignore = model1.E model1.nu -#@ benchmark-beaker = 10.73 -#@ benchmark-statler = 8.17 -#@ benchmark-github = 11.35 -#@ -#@ [neo-hookean-2d] -#@ exe=mechanics -#@ dim=2 -#@ check-file = reference/neo-hookean-2d.dat -#@ args = timestep=0.01 -#@ args = alamo.program.mechanics.model=finite.neohookean -#@ args = model1.mu=3.0 -#@ args = model1.kappa=6.5 -#@ args = bc.tensiontest.disp=(0,1:0,1) -#@ args = solver.nriters=10 -#@ ignore = model1.E model1.nu -#@ coverage = true -#@ -#@ [pseudolinear-cubic] -#@ exe=mechanics -#@ dim = 3 -#@ check-file = reference/pseudolinear-cubic.dat -#@ args = alamo.program.mechanics.model=finite.pseudolinear.cubic -#@ args = model1.C11=168.3 -#@ args = model1.C12=1.221 -#@ args = model1.C44=0.757 -#@ ignore = model1.E model1.nu -#@ coverage=true -#@ +#@ #@[linear - isotropic] #@ exe = mechanics #@ dim = 3 #@ check = false #@ #@[linear - cubic] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / linear - cubic.dat #@ args = alamo.program.mechanics.model = linear.cubic #@ args = model1.C11 = 168.3 #@ args = model1.C12 = 1.221 #@ args = model1.C44 = 0.757 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[affine - cubic] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / affine - cubic.dat #@ args = alamo.program.mechanics.model = affine.cubic #@ args = model1.C11 = 168.3 #@ args = model1.C12 = 1.221 #@ args = model1.C44 = 0.757 #@ args = model1.F0 = 0.001 0 0 0 0.001 0 0 0 0.001 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@ #@[affine - hexagonal] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / affine - hexagonal.dat #@ args = alamo.program.mechanics.model = affine.hexagonal #@ args = model1.C11 = 0.597 #@ args = model1.C12 = 0.262 #@ args = model1.C13 = 0.217 #@ args = model1.C33 = 0.617 #@ args = model1.C44 = 0.164 #@ args = model1.F0 = 0.001 0 0 0 0.001 0 0 0 0.001 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[j2] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / j2.dat #@ args = timestep = 0.001 #@ args = alamo.program.mechanics.model = affine.j2 #@ args = model1.E = 210 #@ args = model1.nu = 0.3 #@ args = model1.sigma0 = 0.2 #@ args = bc.tensiontest.disp = (0, 0.25, 0.75, 1.0 : 0, 0.002, -0.002, 0) #@ benchmark - beaker = 10.42 #@ benchmark - statler = 7.02 #@ benchmark - github = 11.09 #@ #@[affine - isotropic] #@ exe = mechanics #@ dim = 3 #@ check = false #@ args = alamo.program.mechanics.model = affine.isotropic #@ args = model1.F0 = 0.001 0 0 0 0.001 0 0 0 0.001 #@ args = model1.E = 210 #@ args = model1.nu = 0.3 #@ coverage = true #@ #@[neo - hookean] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / neo - hookean.dat #@ args = timestep = 0.01 #@ args = alamo.program.mechanics.model = finite.neohookean #@ args = model1.mu = 3.0 #@ args = model1.kappa = 6.5 #@ args = bc.tensiontest.disp = (0, 1 : 0, 1) #@ args = solver.nriters = 10 #@ ignore = model1.E model1.nu #@ benchmark - beaker = 10.73 #@ benchmark - statler = 8.17 #@ benchmark - github = 11.35 #@ #@[neo - hookean - 2d] #@ exe = mechanics #@ dim = 2 #@ check - file = reference / neo - hookean - 2d.dat #@ args = timestep = 0.01 #@ args = alamo.program.mechanics.model = finite.neohookean #@ args = model1.mu = 3.0 #@ args = model1.kappa = 6.5 #@ args = bc.tensiontest.disp = (0, 1 : 0, 1) #@ args = solver.nriters = 10 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[pseudolinear - cubic] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / pseudolinear - cubic.dat #@ args = alamo.program.mechanics.model = finite.pseudolinear.cubic #@ args = model1.C11 = 168.3 #@ args = model1.C12 = 1.221 #@ args = model1.C44 = 0.757 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[adhesion] #@ exe = mechanics #@ check = false #@ dim = 3 #@ args = timestep = 0.01 #@ args = alamo.program.mechanics.model = finite.adhesion #@ args = model1.mu = 3.0 #@ args = model1.kappa = 6.5 #@ args = model1.zeta = 0.01 #@ args = model1.d = 1 #@ args = bc.tensiontest.disp = (0, 1 : 0, 1) #@ args = solver.nriters = 10 #@ ignore = model1.E model1.nu #@ benchmark - beaker = 10.73 #@ benchmark - statler = 8.17 #@ benchmark - github = 11.35 -alamo.program = mechanics + alamo.program + = mechanics + plot_file + = tests / Solid / output -plot_file = tests/Solid/output + type + = static -type=static + timestep + = 0.01 stop_time = 0.01 -timestep = 0.01 -stop_time = 1.0 +#amr parameters + amr.plot_dt + = 0.1 amr.max_level = 0 amr.n_cell = 4 4 4 amr.blocking_factor = 2 -# amr parameters -amr.plot_dt = 0.1 -amr.max_level = 0 -amr.n_cell = 4 4 4 -amr.blocking_factor = 2 + amr.thermo.int + = 1 amr.thermo.plot_int = 1 -amr.thermo.int = 1 -amr.thermo.plot_int = 1 +#geometry + geometry.prob_lo + = 0 0 0 geometry.prob_hi = 1 1 1 -# geometry -geometry.prob_lo = 0 0 0 -geometry.prob_hi = 1 1 1 + nmodels + = 1 model1.E = 210 model1.nu = 0.3 -nmodels = 1 -model1.E=210 -model1.nu=0.3 + solver.verbose + = 3 solver.nriters = 1 solver.max_iter = 30 #30 -solver.verbose = 3 -solver.nriters = 1 -solver.max_iter = 30 #30 - -bc.type = tensiontest -bc.tensiontest.type = uniaxial_stress -bc.tensiontest.disp=(0,1:0,0.01) + bc.type + = tensiontest + bc.tensiontest.type + = uniaxial_stress + bc.tensiontest.disp + = (0, 1 : 0, 0.01) From be5570f9f018fc50d6d0dd356797d54f491b3ac4 Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 3 Apr 2025 11:49:39 -0500 Subject: [PATCH 08/18] Documentation of the adhesion and adding the regression tests and clang format --- tests/Solid/input | 178 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 151 insertions(+), 27 deletions(-) diff --git a/tests/Solid/input b/tests/Solid/input index 7e7a9c4519..77fd7428d2 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -1,37 +1,161 @@ -#@ #@[linear - isotropic] #@ exe = mechanics #@ dim = 3 #@ check = false #@ #@[linear - cubic] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / linear - cubic.dat #@ args = alamo.program.mechanics.model = linear.cubic #@ args = model1.C11 = 168.3 #@ args = model1.C12 = 1.221 #@ args = model1.C44 = 0.757 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[affine - cubic] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / affine - cubic.dat #@ args = alamo.program.mechanics.model = affine.cubic #@ args = model1.C11 = 168.3 #@ args = model1.C12 = 1.221 #@ args = model1.C44 = 0.757 #@ args = model1.F0 = 0.001 0 0 0 0.001 0 0 0 0.001 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@ #@[affine - hexagonal] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / affine - hexagonal.dat #@ args = alamo.program.mechanics.model = affine.hexagonal #@ args = model1.C11 = 0.597 #@ args = model1.C12 = 0.262 #@ args = model1.C13 = 0.217 #@ args = model1.C33 = 0.617 #@ args = model1.C44 = 0.164 #@ args = model1.F0 = 0.001 0 0 0 0.001 0 0 0 0.001 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[j2] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / j2.dat #@ args = timestep = 0.001 #@ args = alamo.program.mechanics.model = affine.j2 #@ args = model1.E = 210 #@ args = model1.nu = 0.3 #@ args = model1.sigma0 = 0.2 #@ args = bc.tensiontest.disp = (0, 0.25, 0.75, 1.0 : 0, 0.002, -0.002, 0) #@ benchmark - beaker = 10.42 #@ benchmark - statler = 7.02 #@ benchmark - github = 11.09 #@ #@[affine - isotropic] #@ exe = mechanics #@ dim = 3 #@ check = false #@ args = alamo.program.mechanics.model = affine.isotropic #@ args = model1.F0 = 0.001 0 0 0 0.001 0 0 0 0.001 #@ args = model1.E = 210 #@ args = model1.nu = 0.3 #@ coverage = true #@ #@[neo - hookean] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / neo - hookean.dat #@ args = timestep = 0.01 #@ args = alamo.program.mechanics.model = finite.neohookean #@ args = model1.mu = 3.0 #@ args = model1.kappa = 6.5 #@ args = bc.tensiontest.disp = (0, 1 : 0, 1) #@ args = solver.nriters = 10 #@ ignore = model1.E model1.nu #@ benchmark - beaker = 10.73 #@ benchmark - statler = 8.17 #@ benchmark - github = 11.35 #@ #@[neo - hookean - 2d] #@ exe = mechanics #@ dim = 2 #@ check - file = reference / neo - hookean - 2d.dat #@ args = timestep = 0.01 #@ args = alamo.program.mechanics.model = finite.neohookean #@ args = model1.mu = 3.0 #@ args = model1.kappa = 6.5 #@ args = bc.tensiontest.disp = (0, 1 : 0, 1) #@ args = solver.nriters = 10 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[pseudolinear - cubic] #@ exe = mechanics #@ dim = 3 #@ check - file = reference / pseudolinear - cubic.dat #@ args = alamo.program.mechanics.model = finite.pseudolinear.cubic #@ args = model1.C11 = 168.3 #@ args = model1.C12 = 1.221 #@ args = model1.C44 = 0.757 #@ ignore = model1.E model1.nu #@ coverage = true #@ #@[adhesion] #@ exe = mechanics #@ check = false #@ dim = 3 #@ args = timestep = 0.01 #@ args = alamo.program.mechanics.model = finite.adhesion #@ args = model1.mu = 3.0 #@ args = model1.kappa = 6.5 #@ args = model1.zeta = 0.01 #@ args = model1.d = 1 #@ args = bc.tensiontest.disp = (0, 1 : 0, 1) #@ args = solver.nriters = 10 #@ ignore = model1.E model1.nu #@ benchmark - beaker = 10.73 #@ benchmark - statler = 8.17 #@ benchmark - github = 11.35 +#@ +#@ [linear-isotropic] +#@ exe=mechanics +#@ dim = 3 +#@ check = false +#@ +#@ [linear-cubic] +#@ exe=mechanics +#@ dim = 3 +#@ check-file = reference/linear-cubic.dat +#@ args = alamo.program.mechanics.model=linear.cubic +#@ args = model1.C11 = 168.3 +#@ args = model1.C12 = 1.221 +#@ args = model1.C44 = 0.757 +#@ ignore = model1.E model1.nu +#@ coverage = true +#@ +#@ [affine-cubic] +#@ exe=mechanics +#@ dim = 3 +#@ check-file = reference/affine-cubic.dat +#@ args = alamo.program.mechanics.model=affine.cubic +#@ args = model1.C11 = 168.3 +#@ args = model1.C12 = 1.221 +#@ args = model1.C44 = 0.757 +#@ args = model1.F0=0.001 0 0 0 0.001 0 0 0 0.001 +#@ ignore = model1.E model1.nu +#@ coverage = true +#@ +#@ +#@ [affine-hexagonal] +#@ exe=mechanics +#@ dim = 3 +#@ check-file = reference/affine-hexagonal.dat +#@ args = alamo.program.mechanics.model=affine.hexagonal +#@ args = model1.C11 = 0.597 +#@ args = model1.C12 = 0.262 +#@ args = model1.C13 = 0.217 +#@ args = model1.C33 = 0.617 +#@ args = model1.C44 = 0.164 +#@ args = model1.F0=0.001 0 0 0 0.001 0 0 0 0.001 +#@ ignore = model1.E model1.nu +#@ coverage = true +#@ +#@ [j2] +#@ exe=mechanics +#@ dim = 3 +#@ check-file = reference/j2.dat +#@ args = timestep = 0.001 +#@ args = alamo.program.mechanics.model=affine.j2 +#@ args = model1.E=210 +#@ args = model1.nu=0.3 +#@ args = model1.sigma0=0.2 +#@ args = bc.tensiontest.disp = (0,0.25,0.75,1.0:0,0.002,-0.002,0) +#@ benchmark-beaker = 10.42 +#@ benchmark-statler = 7.02 +#@ benchmark-github = 11.09 +#@ +#@ [affine-isotropic] +#@ exe=mechanics +#@ dim = 3 +#@ check = false +#@ args = alamo.program.mechanics.model=affine.isotropic +#@ args = model1.F0=0.001 0 0 0 0.001 0 0 0 0.001 +#@ args = model1.E=210 +#@ args = model1.nu=0.3 +#@ coverage=true +#@ +#@ [neo-hookean] +#@ exe=mechanics +#@ dim=3 +#@ check-file = reference/neo-hookean.dat +#@ args = timestep=0.01 +#@ args = alamo.program.mechanics.model=finite.neohookean +#@ args = model1.mu=3.0 +#@ args = model1.kappa=6.5 +#@ args = bc.tensiontest.disp=(0,1:0,1) +#@ args = solver.nriters=10 +#@ ignore = model1.E model1.nu +#@ benchmark-beaker = 10.73 +#@ benchmark-statler = 8.17 +#@ benchmark-github = 11.35 +#@ +#@ [neo-hookean-2d] +#@ exe=mechanics +#@ dim=2 +#@ check-file = reference/neo-hookean-2d.dat +#@ args = timestep=0.01 +#@ args = alamo.program.mechanics.model=finite.neohookean +#@ args = model1.mu=3.0 +#@ args = model1.kappa=6.5 +#@ args = bc.tensiontest.disp=(0,1:0,1) +#@ args = solver.nriters=10 +#@ ignore = model1.E model1.nu +#@ coverage = true +#@ +#@ [pseudolinear-cubic] +#@ exe=mechanics +#@ dim = 3 +#@ check-file = reference/pseudolinear-cubic.dat +#@ args = alamo.program.mechanics.model=finite.pseudolinear.cubic +#@ args = model1.C11=168.3 +#@ args = model1.C12=1.221 +#@ args = model1.C44=0.757 +#@ ignore = model1.E model1.nu +#@ coverage=true +#@ +#@ [adhesion] +#@ exe=mechanics +#@ check = false +#@ dim=3 +#@ args = timestep=0.01 +#@ args = alamo.program.mechanics.model=finite.adhesion +#@ args = model1.mu=3.0 +#@ args = model1.kappa=6.5 +#@ args = model1.zeta=0.01 +#@ args = model1.d=1 +#@ args = bc.tensiontest.disp=(0,1:0,1) +#@ args = solver.nriters=10 +#@ ignore = model1.E model1.nu +#@ benchmark-beaker = 10.73 +#@ benchmark-statler = 8.17 +#@ benchmark-github = 11.35 - alamo.program - = mechanics - plot_file - = tests / Solid / output - type - = static - timestep - = 0.01 stop_time = 0.01 +alamo.program = mechanics -#amr parameters - amr.plot_dt - = 0.1 amr.max_level = 0 amr.n_cell = 4 4 4 amr.blocking_factor = 2 - amr.thermo.int - = 1 amr.thermo.plot_int = 1 +plot_file = tests/Solid/output -#geometry - geometry.prob_lo - = 0 0 0 geometry.prob_hi = 1 1 1 +type=static - nmodels - = 1 model1.E = 210 model1.nu = 0.3 +timestep = 0.01 +stop_time = 0.01 - solver.verbose - = 3 solver.nriters = 1 solver.max_iter = 30 #30 +# amr parameters +amr.plot_dt = 0.1 +amr.max_level = 0 +amr.n_cell = 4 4 4 +amr.blocking_factor = 2 - bc.type - = tensiontest - bc.tensiontest.type - = uniaxial_stress - bc.tensiontest.disp - = (0, 1 : 0, 0.01) +amr.thermo.int = 1 +amr.thermo.plot_int = 1 + +# geometry +geometry.prob_lo = 0 0 0 +geometry.prob_hi = 1 1 1 + +nmodels = 1 +model1.E=210 +model1.nu=0.3 + +solver.verbose = 3 +solver.nriters = 1 +solver.max_iter = 30 #30 + +bc.type = tensiontest +bc.tensiontest.type = uniaxial_stress +bc.tensiontest.disp=(0,1:0,0.01) From 58ace0e2266e87d299561ba9e70ba8ec2cb2345c Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 3 Apr 2025 12:06:14 -0500 Subject: [PATCH 09/18] Documentation of Adhesion model and adding regression tests and clang format --- tests/Solid/input | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/Solid/input b/tests/Solid/input index 77fd7428d2..1733df872a 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -114,13 +114,16 @@ #@ args = model1.mu=3.0 #@ args = model1.kappa=6.5 #@ args = model1.zeta=0.01 -#@ args = model1.d=1 +#@ args = model1.d=0.01 +#@ args = model1.n=1 #@ args = bc.tensiontest.disp=(0,1:0,1) #@ args = solver.nriters=10 #@ ignore = model1.E model1.nu #@ benchmark-beaker = 10.73 #@ benchmark-statler = 8.17 #@ benchmark-github = 11.35 +#@ + From 9c944f4fd1198f280b0db6d5dd1ffb969ded9b6d Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Thu, 3 Apr 2025 12:59:29 -0500 Subject: [PATCH 10/18] Documentation of Adhesion model and adding regression tests in clang format --- src/Model/Solid/Finite/Adhesion.H | 2 +- tests/Solid/input | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 1c35a7bc30..2813facdc2 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -290,4 +290,4 @@ public: } // namespace Solid } // namespace Model -#endif // MODEL_SOLID_FINITE_ADHESION_H_ +#endif // MODEL_SOLID_FINITE_ADHESION_H_ \ No newline at end of file diff --git a/tests/Solid/input b/tests/Solid/input index 1733df872a..03b7f39066 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -115,7 +115,6 @@ #@ args = model1.kappa=6.5 #@ args = model1.zeta=0.01 #@ args = model1.d=0.01 -#@ args = model1.n=1 #@ args = bc.tensiontest.disp=(0,1:0,1) #@ args = solver.nriters=10 #@ ignore = model1.E model1.nu From bd171a6c8161239efd5be36dbfc958b23770bfaf Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Mon, 14 Apr 2025 13:29:54 -0500 Subject: [PATCH 11/18] Adhesion model with the added test --- src/Model/Solid/Finite/Adhesion.H | 57 +++++----------- tests/Solid/input | 12 ++-- tests/Solid/reference/adhesion.dat | 101 +++++++++++++++++++++++++++++ 3 files changed, 126 insertions(+), 44 deletions(-) create mode 100644 tests/Solid/reference/adhesion.dat diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 2813facdc2..3b82d42cd2 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -1,3 +1,13 @@ +#ifndef MODEL_SOLID_FINITE_ADHESION_H_ +#define MODEL_SOLID_FINITE_ADHESION_H_ + +#include "IO/ParmParse.H" +#include "Model/Solid/Solid.H" +#include "Set/Set.H" +#include // For std::array to represent fourth-order tensors +#include // For std::pow and std::abs +#include // For Eigen::Matrix3d + // This models a finite-strain, solid with adhesion behavior governed by // a combination of elastic energy and an adhesion decay term. // @@ -30,16 +40,6 @@ // 2. Non-degradable (adhesion term decaying with volume expansion) // This model is rate-independent, reversible, and purely elastic. -#ifndef MODEL_SOLID_FINITE_ADHESION_H_ -#define MODEL_SOLID_FINITE_ADHESION_H_ - -#include "IO/ParmParse.H" -#include "Model/Solid/Solid.H" -#include "Set/Set.H" -#include // For std::array to represent fourth-order tensors -#include // For std::pow and std::abs -#include // For Eigen::Matrix3d - namespace Set { @@ -118,36 +118,16 @@ public: return adhesion; } - // Updated Parse function: - // If "E" and "nu" (and optionally "F0") are provided, use these to compute mu and kappa. - // Otherwise, query for "mu" and "kappa" directly. + // Parse function: Reads material parameters from the input file. + // Here, we assume the input file provides: mu, kappa, d, zeta, and n. static void Parse(Adhesion &value, IO::ParmParse &pp) { - if (pp.contains("E") && pp.contains("nu")) - { - Set::Scalar E, nu; - pp_query("E", E); // Elastic modulus - pp_query("nu", nu); // Poisson's ratio - // Map E and nu to shear and bulk moduli (using a convention similar to NeoHookean :contentReference[oaicite:0]{index=0}) - value.mu = E / (2.0 * (1.0 + nu)); - value.kappa = E / (3.0 * (1.0 - 2.0 * nu)); - // Consume F0 if provided, even though Adhesion does not use it. - if (pp.contains("F0")) - { - Set::Matrix dummy; - pp_queryarr("F0", dummy); // F0 - } - } - else - { - pp_query("mu", value.mu); // Shear modulus - pp_query("kappa", value.kappa); // Bulk modulus - } - // Query the remaining Adhesion-specific parameters. - pp_query("d", value.d); // Adhesion strength - pp_query("zeta", value.zeta); // Adhesion decay parameter - pp_query("n", value.n); // Decay exponent + pp_query("mu", value.mu); + pp_query("kappa", value.kappa); + pp_query("d", value.d); + pp_query("zeta", value.zeta); + pp_query("n", value.n); } // Strain energy function: @@ -190,7 +170,6 @@ public: Set::Scalar J = F.determinant(); Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); Eigen::Matrix3d FinvT = F.inverse().transpose(); - // Derivative from degradable elastic part and bulk contribution Eigen::Matrix3d stress = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3.0) * FinvT / J23) + d * kappa * (J - 1.0) * J * FinvT; @@ -290,4 +269,4 @@ public: } // namespace Solid } // namespace Model -#endif // MODEL_SOLID_FINITE_ADHESION_H_ \ No newline at end of file +#endif // MODEL_SOLID_FINITE_ADHESION_H_ diff --git a/tests/Solid/input b/tests/Solid/input index 03b7f39066..cb86578f7f 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -107,21 +107,23 @@ #@ #@ [adhesion] #@ exe=mechanics -#@ check = false #@ dim=3 +#@ check = true +#@ check-file = reference/adhesion.dat #@ args = timestep=0.01 #@ args = alamo.program.mechanics.model=finite.adhesion #@ args = model1.mu=3.0 #@ args = model1.kappa=6.5 -#@ args = model1.zeta=0.01 -#@ args = model1.d=0.01 +#@ args = model1.d=1.0 +#@ args = model1.zeta=0.5 +#@ args = model1.n=2.0 #@ args = bc.tensiontest.disp=(0,1:0,1) #@ args = solver.nriters=10 #@ ignore = model1.E model1.nu #@ benchmark-beaker = 10.73 #@ benchmark-statler = 8.17 #@ benchmark-github = 11.35 -#@ +#@ @@ -135,7 +137,7 @@ plot_file = tests/Solid/output type=static timestep = 0.01 -stop_time = 0.01 +stop_time = 1 # amr parameters amr.plot_dt = 0.1 diff --git a/tests/Solid/reference/adhesion.dat b/tests/Solid/reference/adhesion.dat new file mode 100644 index 0000000000..6dc6910767 --- /dev/null +++ b/tests/Solid/reference/adhesion.dat @@ -0,0 +1,101 @@ +time disp_xhi_x disp_xhi_y disp_yhi_x disp_yhi_y trac_xhi_x trac_xhi_y trac_yhi_x trac_yhi_y +0 0 -0.123824 0 0 1.15602 -4.01326e-11 0 0 +0.01 0.01 -0.12456 0 0 1.2205 -2.86443e-10 0 0 +0.02 0.02 -0.125279 0 0 1.28357 -3.93683e-10 0 0 +0.03 0.03 -0.125982 0 0 1.34527 -3.4169e-10 0 0 +0.04 0.04 -0.126672 0 0 1.40568 -3.48891e-10 0 0 +0.05 0.05 -0.127348 0 0 1.46483 -3.40385e-10 0 0 +0.06 0.06 -0.128012 0 0 1.52278 -2.15224e-10 0 0 +0.07 0.07 -0.128666 0 0 1.57958 -2.15809e-10 0 0 +0.08 0.08 -0.129309 0 0 1.63528 -2.85579e-10 0 0 +0.09 0.09 -0.129944 0 0 1.68991 -4.1995e-10 0 0 +0.1 0.1 -0.130571 0 0 1.74353 -1.97542e-12 0 0 +0.11 0.11 -0.13119 0 0 1.79616 -7.2879e-12 0 0 +0.12 0.12 -0.131804 0 0 1.84785 -4.12464e-11 0 0 +0.13 0.13 -0.132411 0 0 1.89863 -1.15129e-13 0 0 +0.14 0.14 -0.133013 0 0 1.94854 -2.45439e-10 0 0 +0.15 0.15 -0.133611 0 0 1.99761 -2.35392e-10 0 0 +0.16 0.16 -0.134205 0 0 2.04587 -4.76978e-11 0 0 +0.17 0.17 -0.134795 0 0 2.09335 1.49388e-10 0 0 +0.18 0.18 -0.135382 0 0 2.14007 -4.91033e-12 0 0 +0.19 0.19 -0.135967 0 0 2.18607 -6.66379e-10 0 0 +0.2 0.2 -0.136549 0 0 2.23137 -8.69189e-11 0 0 +0.21 0.21 -0.137129 0 0 2.27599 4.86789e-10 0 0 +0.22 0.22 -0.137708 0 0 2.31996 -5.14368e-10 0 0 +0.23 0.23 -0.138285 0 0 2.3633 -1.30087e-12 0 0 +0.24 0.24 -0.138861 0 0 2.40602 -1.76863e-10 0 0 +0.25 0.25 -0.139436 0 0 2.44816 -3.11804e-10 0 0 +0.26 0.26 -0.14001 0 0 2.48972 -3.65531e-10 0 0 +0.27 0.27 -0.140583 0 0 2.53073 -2.95396e-10 0 0 +0.28 0.28 -0.141156 0 0 2.5712 3.44428e-11 0 0 +0.29 0.29 -0.141729 0 0 2.61115 -3.22309e-11 0 0 +0.3 0.3 -0.142301 0 0 2.6506 -1.88363e-10 0 0 +0.31 0.31 -0.142873 0 0 2.68956 -2.5663e-10 0 0 +0.32 0.32 -0.143445 0 0 2.72804 -2.13702e-12 0 0 +0.33 0.33 -0.144017 0 0 2.76606 2.24091e-10 0 0 +0.34 0.34 -0.144588 0 0 2.80363 -2.44375e-10 0 0 +0.35 0.35 -0.14516 0 0 2.84077 1.04304e-11 0 0 +0.36 0.36 -0.145731 0 0 2.87748 -2.56168e-10 0 0 +0.37 0.37 -0.146303 0 0 2.91379 -2.05313e-10 0 0 +0.38 0.38 -0.146874 0 0 2.94969 -2.77199e-10 0 0 +0.39 0.39 -0.147446 0 0 2.9852 7.41536e-12 0 0 +0.4 0.4 -0.148017 0 0 3.02033 -8.64127e-12 0 0 +0.41 0.41 -0.148588 0 0 3.0551 -3.00506e-10 0 0 +0.42 0.42 -0.149159 0 0 3.0895 4.60683e-11 0 0 +0.43 0.43 -0.14973 0 0 3.12355 -3.36724e-10 0 0 +0.44 0.44 -0.150301 0 0 3.15726 -5.79441e-11 0 0 +0.45 0.45 -0.150872 0 0 3.19063 -2.86736e-10 0 0 +0.46 0.46 -0.151442 0 0 3.22368 2.43842e-10 0 0 +0.47 0.47 -0.152012 0 0 3.25642 2.86899e-12 0 0 +0.48 0.48 -0.152582 0 0 3.28884 2.16865e-10 0 0 +0.49 0.49 -0.153151 0 0 3.32096 3.02515e-12 0 0 +0.5 0.5 -0.15372 0 0 3.35278 1.91218e-13 0 0 +0.51 0.51 -0.154288 0 0 3.38431 2.53284e-11 0 0 +0.52 0.52 -0.154856 0 0 3.41556 -2.02127e-12 0 0 +0.53 0.53 -0.155423 0 0 3.44654 2.64056e-11 0 0 +0.54 0.54 -0.15599 0 0 3.47724 -2.6267e-10 0 0 +0.55 0.55 -0.156556 0 0 3.50768 1.40424e-10 0 0 +0.56 0.56 -0.157122 0 0 3.53787 -1.07387e-10 0 0 +0.57 0.57 -0.157686 0 0 3.5678 -1.00735e-10 0 0 +0.58 0.58 -0.15825 0 0 3.59748 -1.29752e-10 0 0 +0.59 0.59 -0.158813 0 0 3.62692 -8.31481e-12 0 0 +0.6 0.6 -0.159375 0 0 3.65612 -9.11869e-12 0 0 +0.61 0.61 -0.159936 0 0 3.68509 -9.81089e-12 0 0 +0.62 0.62 -0.160497 0 0 3.71384 -1.25632e-11 0 0 +0.63 0.63 -0.161056 0 0 3.74236 -1.38536e-11 0 0 +0.64 0.64 -0.161614 0 0 3.77066 -1.15191e-11 0 0 +0.65 0.65 -0.162171 0 0 3.79875 -1.08418e-11 0 0 +0.66 0.66 -0.162727 0 0 3.82663 -1.18489e-11 0 0 +0.67 0.67 -0.163282 0 0 3.8543 -8.40617e-12 0 0 +0.68 0.68 -0.163835 0 0 3.88177 -1.031e-11 0 0 +0.69 0.69 -0.164388 0 0 3.90904 8.89774e-12 0 0 +0.7 0.7 -0.164939 0 0 3.93612 -1.94457e-11 0 0 +0.71 0.71 -0.165489 0 0 3.96301 -1.72192e-11 0 0 +0.72 0.72 -0.166037 0 0 3.98971 3.39921e-11 0 0 +0.73 0.73 -0.166584 0 0 4.01623 1.45924e-11 0 0 +0.74 0.74 -0.16713 0 0 4.04257 -2.67873e-12 0 0 +0.75 0.75 -0.167675 0 0 4.06873 -3.08766e-10 0 0 +0.76 0.76 -0.168217 0 0 4.09471 2.23796e-10 0 0 +0.77 0.77 -0.168759 0 0 4.12053 3.09863e-11 0 0 +0.78 0.78 -0.169299 0 0 4.14618 1.1446e-11 0 0 +0.79 0.79 -0.169837 0 0 4.17166 -2.76966e-10 0 0 +0.8 0.8 -0.170374 0 0 4.19698 3.90644e-11 0 0 +0.81 0.81 -0.170909 0 0 4.22215 -2.30649e-13 0 0 +0.82 0.82 -0.171443 0 0 4.24715 1.16604e-10 0 0 +0.83 0.83 -0.171975 0 0 4.27201 9.44755e-11 0 0 +0.84 0.84 -0.172506 0 0 4.29671 -4.09357e-13 0 0 +0.85 0.85 -0.173035 0 0 4.32126 -3.04696e-12 0 0 +0.86 0.86 -0.173562 0 0 4.34567 1.19393e-10 0 0 +0.87 0.87 -0.174088 0 0 4.36993 1.39926e-10 0 0 +0.88 0.88 -0.174611 0 0 4.39406 1.81844e-10 0 0 +0.89 0.89 -0.175134 0 0 4.41804 -6.50169e-11 0 0 +0.9 0.9 -0.175654 0 0 4.44188 3.89623e-11 0 0 +0.91 0.91 -0.176173 0 0 4.4656 -1.05935e-11 0 0 +0.92 0.92 -0.17669 0 0 4.48917 -1.63238e-11 0 0 +0.93 0.93 -0.177205 0 0 4.51262 -2.55951e-12 0 0 +0.94 0.94 -0.177719 0 0 4.53594 7.32922e-12 0 0 +0.95 0.95 -0.17823 0 0 4.55914 -1.55842e-11 0 0 +0.96 0.96 -0.17874 0 0 4.58221 -8.00096e-12 0 0 +0.97 0.97 -0.179248 0 0 4.60515 -1.05624e-11 0 0 +0.98 0.98 -0.179755 0 0 4.62798 1.43395e-10 0 0 +0.99 0.99 -0.180259 0 0 4.65069 -9.58419e-12 0 0 From 841ef64d3b1ee3af5aa991d29648471908bdeef5 Mon Sep 17 00:00:00 2001 From: Nima Sina Date: Mon, 14 Apr 2025 13:40:59 -0500 Subject: [PATCH 12/18] Adhesion model with added test --- src/Model/Solid/Finite/Adhesion.H | 36 +++++++++++++++---------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 3b82d42cd2..82a82b4c44 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -97,11 +97,11 @@ public: Zero() { Adhesion adhesion; - adhesion.d = 0; - adhesion.mu = 0; - adhesion.kappa = 0; - adhesion.zeta = 0; - adhesion.n = 0; + adhesion.d = 0; ///< Adhesion strength (degradable part coefficient) + adhesion.mu = 0; ///< Shear modulus + adhesion.kappa = 0; ///< Bulk modulus + adhesion.zeta = 0; ///< Adhesion decay parameter (non-degradable part) + adhesion.n = 0; ///< Decay exponent return adhesion; } @@ -110,11 +110,11 @@ public: Random() { Adhesion adhesion; - adhesion.d = Util::Random(); - adhesion.mu = Util::Random(); - adhesion.kappa = Util::Random(); - adhesion.zeta = Util::Random(); - adhesion.n = Util::Random(); + adhesion.d = Util::Random(); ///< Adhesion strength (degradable part coefficient) + adhesion.mu = Util::Random(); ///< Shear modulus + adhesion.kappa = Util::Random(); ///< Bulk modulus + adhesion.zeta = Util::Random(); ///< Adhesion decay parameter (non-degradable part) + adhesion.n = Util::Random(); ///< Decay exponent return adhesion; } @@ -123,11 +123,11 @@ public: static void Parse(Adhesion &value, IO::ParmParse &pp) { - pp_query("mu", value.mu); - pp_query("kappa", value.kappa); - pp_query("d", value.d); - pp_query("zeta", value.zeta); - pp_query("n", value.n); + pp_query("mu", value.mu); ///< Shear modulus + pp_query("kappa", value.kappa); ///< Bulk modulus + pp_query("d", value.d); ///< Adhesion strength (degradable part coefficient) + pp_query("zeta", value.zeta); ///< Adhesion decay parameter (non-degradable part) + pp_query("n", value.n); ///< Decay exponent } // Strain energy function: @@ -172,7 +172,7 @@ public: Eigen::Matrix3d FinvT = F.inverse().transpose(); // Derivative from degradable elastic part and bulk contribution Eigen::Matrix3d stress = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3.0) * FinvT / J23) - + d * kappa * (J - 1.0) * J * FinvT; + + d * kappa * (J - 1.0) * J * FinvT; // Adhesion term contribution stress -= zeta * n * std::pow(J, -n) * FinvT; #if AMREX_SPACEDIM == 2 @@ -228,10 +228,10 @@ public: t1 += (1.0 / 3.0) * (trace_FTF - 3.0) * FinvT(i, l) * FinvT(k, j); // Bulk-related term (kappa) t2 += (2.0 * J * J - J) * FinvT(i, j) * FinvT(k, l) - + (J - J * J) * FinvT(i, l) * FinvT(k, j); + + (J - J * J) * FinvT(i, l) * FinvT(k, j); // Adhesion-related term (zeta) t3 += (n * n) * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) - + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); + + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); stiffness(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; } } From 053f43fa8506279d7ae8f3c3418b9c706e0ea938 Mon Sep 17 00:00:00 2001 From: Brandon Runnels Date: Mon, 14 Apr 2025 14:29:40 -0500 Subject: [PATCH 13/18] Final tweaks --- scripts/runtests.py | 2 +- tests/Solid/input | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/scripts/runtests.py b/scripts/runtests.py index 9272c0a4be..d28c12573d 100755 --- a/scripts/runtests.py +++ b/scripts/runtests.py @@ -629,4 +629,4 @@ class stats: return_code += stats.timeouts # Return nonzero only if no tests failed or were unexpectedly skipped -exit(return_code) \ No newline at end of file +exit(return_code) diff --git a/tests/Solid/input b/tests/Solid/input index cb86578f7f..969d55bdee 100644 --- a/tests/Solid/input +++ b/tests/Solid/input @@ -125,10 +125,6 @@ #@ benchmark-github = 11.35 #@ - - - - alamo.program = mechanics @@ -137,7 +133,7 @@ plot_file = tests/Solid/output type=static timestep = 0.01 -stop_time = 1 +stop_time = 1.0 # amr parameters amr.plot_dt = 0.1 From ef83cbf5b70065413ae5596832cb6cafbe75e7f0 Mon Sep 17 00:00:00 2001 From: Brandon Runnels Date: Mon, 14 Apr 2025 14:36:54 -0500 Subject: [PATCH 14/18] removing extra code and bringing into compliance --- src/Model/Solid/Finite/Adhesion.H | 49 +++++++------------------------ 1 file changed, 10 insertions(+), 39 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 82a82b4c44..511f81596a 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -40,37 +40,6 @@ // 2. Non-degradable (adhesion term decaying with volume expansion) // This model is rate-independent, reversible, and purely elastic. -namespace Set -{ - -// Outer product for two 3x3 matrices producing a 3x3x3x3 tensor -inline std::array, 3>, 3>, 3> -Outer(const Eigen::Matrix3d &A, const Eigen::Matrix3d &B) -{ - std::array, 3>, 3>, 3> result = {}; - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - for (int k = 0; k < 3; ++k) - for (int l = 0; l < 3; ++l) - result[i][j][k][l] = A(i, j) * B(k, l); - return result; -} - -// Derivative of inverse transpose: d(F^-T) / dF -inline std::array, 3>, 3>, 3> -Derivative(const Eigen::Matrix3d &FinvT) -{ - std::array, 3>, 3>, 3> result = {}; - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - for (int k = 0; k < 3; ++k) - for (int l = 0; l < 3; ++l) - result[i][j][k][l] = -FinvT(i, j) * FinvT(j, i); - return result; -} - -} // namespace Set - namespace Model { namespace Solid @@ -123,11 +92,16 @@ public: static void Parse(Adhesion &value, IO::ParmParse &pp) { - pp_query("mu", value.mu); ///< Shear modulus - pp_query("kappa", value.kappa); ///< Bulk modulus - pp_query("d", value.d); ///< Adhesion strength (degradable part coefficient) - pp_query("zeta", value.zeta); ///< Adhesion decay parameter (non-degradable part) - pp_query("n", value.n); ///< Decay exponent + // Shear modulus + pp_query_required("mu", value.mu); + // Bulk modulus + pp_query_required("kappa", value.kappa); + // Initial value for damage + pp_query_default("d", value.d, 1.0); + // Adhesion decay parameter (non-degradable part) + pp_query_default("zeta", value.zeta, 0.0); + // Decay exponent + pp_query_default("n", value.n,1.0); } // Strain energy function: @@ -205,9 +179,6 @@ public: Eigen::Matrix3d FinvT = F.inverse().transpose(); Set::Scalar trace_FTF = (F.transpose() * F).trace(); - // Precompute derivative helper for FinvT - auto Deriv = ::Set::Derivative(FinvT); - Set::Matrix4<3, Set::Sym::Major> stiffness; for (int i = 0; i < 3; ++i) { From 0a436fecaaa1cb450aa03f64b8268aa0672532f0 Mon Sep 17 00:00:00 2001 From: Brandon Runnels Date: Mon, 14 Apr 2025 16:07:58 -0500 Subject: [PATCH 15/18] integrating with nh --- input | 53 +++++++++ src/Model/Solid/Finite/Adhesion.H | 176 ++++++++++-------------------- 2 files changed, 110 insertions(+), 119 deletions(-) create mode 100644 input diff --git a/input b/input new file mode 100644 index 0000000000..9619e4bd5f --- /dev/null +++ b/input @@ -0,0 +1,53 @@ +alamo.program = mechanics +alamo.program.mechanics.model = finite.adhesion + +plot_file = output + +# this is not a time integration, so do +# exactly one timestep and then quit +timestep = 0.1 +stop_time = 1.0 + +# amr parameters +amr.plot_int = 1 +amr.max_level = 5 +amr.n_cell = 32 32 32 +amr.blocking_factor = 8 +amr.regrid_int = 1 +amr.grid_eff = 1.0 +amr.cell.all = 1 + +# geometry +geometry.prob_lo = -8 -8 -8 +geometry.prob_hi = 8 8 8 +geometry.is_periodic = 0 0 0 + +#ic.type = ellipse +#ic.ellipse.a = 1.0 0.75 0.5 # ellipse radii +#ic.ellipse.x0 = 0 0 0 # location of ellipse center +#ic.ellipse.eps = 0.1 # diffuse boundary + +ic.type = expression +ic.expression.constant.eps = 0.1 +ic.expression.region0 = "0.5 + 0.5*erf((x^2 + y^2 - 1.0)/eps)" +ic.expression.region1 = "1 - (0.5 + 0.5*erf((x^2 + y^2 - 1.0)/eps))" +#ic.expression.region2 = "0.0" + +# elastic moduli +nmodels = 2 +model1.mu = 3.0 +model1.kappa = 6.5 +model2.mu = 6.0 +model2.kappa = 0.3 +#model3.mu = 210 +#model3.kappa = 0.3 + +solver.verbose = 3 +solver.nriters = 1 +#solver.max_iter = 20 +#solver.fixed_iter = 10 +#elasticop.small = 0.1 + +bc.type = tensiontest +bc.tensiontest.type = uniaxial_stress +bc.tensiontest.disp=(0,1:0,0.1) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 511f81596a..95b28fbdae 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -1,13 +1,3 @@ -#ifndef MODEL_SOLID_FINITE_ADHESION_H_ -#define MODEL_SOLID_FINITE_ADHESION_H_ - -#include "IO/ParmParse.H" -#include "Model/Solid/Solid.H" -#include "Set/Set.H" -#include // For std::array to represent fourth-order tensors -#include // For std::pow and std::abs -#include // For Eigen::Matrix3d - // This models a finite-strain, solid with adhesion behavior governed by // a combination of elastic energy and an adhesion decay term. // @@ -40,6 +30,17 @@ // 2. Non-degradable (adhesion term decaying with volume expansion) // This model is rate-independent, reversible, and purely elastic. +#ifndef MODEL_SOLID_FINITE_ADHESION_H_ +#define MODEL_SOLID_FINITE_ADHESION_H_ + +#include "IO/ParmParse.H" +#include "Model/Solid/Finite/NeoHookean.H" +#include "Model/Solid/Solid.H" +#include "Set/Set.H" +#include // For std::array to represent fourth-order tensors +#include // For std::pow and std::abs +#include // For Eigen::Matrix3d + namespace Model { namespace Solid @@ -47,30 +48,30 @@ namespace Solid namespace Finite { -class Adhesion : public Solid +class Adhesion : public NeoHookean { public: + static constexpr KinematicVariable kinvar = KinematicVariable::F; + Adhesion() {} - Adhesion(Solid base) : Solid(base) {} + Adhesion(Solid base) : NeoHookean(base) {} virtual ~Adhesion() {} // Material parameters (declared in order) - Set::Scalar d = NAN; ///< Adhesion strength (degradable part coefficient) - Set::Scalar mu = NAN; ///< Shear modulus - Set::Scalar kappa = NAN; ///< Bulk modulus - Set::Scalar zeta = NAN; ///< Adhesion decay parameter (non-degradable part) - Set::Scalar n = NAN; ///< Decay exponent + Set::Scalar d = NAN; ///< Adhesion strength (degradable part coefficient) + Set::Scalar zeta = NAN; ///< Adhesion decay parameter (non-degradable part) + Set::Scalar n = NAN; ///< Decay exponent // Return a zero-parameter adhesion model. static Adhesion Zero() { Adhesion adhesion; - adhesion.d = 0; ///< Adhesion strength (degradable part coefficient) - adhesion.mu = 0; ///< Shear modulus - adhesion.kappa = 0; ///< Bulk modulus - adhesion.zeta = 0; ///< Adhesion decay parameter (non-degradable part) - adhesion.n = 0; ///< Decay exponent + adhesion.d = 0.0; ///< Adhesion strength (degradable part coefficient) + adhesion.zeta = 0.0; ///< Adhesion decay parameter (non-degradable part) + adhesion.n = 0.0; ///< Decay exponent + adhesion.mu = 0.0; + adhesion.kappa = 0.0; return adhesion; } @@ -79,11 +80,11 @@ public: Random() { Adhesion adhesion; - adhesion.d = Util::Random(); ///< Adhesion strength (degradable part coefficient) - adhesion.mu = Util::Random(); ///< Shear modulus - adhesion.kappa = Util::Random(); ///< Bulk modulus - adhesion.zeta = Util::Random(); ///< Adhesion decay parameter (non-degradable part) - adhesion.n = Util::Random(); ///< Decay exponent + adhesion.d =1.0 ; // Util::Random(); ///< Adhesion strength (degradable part coefficient) + adhesion.zeta = 1.0; ///< Adhesion decay parameter (non-degradable part) + adhesion.n = Util::Random(); ///< Decay exponent + adhesion.mu = Util::Random(); + adhesion.kappa = Util::Random();// Util::Random(); return adhesion; } @@ -92,133 +93,70 @@ public: static void Parse(Adhesion &value, IO::ParmParse &pp) { - // Shear modulus - pp_query_required("mu", value.mu); - // Bulk modulus - pp_query_required("kappa", value.kappa); + pp.queryclass(value); // Initial value for damage pp_query_default("d", value.d, 1.0); // Adhesion decay parameter (non-degradable part) pp_query_default("zeta", value.zeta, 0.0); // Decay exponent - pp_query_default("n", value.n,1.0); + pp_query_default("n", value.n, 1.0); } // Strain energy function: // W = d * [ 0.5*mu*(trace(F^T F)-3)/J^(2/3) + 0.5*kappa*(J-1)^2 ] + zeta/J^n Set::Scalar - W(const Set::Matrix &a_F) const override + W(const Set::Matrix &F) const override { -#if AMREX_SPACEDIM == 2 - Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); - F(0, 0) = a_F(0, 0); - F(0, 1) = a_F(0, 1); - F(1, 0) = a_F(1, 0); - F(1, 1) = a_F(1, 1); -#elif AMREX_SPACEDIM == 3 - Eigen::Matrix3d F = a_F; -#endif + Set::Scalar w = (1.0 - d) * NeoHookean::W(F); Set::Scalar J = F.determinant(); - Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); - Set::Scalar energy = 0.0; - // Degradable part: shear and volumetric terms (scaled by d) - energy += d * (0.5 * mu * ((F.transpose() * F).trace() - 3.0) / J23 + 0.5 * kappa * (J - 1.0) * (J - 1.0)); // Non-degradable (adhesion) part: decays with volume expansion - energy += zeta / std::pow(J, n); - return energy; + w += zeta * std::pow(J, -n); + return w; } // First Piola-Kirchhoff stress (derivative of W) Set::Matrix - DW(const Set::Matrix &a_F) const override + DW(const Set::Matrix &F) const override { -#if AMREX_SPACEDIM == 2 - Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); - F(0, 0) = a_F(0, 0); - F(0, 1) = a_F(0, 1); - F(1, 0) = a_F(1, 0); - F(1, 1) = a_F(1, 1); -#elif AMREX_SPACEDIM == 3 - Eigen::Matrix3d F = a_F; -#endif + Set::Matrix dw = (1.0 - d) * NeoHookean::DW(F); Set::Scalar J = F.determinant(); - Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); - Eigen::Matrix3d FinvT = F.inverse().transpose(); - // Derivative from degradable elastic part and bulk contribution - Eigen::Matrix3d stress = d * mu * (F / J23 - (1.0 / 3.0) * (((F.transpose() * F).trace()) - 3.0) * FinvT / J23) - + d * kappa * (J - 1.0) * J * FinvT; - // Adhesion term contribution - stress -= zeta * n * std::pow(J, -n) * FinvT; -#if AMREX_SPACEDIM == 2 - Set::Matrix r_stress; - r_stress(0, 0) = stress(0, 0); - r_stress(0, 1) = stress(0, 1); - r_stress(1, 0) = stress(1, 0); - r_stress(1, 1) = stress(1, 1); - return r_stress; -#elif AMREX_SPACEDIM == 3 - return stress; -#endif + Set::Matrix FinvT = F.inverse().transpose(); + dw += zeta * (-n * std::pow(J, -n) * FinvT); + return dw; } // Tangent stiffness tensor (second derivative of W) Set::Matrix4 - DDW(const Set::Matrix &a_F) const override + DDW(const Set::Matrix &F) const override { -#if AMREX_SPACEDIM == 2 - Eigen::Matrix3d F = Eigen::Matrix3d::Identity(); - F(0, 0) = a_F(0, 0); - F(0, 1) = a_F(0, 1); - F(1, 0) = a_F(1, 0); - F(1, 1) = a_F(1, 1); -#elif AMREX_SPACEDIM == 3 - Eigen::Matrix3d F = a_F; -#endif + Set::Matrix4 + ddw = (1.0 - d) * NeoHookean::DDW(F); + Set::Scalar J = F.determinant(); - Set::Scalar J23 = std::pow(std::abs(J), 2.0 / 3.0); - Eigen::Matrix3d FinvT = F.inverse().transpose(); - Set::Scalar trace_FTF = (F.transpose() * F).trace(); + Set::Matrix FinvT = F.inverse().transpose(); - Set::Matrix4<3, Set::Sym::Major> stiffness; - for (int i = 0; i < 3; ++i) + for (int i = 0; i < AMREX_SPACEDIM; i++) { - for (int j = 0; j < 3; ++j) + for (int j = 0; j < AMREX_SPACEDIM; j++) { - for (int k = 0; k < 3; ++k) + for (int k = 0; k < AMREX_SPACEDIM; k++) { - for (int l = 0; l < 3; ++l) + for (int l = 0; l < AMREX_SPACEDIM; l++) { - stiffness(i, j, k, l) = 0.0; - Set::Scalar t1 = 0.0, t2 = 0.0, t3 = 0.0; - // Shear-related term (mu) - if (i == k && j == l) - t1 += 1.0; - t1 -= (2.0 / 3.0) * F(i, j) * FinvT(k, l); - t1 -= (2.0 / 3.0) * FinvT(i, j) * F(k, l); - t1 += (2.0 / 9.0) * (trace_FTF - 3.0) * FinvT(i, j) * FinvT(k, l); - t1 += (1.0 / 3.0) * (trace_FTF - 3.0) * FinvT(i, l) * FinvT(k, j); - // Bulk-related term (kappa) - t2 += (2.0 * J * J - J) * FinvT(i, j) * FinvT(k, l) - + (J - J * J) * FinvT(i, l) * FinvT(k, j); // Adhesion-related term (zeta) - t3 += (n * n) * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) - + n * std::pow(J, -n) * FinvT(i, l) * FinvT(k, j); - stiffness(i, j, k, l) = d * mu * t1 / J23 + d * kappa * t2 + zeta * t3; + + // dFinvT(i,j) / dF(k,l) + // = dFinv(j,i) / dF(k,l) + // = -dFinvT(k,j) * dFinvT(iL) + + Set::Scalar t3 = n * n * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) + + n * std::pow(J, -n) * FinvT(k, j) * FinvT(i, l); + ddw(i, j, k, l) += zeta * t3; } } } } -#if AMREX_SPACEDIM == 2 - Set::Matrix4 r_stiffness; - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - for (int k = 0; k < 2; k++) - for (int l = 0; l < 2; l++) - r_stiffness(i, j, k, l) = stiffness(i, j, k, l); - return r_stiffness; -#elif AMREX_SPACEDIM == 3 - return stiffness; -#endif + return ddw; } // Print material parameters From 5c7b42ab3d093a49b2042663f4871129f2c8f7f9 Mon Sep 17 00:00:00 2001 From: Brandon Runnels Date: Mon, 14 Apr 2025 16:43:27 -0500 Subject: [PATCH 16/18] fixed tests (issue was failing to apply symmetry in for loops for matrix4 --- src/Model/Solid/Finite/Adhesion.H | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 95b28fbdae..d9ffb51b12 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -80,8 +80,8 @@ public: Random() { Adhesion adhesion; - adhesion.d =1.0 ; // Util::Random(); ///< Adhesion strength (degradable part coefficient) - adhesion.zeta = 1.0; ///< Adhesion decay parameter (non-degradable part) + adhesion.d = Util::Random() ; // Util::Random(); ///< Adhesion strength (degradable part coefficient) + adhesion.zeta = Util::Random(); ///< Adhesion decay parameter (non-degradable part) adhesion.n = Util::Random(); ///< Decay exponent adhesion.mu = Util::Random(); adhesion.kappa = Util::Random();// Util::Random(); @@ -109,7 +109,6 @@ public: { Set::Scalar w = (1.0 - d) * NeoHookean::W(F); Set::Scalar J = F.determinant(); - // Non-degradable (adhesion) part: decays with volume expansion w += zeta * std::pow(J, -n); return w; } @@ -129,12 +128,18 @@ public: Set::Matrix4 DDW(const Set::Matrix &F) const override { - Set::Matrix4 - ddw = (1.0 - d) * NeoHookean::DDW(F); + Set::Matrix4 ddw = + (1.0 - d) * NeoHookean::DDW(F); + + Util::Message(INFO,zeta); + Util::Message(INFO,n); + Util::Message(INFO,d); Set::Scalar J = F.determinant(); Set::Matrix FinvT = F.inverse().transpose(); + Set::Matrix4 ddw_compress; + for (int i = 0; i < AMREX_SPACEDIM; i++) { for (int j = 0; j < AMREX_SPACEDIM; j++) @@ -151,12 +156,12 @@ public: Set::Scalar t3 = n * n * std::pow(J, -n) * FinvT(i, j) * FinvT(k, l) + n * std::pow(J, -n) * FinvT(k, j) * FinvT(i, l); - ddw(i, j, k, l) += zeta * t3; + ddw_compress(i, j, k, l) = zeta * t3; } } } } - return ddw; + return ddw + ddw_compress; } // Print material parameters From f1f492084b95bffeb28c9df53994241097a975f0 Mon Sep 17 00:00:00 2001 From: Brandon Runnels Date: Mon, 14 Apr 2025 17:12:59 -0500 Subject: [PATCH 17/18] working test case --- input | 29 +++++++++++++++++++---------- src/Model/Solid/Finite/Adhesion.H | 20 ++++++++------------ 2 files changed, 27 insertions(+), 22 deletions(-) diff --git a/input b/input index 9619e4bd5f..8ab7c85b9a 100644 --- a/input +++ b/input @@ -10,7 +10,7 @@ stop_time = 1.0 # amr parameters amr.plot_int = 1 -amr.max_level = 5 +amr.max_level = 6 amr.n_cell = 32 32 32 amr.blocking_factor = 8 amr.regrid_int = 1 @@ -28,26 +28,35 @@ geometry.is_periodic = 0 0 0 #ic.ellipse.eps = 0.1 # diffuse boundary ic.type = expression -ic.expression.constant.eps = 0.1 -ic.expression.region0 = "0.5 + 0.5*erf((x^2 + y^2 - 1.0)/eps)" -ic.expression.region1 = "1 - (0.5 + 0.5*erf((x^2 + y^2 - 1.0)/eps))" +ic.expression.constant.eps = 0.3 +ic.expression.region0 = "(max(erf((x^2 + y^2 - 1.0)/eps),0.0))^2" +ic.expression.region1 = "(max(-erf((x^2 + y^2 - 1.0)/eps),0.0))^2" +ic.expression.region2 = "1.0 - (max(erf((x^2 + y^2 - 1.0)/eps),0.0))^2 - (max(-erf((x^2 + y^2 - 1.0)/eps),0.0))^2" #ic.expression.region2 = "0.0" # elastic moduli nmodels = 2 -model1.mu = 3.0 -model1.kappa = 6.5 -model2.mu = 6.0 -model2.kappa = 0.3 +model1.mu = 6.0 +model1.kappa = 6.0 +model2.mu = 20.0 +model2.kappa = 20 +model3.mu = 6.0 +model3.kappa = 6 +model1.d = 0.0 +model2.d = 0.0 +model3.d = 0.6 +model1.zeta = 0.0 +model2.zeta = 0.0 +model3.zeta = 0.0 #model3.mu = 210 #model3.kappa = 0.3 solver.verbose = 3 solver.nriters = 1 #solver.max_iter = 20 -#solver.fixed_iter = 10 +solver.fixed_iter = 100 #elasticop.small = 0.1 bc.type = tensiontest bc.tensiontest.type = uniaxial_stress -bc.tensiontest.disp=(0,1:0,0.1) +bc.tensiontest.disp=(0,1:0,8) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index d9ffb51b12..858dfcc29b 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -107,34 +107,30 @@ public: Set::Scalar W(const Set::Matrix &F) const override { - Set::Scalar w = (1.0 - d) * NeoHookean::W(F); + Set::Scalar w_nh = (1.0 - d) * NeoHookean::W(F); Set::Scalar J = F.determinant(); - w += zeta * std::pow(J, -n); - return w; + Set::Scalar w_compress = zeta * std::pow(J, -n); + return w_nh + w_compress; } // First Piola-Kirchhoff stress (derivative of W) Set::Matrix DW(const Set::Matrix &F) const override { - Set::Matrix dw = (1.0 - d) * NeoHookean::DW(F); + Set::Matrix dw_nh = (1.0 - d) * NeoHookean::DW(F); Set::Scalar J = F.determinant(); Set::Matrix FinvT = F.inverse().transpose(); - dw += zeta * (-n * std::pow(J, -n) * FinvT); - return dw; + Set::Matrix dw_compress = zeta * (-n * std::pow(J, -n) * FinvT); + return dw_nh + dw_compress; } // Tangent stiffness tensor (second derivative of W) Set::Matrix4 DDW(const Set::Matrix &F) const override { - Set::Matrix4 ddw = + Set::Matrix4 ddw_nh = (1.0 - d) * NeoHookean::DDW(F); - Util::Message(INFO,zeta); - Util::Message(INFO,n); - Util::Message(INFO,d); - Set::Scalar J = F.determinant(); Set::Matrix FinvT = F.inverse().transpose(); @@ -161,7 +157,7 @@ public: } } } - return ddw + ddw_compress; + return ddw_nh + ddw_compress; } // Print material parameters From 3e0064247666360e9d13b028d816995bf0c4f386 Mon Sep 17 00:00:00 2001 From: Brandon Runnels Date: Thu, 24 Apr 2025 17:00:40 -0500 Subject: [PATCH 18/18] working example --- input | 37 ++++++++++-------- src/Model/Solid/Finite/Adhesion.H | 64 ++++++++++++++++++++++++++++++- 2 files changed, 84 insertions(+), 17 deletions(-) diff --git a/input b/input index 8ab7c85b9a..5ed36c65f5 100644 --- a/input +++ b/input @@ -5,8 +5,8 @@ plot_file = output # this is not a time integration, so do # exactly one timestep and then quit -timestep = 0.1 -stop_time = 1.0 +timestep = 0.01 +stop_time = 2.0 # amr parameters amr.plot_int = 1 @@ -16,6 +16,7 @@ amr.blocking_factor = 8 amr.regrid_int = 1 amr.grid_eff = 1.0 amr.cell.all = 1 +amr.thermo.plot_int = 1 # geometry geometry.prob_lo = -8 -8 -8 @@ -29,34 +30,40 @@ geometry.is_periodic = 0 0 0 ic.type = expression ic.expression.constant.eps = 0.3 -ic.expression.region0 = "(max(erf((x^2 + y^2 - 1.0)/eps),0.0))^2" -ic.expression.region1 = "(max(-erf((x^2 + y^2 - 1.0)/eps),0.0))^2" -ic.expression.region2 = "1.0 - (max(erf((x^2 + y^2 - 1.0)/eps),0.0))^2 - (max(-erf((x^2 + y^2 - 1.0)/eps),0.0))^2" +ic.expression.constant.R = 4.0 +ic.expression.region0 = " (max( erf((x^2 + y^2 - R)/eps),0.0))^2" +ic.expression.region1 = " (max(-erf((x^2 + y^2 - R)/eps),0.0))^2" +ic.expression.region2 = "1.0 - (max( erf((x^2 + y^2 - R)/eps),0.0))^2 - (max(-erf((x^2 + y^2 - R)/eps),0.0))^2" #ic.expression.region2 = "0.0" # elastic moduli nmodels = 2 -model1.mu = 6.0 -model1.kappa = 6.0 -model2.mu = 20.0 -model2.kappa = 20 -model3.mu = 6.0 -model3.kappa = 6 +model1.mu = 2.0 +model1.kappa = 3.0 +model2.mu = 10.0 +model2.kappa = 30.0 +model3.mu = 2.0 +model3.kappa = 3.0 model1.d = 0.0 model2.d = 0.0 -model3.d = 0.6 +model3.d = 0.0 model1.zeta = 0.0 model2.zeta = 0.0 -model3.zeta = 0.0 +model3.zeta = 1.0 +model1.L = 0.0 +model2.L = 0.0 +model3.L = 0.2 #model3.mu = 210 #model3.kappa = 0.3 solver.verbose = 3 -solver.nriters = 1 +solver.nriters = 10 #solver.max_iter = 20 solver.fixed_iter = 100 #elasticop.small = 0.1 +print_model = 1 + bc.type = tensiontest bc.tensiontest.type = uniaxial_stress -bc.tensiontest.disp=(0,1:0,8) +bc.tensiontest.disp=(0,1,2:0,8,0) diff --git a/src/Model/Solid/Finite/Adhesion.H b/src/Model/Solid/Finite/Adhesion.H index 858dfcc29b..e8c7788642 100644 --- a/src/Model/Solid/Finite/Adhesion.H +++ b/src/Model/Solid/Finite/Adhesion.H @@ -61,6 +61,7 @@ public: Set::Scalar d = NAN; ///< Adhesion strength (degradable part coefficient) Set::Scalar zeta = NAN; ///< Adhesion decay parameter (non-degradable part) Set::Scalar n = NAN; ///< Decay exponent + Set::Scalar L = NAN; ///< Mobility // Return a zero-parameter adhesion model. static Adhesion @@ -70,6 +71,7 @@ public: adhesion.d = 0.0; ///< Adhesion strength (degradable part coefficient) adhesion.zeta = 0.0; ///< Adhesion decay parameter (non-degradable part) adhesion.n = 0.0; ///< Decay exponent + adhesion.L = 0.0; adhesion.mu = 0.0; adhesion.kappa = 0.0; return adhesion; @@ -85,6 +87,7 @@ public: adhesion.n = Util::Random(); ///< Decay exponent adhesion.mu = Util::Random(); adhesion.kappa = Util::Random();// Util::Random(); + adhesion.L = Util::Random(); return adhesion; } @@ -95,11 +98,13 @@ public: { pp.queryclass(value); // Initial value for damage - pp_query_default("d", value.d, 1.0); + pp_query_default("d", value.d, 0.0); // Adhesion decay parameter (non-degradable part) pp_query_default("zeta", value.zeta, 0.0); // Decay exponent pp_query_default("n", value.n, 1.0); + // Mobility + pp_query_default("L", value.L, 0.0); } // Strain energy function: @@ -168,8 +173,21 @@ public: << ", zeta = " << zeta << ", n = " << n; } + virtual void Advance(Set::Scalar dt, Set::Matrix strain, + Set::Matrix P, Set::Scalar /*time*/) override + { + // Util::Message(INFO,dt); + // Util::Message(INFO,L); + // Util::Message(INFO,(P.transpose()*P).trace()); + Set::Scalar dd = std::max(0.0,(P.transpose()*P).trace() * (0.8 - d)); + // Util::Message(INFO,dd); + d += dd * L * dt; + } + + + #define OP_CLASS Adhesion -#define OP_VARS X(d) X(mu) X(kappa) X(zeta) X(n) +#define OP_VARS X(d) X(mu) X(kappa) X(zeta) X(n) X(L) #include "Model/Solid/InClassOperators.H" }; @@ -179,4 +197,46 @@ public: } // namespace Solid } // namespace Model +template<> +ALAMO_SINGLE_DEFINITION +int Set::Field::NComp() const +{ + return 4; +} +template<> +ALAMO_SINGLE_DEFINITION +std::string Set::Field::Name(int i) const +{ + if (i == 0) return name + "_mu"; + if (i == 1) return name + "_kapp"; + if (i == 2) return name + "_d"; + if (i == 3) return name + "_L"; + return name; +} +template<> +ALAMO_SINGLE_DEFINITION +void Set::Field::Copy(int a_lev, amrex::MultiFab& a_dst, int a_dstcomp, int a_nghost) const +{ + for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost)); + if (bx.ok()) + { + amrex::Array4 const& src = ((*this)[a_lev])->array(mfi); + amrex::Array4 const& dst = a_dst.array(mfi); + for (int n = 0; n < AMREX_SPACEDIM * AMREX_SPACEDIM; n++) + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + dst(i, j, k, a_dstcomp + 0) = src(i, j, k).mu; + dst(i, j, k, a_dstcomp + 1) = src(i, j, k).kappa; + dst(i, j, k, a_dstcomp + 2) = src(i, j, k).d; + dst(i, j, k, a_dstcomp + 3) = src(i, j, k).L; + }); + } + } + } +} + + + #endif // MODEL_SOLID_FINITE_ADHESION_H_