From 9b22f76dbcdcb296524367cca3d51ea0fbb5810e Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Thu, 24 Jul 2025 12:19:38 -0700 Subject: [PATCH 01/11] SmoothMortar --- src/serac/physics/contact/contact_config.hpp | 12 +- .../physics/contact/contact_interaction.cpp | 12 +- src/serac/physics/tests/CMakeLists.txt | 1 + src/serac/physics/tests/contact_patch_2D.cpp | 162 ++++++++++++++++++ tribol | 2 +- 5 files changed, 185 insertions(+), 4 deletions(-) create mode 100644 src/serac/physics/tests/contact_patch_2D.cpp diff --git a/src/serac/physics/contact/contact_config.hpp b/src/serac/physics/contact/contact_config.hpp index 421121712..a05445d4f 100644 --- a/src/serac/physics/contact/contact_config.hpp +++ b/src/serac/physics/contact/contact_config.hpp @@ -19,16 +19,22 @@ namespace serac { */ enum class ContactMethod { - SingleMortar /**< Puso and Laursen 2004 */ + SingleMortar, /**< Puso and Laursen 2004 */ + SmoothMortar /**< Pointwise gap enforced smoothed contact method */ }; + + + + /** * @brief Describes how to enforce the contact constraint equations */ enum class ContactEnforcement { Penalty, /**< Equal penalty applied to all constrained dofs */ - LagrangeMultiplier /**< Solve for exact pressures to satisfy constraints */ + LagrangeMultiplier, /**< Solve for exact pressures to satisfy constraints */ + NotRequired }; /** @@ -65,6 +71,8 @@ struct ContactOptions { /// Penalty parameter (only used when enforcement == ContactEnforcement::Penalty) double penalty = 1.0e3; + double penalty2 = 0.0; + /// The method to use for Jacobian calculations ContactJacobian jacobian = ContactJacobian::Approximate; }; diff --git a/src/serac/physics/contact/contact_interaction.cpp b/src/serac/physics/contact/contact_interaction.cpp index 849106448..68b9b3839 100644 --- a/src/serac/physics/contact/contact_interaction.cpp +++ b/src/serac/physics/contact/contact_interaction.cpp @@ -5,6 +5,8 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include "serac/physics/contact/contact_interaction.hpp" +#include "tribol/common/Parameters.hpp" +// #include "gtest/gtest.h" #ifdef SERAC_USE_TRIBOL @@ -63,6 +65,13 @@ ContactInteraction::ContactInteraction(int interaction_id, const mfem::ParMesh& mfem::FiniteElementSpace::MarkerToList(tdof_markers, inactive_tdofs_); } + if(getContactOptions().method == ContactMethod::SmoothMortar) + { + contact_opts_.enforcement = ContactEnforcement::NotRequired; + tribol::setMfemKinematicConstantPenalty(interaction_id, contact_opts_.penalty, contact_opts_.penalty2); + contact_opts_.type = ContactType::TiedNormal; + } + // set up Tribol to compute exact Jacobian if requested if (getContactOptions().jacobian == ContactJacobian::Exact) { tribol::enableEnzyme(interaction_id, true); @@ -155,7 +164,8 @@ tribol::ContactMethod ContactInteraction::getMethod() const switch (contact_opts_.method) { case ContactMethod::SingleMortar: return tribol::SINGLE_MORTAR; - break; + case ContactMethod::SmoothMortar: + return tribol::SMOOTH_MORTAR; default: SLIC_ERROR_ROOT("Unsupported contact method."); // return something so we don't get an error diff --git a/src/serac/physics/tests/CMakeLists.txt b/src/serac/physics/tests/CMakeLists.txt index afb44edce..9451ab1d6 100644 --- a/src/serac/physics/tests/CMakeLists.txt +++ b/src/serac/physics/tests/CMakeLists.txt @@ -78,6 +78,7 @@ set(physics_parallel_test_sources ) set(physics_parallel_tribol_test_sources contact_patch.cpp + contact_patch_2D.cpp contact_patch_tied.cpp contact_beam.cpp ) diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp new file mode 100644 index 000000000..360edc428 --- /dev/null +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -0,0 +1,162 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include "serac/physics/solid_mechanics_contact.hpp" + +#include +#include +#include +#include +#include +#include "axom/slic/core/SimpleLogger.hpp" +#include +#include "mfem.hpp" + +#include "shared/mesh/MeshBuilder.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" +#include "serac/physics/boundary_conditions/components.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/mesh.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/serac_config.hpp" +#include "serac/infrastructure/application_manager.hpp" + +namespace serac { + +class ContactTest : public testing::TestWithParam> {}; + +TEST_P(ContactTest, patch) +{ + // NOTE: p must be equal to 1 for now + constexpr int p = 1; + constexpr int dim = 2; + + MPI_Barrier(MPI_COMM_WORLD); + + // Create DataStore + std::string name = "contact_patch_" + std::get<2>(GetParam()); + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, name + "_data"); + + // Construct the appropriate dimension mesh and give it to the data store + + auto mesh = std::make_shared(shared::MeshBuilder::Unify({ + shared::MeshBuilder::SquareMesh(1 , 1).translate({0.0, 0.95}).bdrAttribInfo() + .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), + shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); + + mfem::VisItDataCollection visit_dc("contact_patch_visit", &mesh->mfemParMesh()); + + visit_dc.SetPrefixPath("visit_out"); + visit_dc.Save(); + + + + + mesh->addDomainOfBoundaryElements("x0_faces", serac::by_attr(7)); + mesh->addDomainOfBoundaryElements("y0_faces", serac::by_attr(8)); + mesh->addDomainOfBoundaryElements("Ymax_face", serac::by_attr(9)); + + // TODO: investigate performance with Petsc + // #ifdef SERAC_USE_PETSC + // LinearSolverOptions linear_options{ + // .linear_solver = LinearSolver::PetscGMRES, + // .preconditioner = Preconditioner::Petsc, + // .petsc_preconditioner = PetscPCType::HMG, + // .absolute_tol = 1e-16, + // .print_level = 1, + // }; + // #elif defined(MFEM_USE_STRUMPACK) +#ifdef MFEM_USE_STRUMPACK + LinearSolverOptions linear_options{.linear_solver = LinearSolver::Strumpack, .print_level = 0}; +#else + LinearSolverOptions linear_options{}; + SLIC_INFO_ROOT("Contact requires MFEM built with strumpack."); + return; +#endif + + NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1.0e-13, + .absolute_tol = 1.0e-13, + .max_iterations = 20, + .print_level = 1}; + + ContactOptions contact_options{.method = ContactMethod::SmoothMortar, + .enforcement = std::get<0>(GetParam()), + .type = ContactType::Frictionless, + .penalty = 8.0e2, + .jacobian = std::get<1>(GetParam())}; + + SolidMechanicsContact solid_solver(nonlinear_options, linear_options, + solid_mechanics::default_quasistatic_options, name, mesh); + + double K = 10.0; + double G = 0.25; + solid_mechanics::NeoHookean mat{1.0, K, G}; + solid_solver.setMaterial(mat, mesh->entireBody()); + + // Define the function for the initial displacement and boundary condition + // constexpr int dim = 2; + auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.25}}; }; + + // Define a boundary attribute set and specify initial / boundary conditions + solid_solver.setFixedBCs(mesh->domain("x0_faces"), Component::X); + solid_solver.setFixedBCs(mesh->domain("y0_faces"), Component::Y); + solid_solver.setDisplacementBCs(applied_disp_function, mesh->domain("Ymax_face"), Component::Y); + + // Add the contact interaction + solid_solver.addContactInteraction(0, {6}, {5}, contact_options); + + // Finalize the data structures + solid_solver.completeSetup(); + + std::string paraview_name = name + "_paraview"; + solid_solver.outputStateToDisk(paraview_name); + + // Perform the quasi-static solve + double dt = 1.0; + solid_solver.advanceTimestep(dt); + // solid_solver.advanceTimestep(dt); + + // Output the sidre-based plot files + solid_solver.outputStateToDisk(paraview_name); + + // Check the l2 norm of the displacement dofs + auto c = (3.0 * K - 2.0 * G) / (3.0 * K + G); + mfem::VectorFunctionCoefficient elasticity_sol_coeff(2, [c](const mfem::Vector& x, mfem::Vector& u) { + u[0] = 0.25 * 0.01 * c * x[0]; + u[1] = 0.25 * 0.01 * c * x[1]; + // u[2] = -0.5 * 0.01 * x[2]; + }); + mfem::ParFiniteElementSpace elasticity_fes(solid_solver.reactions().space()); + mfem::ParGridFunction elasticity_sol(&elasticity_fes); + elasticity_sol.ProjectCoefficient(elasticity_sol_coeff); + mfem::ParGridFunction approx_error(elasticity_sol); + approx_error -= solid_solver.displacement().gridFunction(); + auto approx_error_l2 = mfem::ParNormlp(approx_error, 2, MPI_COMM_WORLD); + EXPECT_NEAR(0.0, approx_error_l2, 1.0e-3); +} + +INSTANTIATE_TEST_SUITE_P( + tribol, ContactTest, + testing::Values(std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Approximate, "penalty_approxJ"), + std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactJacobian::Approximate, + "lagrange_multiplier_approxJ"), + std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"), + std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactJacobian::Exact, + "lagrange_multiplier_exactJ"))); + +} // namespace serac + +int main(int argc, char* argv[]) +{ + + + + testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/tribol b/tribol index 8351aa4e1..503df36b4 160000 --- a/tribol +++ b/tribol @@ -1 +1 @@ -Subproject commit 8351aa4e104fce00ee25c760a5c3f3c882557d50 +Subproject commit 503df36b4352c2146fcd75a7ff20ba4a7be05395 From f51d6d1899beb3b15876f0714911503b7d6998f9 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Mon, 28 Jul 2025 15:43:42 -0700 Subject: [PATCH 02/11] Patch test --- .../physics/materials/solid_material.hpp | 1 + src/serac/physics/tests/contact_patch_2D.cpp | 21 ++++++++++--------- tribol | 2 +- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index 73508559e..4cebb94f9 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -113,6 +113,7 @@ struct NeoHookean { constexpr auto I = Identity(); auto lambda = K - (2.0 / 3.0) * G; auto B_minus_I = dot(du_dX, transpose(du_dX)) + transpose(du_dX) + du_dX; + std::cout << "du_dx: " << du_dX << std::endl; //Added debugging code auto logJ = log1p(detApIm1(du_dX)); // Kirchoff stress, in form that avoids cancellation error when F is near I diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp index 360edc428..8e26b2c0e 100644 --- a/src/serac/physics/tests/contact_patch_2D.cpp +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -6,6 +6,7 @@ #include "serac/physics/solid_mechanics_contact.hpp" +#include #include #include #include @@ -23,6 +24,9 @@ #include "serac/physics/materials/solid_material.hpp" #include "serac/serac_config.hpp" #include "serac/infrastructure/application_manager.hpp" +#include + + namespace serac { @@ -87,7 +91,7 @@ TEST_P(ContactTest, patch) ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement = std::get<0>(GetParam()), .type = ContactType::Frictionless, - .penalty = 8.0e2, + .penalty = 1.0, .jacobian = std::get<1>(GetParam())}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, @@ -100,7 +104,7 @@ TEST_P(ContactTest, patch) // Define the function for the initial displacement and boundary condition // constexpr int dim = 2; - auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.25}}; }; + auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.01}}; }; // Define a boundary attribute set and specify initial / boundary conditions solid_solver.setFixedBCs(mesh->domain("x0_faces"), Component::X); @@ -137,24 +141,21 @@ TEST_P(ContactTest, patch) mfem::ParGridFunction approx_error(elasticity_sol); approx_error -= solid_solver.displacement().gridFunction(); auto approx_error_l2 = mfem::ParNormlp(approx_error, 2, MPI_COMM_WORLD); - EXPECT_NEAR(0.0, approx_error_l2, 1.0e-3); + EXPECT_NEAR(0.0, approx_error_l2, 1.0e-2); } INSTANTIATE_TEST_SUITE_P( tribol, ContactTest, - testing::Values(std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Approximate, "penalty_approxJ"), - std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactJacobian::Approximate, - "lagrange_multiplier_approxJ"), - std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"), - std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactJacobian::Exact, - "lagrange_multiplier_exactJ"))); + testing::Values( + std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ") +)); } // namespace serac int main(int argc, char* argv[]) { - +feenableexcept(FE_INVALID | FE_OVERFLOW); testing::InitGoogleTest(&argc, argv); serac::ApplicationManager applicationManager(argc, argv); diff --git a/tribol b/tribol index 503df36b4..1a6130b3c 160000 --- a/tribol +++ b/tribol @@ -1 +1 @@ -Subproject commit 503df36b4352c2146fcd75a7ff20ba4a7be05395 +Subproject commit 1a6130b3ced91380e670d428c4492e98530f036e From ab8ff101a0dcc3cd5d9624e98b57e1dd946ecff8 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Wed, 30 Jul 2025 16:44:13 -0700 Subject: [PATCH 03/11] patch test --- src/serac/numerics/equation_solver.cpp | 1 + src/serac/physics/materials/solid_material.hpp | 1 - src/serac/physics/tests/contact_patch_2D.cpp | 18 ++++++++++++------ 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp index 0af69dcff..4c51a66bb 100644 --- a/src/serac/numerics/equation_solver.cpp +++ b/src/serac/numerics/equation_solver.cpp @@ -75,6 +75,7 @@ class NewtonSolver : public mfem::NewtonSolver { { SERAC_MARK_FUNCTION; prec->Mult(r_, c_); // c = [DF(x_i)]^{-1} [F(x_i)-b] + std::cout << "[DEBUG] norm c: " << c_.Norml2() << " [DEBUG] norm r: " << r_.Norml2() << std::endl; } /// @overload diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index 4cebb94f9..73508559e 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -113,7 +113,6 @@ struct NeoHookean { constexpr auto I = Identity(); auto lambda = K - (2.0 / 3.0) * G; auto B_minus_I = dot(du_dX, transpose(du_dX)) + transpose(du_dX) + du_dX; - std::cout << "du_dx: " << du_dX << std::endl; //Added debugging code auto logJ = log1p(detApIm1(du_dX)); // Kirchoff stress, in form that avoids cancellation error when F is near I diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp index 8e26b2c0e..29e29b7fc 100644 --- a/src/serac/physics/tests/contact_patch_2D.cpp +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -48,7 +48,7 @@ TEST_P(ContactTest, patch) // Construct the appropriate dimension mesh and give it to the data store auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(1 , 1).translate({0.0, 0.95}).bdrAttribInfo() + shared::MeshBuilder::SquareMesh(1 , 1).translate({0.5, 1.0}).bdrAttribInfo() .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); @@ -91,7 +91,8 @@ TEST_P(ContactTest, patch) ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement = std::get<0>(GetParam()), .type = ContactType::Frictionless, - .penalty = 1.0, + .penalty = 1e-2, + .penalty2 = 1e-2, .jacobian = std::get<1>(GetParam())}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, @@ -104,7 +105,7 @@ TEST_P(ContactTest, patch) // Define the function for the initial displacement and boundary condition // constexpr int dim = 2; - auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.01}}; }; + auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.1}}; }; // Define a boundary attribute set and specify initial / boundary conditions solid_solver.setFixedBCs(mesh->domain("x0_faces"), Component::X); @@ -123,19 +124,24 @@ TEST_P(ContactTest, patch) // Perform the quasi-static solve double dt = 1.0; solid_solver.advanceTimestep(dt); + solid_solver.advanceTimestep(dt); + // solid_solver.advanceTimestep(dt); // solid_solver.advanceTimestep(dt); // Output the sidre-based plot files solid_solver.outputStateToDisk(paraview_name); // Check the l2 norm of the displacement dofs + // auto c = 1.0; auto c = (3.0 * K - 2.0 * G) / (3.0 * K + G); mfem::VectorFunctionCoefficient elasticity_sol_coeff(2, [c](const mfem::Vector& x, mfem::Vector& u) { - u[0] = 0.25 * 0.01 * c * x[0]; - u[1] = 0.25 * 0.01 * c * x[1]; + // u[0] = 0.0; + // u[1] = -0.01 * c * x[1]; + u[0] = c * 0.1 * x[0]; + u[1] = -0.05 * x[1]; // u[2] = -0.5 * 0.01 * x[2]; }); - mfem::ParFiniteElementSpace elasticity_fes(solid_solver.reactions().space()); + mfem::ParFiniteElementSpace elasticity_fes(solid_solver.displacement().space()); mfem::ParGridFunction elasticity_sol(&elasticity_fes); elasticity_sol.ProjectCoefficient(elasticity_sol_coeff); mfem::ParGridFunction approx_error(elasticity_sol); From 773ef0f28b90c328ec659d07d98b5a3b427c5efc Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Fri, 1 Aug 2025 16:43:26 -0700 Subject: [PATCH 04/11] Tribol FD Tests --- src/serac/numerics/equation_solver.cpp | 4 + .../physics/tests/contact_finite_diff.cpp | 80 +++++++++++-------- src/serac/physics/tests/contact_patch_2D.cpp | 10 +-- .../physics/tests/tribol_finite_diff.cpp | 61 ++++++++++---- 4 files changed, 100 insertions(+), 55 deletions(-) diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp index 4c51a66bb..884b4531d 100644 --- a/src/serac/numerics/equation_solver.cpp +++ b/src/serac/numerics/equation_solver.cpp @@ -193,6 +193,10 @@ class NewtonSolver : public mfem::NewtonSolver { final_iter = it; final_norm = norm; + // for(int i = 0; i < 16; ++i) { + // std::cout << "X:[" << i + 1 << "] = " << x[i] << std::endl; + // } + if (print_options.summary || (!converged && print_options.warnings) || print_options.first_and_last) { mfem::out << "Newton: Number of iterations: " << final_iter << '\n' << " ||r|| = " << final_norm << '\n'; } diff --git a/src/serac/physics/tests/contact_finite_diff.cpp b/src/serac/physics/tests/contact_finite_diff.cpp index 91fce7e4c..23a8bb913 100644 --- a/src/serac/physics/tests/contact_finite_diff.cpp +++ b/src/serac/physics/tests/contact_finite_diff.cpp @@ -32,7 +32,7 @@ TEST_P(ContactFiniteDiff, patch) { // NOTE: p must be equal to 1 for now constexpr int p = 1; - constexpr int dim = 3; + constexpr int dim = 2; constexpr double eps = 1.0e-7; @@ -48,23 +48,23 @@ TEST_P(ContactFiniteDiff, patch) double shift = eps * 10; // clang-format off auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::CubeMesh(1, 1, 1), - shared::MeshBuilder::CubeMesh(1, 1, 1) + shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.999}).translate({shift, 0.0}).bdrAttribInfo() + .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), + shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5) // shift up height of element - .translate({0.0, 0.0, 0.999}) - // shift x and y so the element edges are not overlapping - .translate({shift, shift, 0.0}) - // change the mesh1 boundary attribute from 1 to 7 - .updateBdrAttrib(1, 7) - // change the mesh1 boundary attribute from 6 to 8 - .updateBdrAttrib(6, 8) + + // // shift x and y so the element edges are not overlapping + + // // change the mesh1 boundary attribute from 1 to 7 + // .updateBdrAttrib(1, 7) + // // change the mesh1 boundary attribute from 6 to 8 + // .updateBdrAttrib(6, 8) }), "patch_mesh", 0, 0); // clang-format on - mesh->addDomainOfBoundaryElements("x0_faces", serac::by_attr(5)); - mesh->addDomainOfBoundaryElements("y0_faces", serac::by_attr(2)); - mesh->addDomainOfBoundaryElements("z0_face", serac::by_attr(1)); - mesh->addDomainOfBoundaryElements("zmax_face", serac::by_attr(8)); + mesh->addDomainOfBoundaryElements("x0_faces", serac::by_attr(7)); + mesh->addDomainOfBoundaryElements("y0_faces", serac::by_attr(8)); + mesh->addDomainOfBoundaryElements("ymax_face", serac::by_attr(9)); #ifdef MFEM_USE_STRUMPACK LinearSolverOptions linear_options{.linear_solver = LinearSolver::Strumpack, .print_level = 0}; @@ -81,10 +81,10 @@ TEST_P(ContactFiniteDiff, patch) .max_iterations = 1, .print_level = 1}; - ContactOptions contact_options{.method = ContactMethod::SingleMortar, + ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement = GetParam().first, - .type = ContactType::TiedNormal, - .penalty = 1.0, + .type = ContactType::Frictionless, + .penalty = 0.1, .jacobian = ContactJacobian::Exact}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, @@ -96,21 +96,19 @@ TEST_P(ContactFiniteDiff, patch) solid_mechanics::NeoHookean mat{1.0, K, G}; solid_solver.setMaterial(mat, mesh->entireBody()); - auto nonzero_disp_bc = [](vec3, double) { return vec3{{0.0, 0.0, 0.0}}; }; + auto nonzero_disp_bc = [](vec2, double) { return vec2{{0.0, 0.0}}; }; // Define a boundary attribute set and specify initial / boundary conditions solid_solver.setFixedBCs(mesh->domain("x0_faces"), Component::X); solid_solver.setFixedBCs(mesh->domain("y0_faces"), Component::Y); - solid_solver.setFixedBCs(mesh->domain("z0_face"), Component::Z); - solid_solver.setDisplacementBCs(nonzero_disp_bc, mesh->domain("zmax_face"), Component::Z); + solid_solver.setDisplacementBCs(nonzero_disp_bc, mesh->domain("ymax_face"), Component::Y); // Create a list of vdofs from Domains auto x0_face_dofs = mesh->domain("x0_faces").dof_list(&solid_solver.displacement().space()); auto y0_face_dofs = mesh->domain("y0_faces").dof_list(&solid_solver.displacement().space()); - auto z0_face_dofs = mesh->domain("z0_face").dof_list(&solid_solver.displacement().space()); - auto zmax_face_dofs = mesh->domain("zmax_face").dof_list(&solid_solver.displacement().space()); + auto ymax_face_dofs = mesh->domain("ymax_face").dof_list(&solid_solver.displacement().space()); mfem::Array bc_vdofs(dim * - (x0_face_dofs.Size() + y0_face_dofs.Size() + z0_face_dofs.Size() + zmax_face_dofs.Size())); + (x0_face_dofs.Size() + y0_face_dofs.Size() + ymax_face_dofs.Size())); int dof_ct = 0; for (int i{0}; i < x0_face_dofs.Size(); ++i) { for (int d{0}; d < dim; ++d) { @@ -122,21 +120,21 @@ TEST_P(ContactFiniteDiff, patch) bc_vdofs[dof_ct++] = solid_solver.displacement().space().DofToVDof(y0_face_dofs[i], d); } } - for (int i{0}; i < z0_face_dofs.Size(); ++i) { + // for (int i{0}; i < z0_face_dofs.Size(); ++i) { + // for (int d{0}; d < dim; ++d) { + // bc_vdofs[dof_ct++] = solid_solver.displacement().space().DofToVDof(z0_face_dofs[i], d); + // } + // } + for (int i{0}; i < ymax_face_dofs.Size(); ++i) { for (int d{0}; d < dim; ++d) { - bc_vdofs[dof_ct++] = solid_solver.displacement().space().DofToVDof(z0_face_dofs[i], d); - } - } - for (int i{0}; i < zmax_face_dofs.Size(); ++i) { - for (int d{0}; d < dim; ++d) { - bc_vdofs[dof_ct++] = solid_solver.displacement().space().DofToVDof(zmax_face_dofs[i], d); + bc_vdofs[dof_ct++] = solid_solver.displacement().space().DofToVDof(ymax_face_dofs[i], d); } } bc_vdofs.Sort(); bc_vdofs.Unique(); // Add the contact interaction - solid_solver.addContactInteraction(0, {6}, {7}, contact_options); + solid_solver.addContactInteraction(0, {6}, {5}, contact_options); // Finalize the data structures solid_solver.completeSetup(); @@ -181,6 +179,17 @@ TEST_P(ContactFiniteDiff, patch) J_fd -= f; J_fd /= eps; merged_sol[j] -= eps; + + for(int m = 0; m < 16; ++m) { + std::cout << "J exact: " << J_exact[m] << std:: endl; +} + for(int m = 0; m < 16; ++m) { + std::cout << "J FD: " << J_fd[m] << std:: endl; +} + + + + // loop through forces (row = k) for (int k{0}; k < merged_sol.Size(); ++k) { if (J_exact[k] != 1.0 && (std::abs(J_exact[k]) > 1.0e-15 || std::abs(J_fd[k]) > 1.0e-15)) { @@ -200,12 +209,17 @@ TEST_P(ContactFiniteDiff, patch) std::cout << "Max diff = " << std::setprecision(15) << max_diff << std::endl; solid_solver.advanceTimestep(dt); + + } } + + + INSTANTIATE_TEST_SUITE_P(tribol, ContactFiniteDiff, - testing::Values(std::make_pair(ContactEnforcement::Penalty, "penalty"), - std::make_pair(ContactEnforcement::LagrangeMultiplier, "lm"))); + testing::Values(std::make_pair(ContactEnforcement::Penalty, "penalty"))); + // std::make_pair(ContactEnforcement::LagrangeMultiplier, "lm"))); } // namespace serac diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp index 29e29b7fc..4a7d80b81 100644 --- a/src/serac/physics/tests/contact_patch_2D.cpp +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -48,7 +48,7 @@ TEST_P(ContactTest, patch) // Construct the appropriate dimension mesh and give it to the data store auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(1 , 1).translate({0.5, 1.0}).bdrAttribInfo() + shared::MeshBuilder::SquareMesh(1 , 1).translate({0.0, 1.0}).bdrAttribInfo() .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); @@ -91,8 +91,8 @@ TEST_P(ContactTest, patch) ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement = std::get<0>(GetParam()), .type = ContactType::Frictionless, - .penalty = 1e-2, - .penalty2 = 1e-2, + .penalty = 0.2, + .penalty2 = 0.0, .jacobian = std::get<1>(GetParam())}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, @@ -137,8 +137,8 @@ TEST_P(ContactTest, patch) mfem::VectorFunctionCoefficient elasticity_sol_coeff(2, [c](const mfem::Vector& x, mfem::Vector& u) { // u[0] = 0.0; // u[1] = -0.01 * c * x[1]; - u[0] = c * 0.1 * x[0]; - u[1] = -0.05 * x[1]; + u[0] = c * -0.1 * x[0]; + u[1] = -0.1 * x[1]; // u[2] = -0.5 * 0.01 * x[2]; }); mfem::ParFiniteElementSpace elasticity_fes(solid_solver.displacement().space()); diff --git a/src/serac/physics/tests/tribol_finite_diff.cpp b/src/serac/physics/tests/tribol_finite_diff.cpp index 2d9b41fec..c59237fd6 100644 --- a/src/serac/physics/tests/tribol_finite_diff.cpp +++ b/src/serac/physics/tests/tribol_finite_diff.cpp @@ -41,30 +41,32 @@ TEST_P(TribolFiniteDiff, patch) // Construct the appropriate dimension mesh and give it to the data store + double shift = eps * 10; // clang-format off auto pmesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::CubeMesh(1, 1, 1), - shared::MeshBuilder::CubeMesh(1, 1, 1) - // shift up 99.9% height of element - .translate({0.0, 0.0, 0.999}) - // shift x and y so the element edges are not overlapping - .translate({shift, shift, 0.0}) - // change the mesh1 boundary attribute from 1 to 7 - .updateBdrAttrib(1, 7) - // change the mesh1 boundary attribute from 6 to 8 - .updateBdrAttrib(6, 8) + shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.999}).translate({shift, 0.0}).bdrAttribInfo() + .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), + shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5) + // shift up height of element + + // // shift x and y so the element edges are not overlapping + + // // change the mesh1 boundary attribute from 1 to 7 + // .updateBdrAttrib(1, 7) + // // change the mesh1 boundary attribute from 6 to 8 + // .updateBdrAttrib(6, 8) }), "patch_mesh", 0, 0); // clang-format on - ContactOptions contact_options{.method = ContactMethod::SingleMortar, + ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement = GetParam().first, - .type = ContactType::TiedNormal, - .penalty = 1.0, + .type = ContactType::Frictionless, + .penalty = 0.1, .jacobian = ContactJacobian::Exact}; ContactData contact_data(pmesh->mfemParMesh()); constexpr int interaction_id = 0; - contact_data.addContactInteraction(interaction_id, {6}, {7}, contact_options); + contact_data.addContactInteraction(interaction_id, {6}, {5}, contact_options); mfem::Vector u(pmesh->mfemParMesh().GetNodes()->Size() + contact_data.getContactInteractions()[0].numPressureDofs()); u = 0.0; @@ -76,14 +78,30 @@ TEST_P(TribolFiniteDiff, patch) double max_diff = 0.0; auto J_op = contact_data.mergedJacobian(); + mfem::Vector u_dot(u.Size()); u_dot = 0.0; // wiggle displacement (col = j) - for (int j{0}; j < u.Size(); ++j) { + + + + // // for (int j{0}; j < u.Size(); ++j) { + for (int j{0}; j < 1; ++j) { + int block_size = 8; + std::cout << "entered loop" << std::endl; u_dot[j] = 1.0; mfem::Vector J_exact(u.Size()); J_exact = 0.0; J_op->Mult(u_dot, J_exact); + // Print the i-th entries from the top-left (0,0) block + std::cout << "Column " << j << " of block (0,0):" << std::endl; + for (int i = 0; i < block_size; ++i) { + std::cout << J_exact[i] << " "; + } + std::cout << std::endl; + + + u_dot[j] = 0.0; u[j] += eps; mfem::Vector J_fd(u.Size()); @@ -107,13 +125,22 @@ TEST_P(TribolFiniteDiff, patch) EXPECT_NEAR(J_exact[k], J_fd[k], eps); } } + + for(int m = 0; m < 16; ++m) { + std::cout << "J exact: " << J_exact[m] << std:: endl; +} + for(int m = 0; m < 16; ++m) { + std::cout << "J FD: " << J_fd[m] << std:: endl; +} } + + std::cout << "Max diff = " << std::setprecision(15) << max_diff << std::endl; } INSTANTIATE_TEST_SUITE_P(tribol, TribolFiniteDiff, - testing::Values(std::make_pair(ContactEnforcement::Penalty, "penalty"), - std::make_pair(ContactEnforcement::LagrangeMultiplier, "lm"))); + testing::Values(std::make_pair(ContactEnforcement::Penalty, "penalty") +)); } // namespace serac From 370ca85bb60eea3f304ec92860c797c434bc2815 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Fri, 1 Aug 2025 16:50:51 -0700 Subject: [PATCH 05/11] Tribol FD Tests --- tribol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tribol b/tribol index 1a6130b3c..2ac69027c 160000 --- a/tribol +++ b/tribol @@ -1 +1 @@ -Subproject commit 1a6130b3ced91380e670d428c4492e98530f036e +Subproject commit 2ac69027ce48b47c1f30432cf60907b929f485f4 From 0c08c319c854da6448fe4ba977035e9f5cfb5d1f Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Wed, 13 Aug 2025 09:23:23 -0700 Subject: [PATCH 06/11] ironing test --- examples/contact/ironing_2D.cpp | 177 ++++++++++++++++++ src/serac/numerics/equation_solver.cpp | 1 + .../physics/tests/contact_finite_diff.cpp | 4 +- src/serac/physics/tests/contact_patch_2D.cpp | 36 ++-- .../physics/tests/tribol_finite_diff.cpp | 4 +- 5 files changed, 195 insertions(+), 27 deletions(-) create mode 100644 examples/contact/ironing_2D.cpp diff --git a/examples/contact/ironing_2D.cpp b/examples/contact/ironing_2D.cpp new file mode 100644 index 000000000..73541c2ac --- /dev/null +++ b/examples/contact/ironing_2D.cpp @@ -0,0 +1,177 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include +#include + +#include "axom/slic.hpp" + +#include "mfem.hpp" + +#include "serac/numerics/solver_config.hpp" +#include "serac/physics/contact/contact_config.hpp" +#include "serac/physics/materials/parameterized_solid_material.hpp" +#include "serac/physics/solid_mechanics.hpp" +#include "serac/physics/solid_mechanics_contact.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/serac.hpp" + +#include "serac/physics/contact/contact_config.hpp" +#include "serac/physics/solid_mechanics_contact.hpp" + +#include +#include +#include +#include +#include +#include +#include "axom/slic/core/SimpleLogger.hpp" +#include "mfem.hpp" + +#include "shared/mesh/MeshBuilder.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" +#include "serac/physics/boundary_conditions/components.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/mesh.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/serac_config.hpp" +#include "serac/infrastructure/application_manager.hpp" +#include + + +int main(int argc, char* argv[]) + { + + //Initialize and automatically finalize MPI and other libraries + serac::ApplicationManager applicationManager(argc, argv); + + // NOTE: p must be equal to 1 to work with Tribol's mortar method + constexpr int p = 1; + + //NOTE: dim must be equal to 2 + constexpr int dim = 2; + + //Create DataStore + std::string name = "contact_ironing_2D_example"; + axom::sidre::DataStore datastore; + serac::StateManager::initialize(datastore, name + "_data"); + + //Construct the appropiate dimension mesh and give it to the data store + // std::string filename = SERAC_REPO_DIR "data/meshes/ironing_2D.mesh"; + // std::shared_ptr mesh = std::make_shared(filename, "ironing_2D_mesh", 2, 0); + + auto mesh = std::make_shared(shared::MeshBuilder::Unify({ + shared::MeshBuilder::SquareMesh(32, 32).updateBdrAttrib(1, 6).updateBdrAttrib(3, 9).bdrAttribInfo().scale({1.0, 1.0}), + shared::MeshBuilder::SquareMesh(32, 32).scale({0.25, 0.25}).translate({0.0, 1.0}).updateBdrAttrib(3, 5).updateBdrAttrib(1,8).updateAttrib(1, 2)}), "iroing_2D_mesh", 0, 0); + + serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level=0}; + + mfem::VisItDataCollection visit_dc("contact_ironing_visit", &mesh->mfemParMesh()); + + visit_dc.SetPrefixPath("visit_out"); + visit_dc.Save(); + + #ifndef MFEM_USE_STRUMPACK + SLIC_INFO_ROOT("Contact requires MFEM built with strumpack."); + return 1; + #endif + + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton, + .relative_tol = 1.0e-8, + .absolute_tol = 1.0e-10, + .max_iterations = 50, + .print_level = 1}; + + serac::ContactOptions contact_options{.method = serac::ContactMethod::SmoothMortar, + .enforcement = serac::ContactEnforcement::Penalty, + .type = serac::ContactType::Frictionless, + .penalty = 100, + .penalty2 = 0.0, + .jacobian = serac::ContactJacobian::Exact}; + + serac::SolidMechanicsContact, serac::L2<0>>> solid_solver( + nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, mesh, + {"bulk_mod", "shear_mod"}); + + + serac::FiniteElementState K_field(serac::StateManager::newState(serac::L2<0>{}, "bulk_mod", mesh->tag())); + + mfem::Vector K_values({10.0, 100.0}); + mfem::PWConstCoefficient K_coeff(K_values); + K_field.project(K_coeff); + solid_solver.setParameter(0, K_field); + + serac::FiniteElementState G_field(serac::StateManager::newState(serac::L2<0>{}, "shear_mod", mesh->tag())); + + mfem::Vector G_values({0.25, 2.5}); + mfem::PWConstCoefficient G_coeff(G_values); + G_field.project(G_coeff); + solid_solver.setParameter(1, G_field); + + serac::solid_mechanics::ParameterizedNeoHookeanSolid mat{1.0, 0.0, 0.0}; + solid_solver.setMaterial(serac::DependsOn<0, 1>{}, mat, mesh->entireBody()); + + //Pass the BC information to the solver object + mesh->addDomainOfBoundaryElements("bottom_of_subtrate", serac::by_attr(6)); + solid_solver.setFixedBCs((mesh->domain("bottom_of_subtrate"))); + + mesh->addDomainOfBoundaryElements("top of indenter", serac::by_attr(5)); + auto applied_displacement = [](serac::tensor, double t) { + constexpr double init_steps = 50.0; + serac::tensor u{}; + // std::cout << "T ========= " << t << std::endl; + if (t <= init_steps + 1.0e-12) { + u[1] = -t * 0.2 / init_steps; + // std::cout << "In IF statement. u[1] = " << u[1] << " and t = " << t << std::endl; + } + else { + u[0] = (t - init_steps) * 0.0005; + u[1] = -0.2; + // std::cout << "in ELSE statement. u[1] = " << u[1] << " and u[0] = " << u[0] << " and t = " << t << std::endl; + } + return u; + }; + + solid_solver.setDisplacementBCs(applied_displacement, mesh->domain("top of indenter")); + // std::cout << "top of indenter size: " << mesh->domain("top of indenter").size() << std::endl; + + + //Add the contact interaction + auto contact_interaction_id = 0; + std::set surface_1_boundary_attributes({8}); + std::set surface_2_boundary_attributes({9}); + solid_solver.addContactInteraction(contact_interaction_id, surface_1_boundary_attributes, surface_2_boundary_attributes, contact_options); + + //Finalize the data structures + solid_solver.completeSetup(); + + std::string visit_name = name + "_visit"; + solid_solver.outputStateToDisk(visit_name); + + solid_solver.completeSetup(); + + //Perform the quasi-static solve + double dt = 1.0; + + for (int i{0}; i < 100; ++i) { + solid_solver.advanceTimestep(dt); + visit_dc.SetCycle(i); + visit_dc.SetTime((i+1)*dt); + visit_dc.Save(); + + //Output the sidre-based plot files + solid_solver.outputStateToDisk(visit_name); + } + + return 0; + } + + + + + + + diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp index 884b4531d..c451b1b53 100644 --- a/src/serac/numerics/equation_solver.cpp +++ b/src/serac/numerics/equation_solver.cpp @@ -87,6 +87,7 @@ class NewtonSolver : public mfem::NewtonSolver { using real_t = mfem::real_t; real_t norm, norm_goal = 0; + norm = initial_norm = evaluateNorm(x, r); if (print_options.first_and_last && !print_options.iterations) { diff --git a/src/serac/physics/tests/contact_finite_diff.cpp b/src/serac/physics/tests/contact_finite_diff.cpp index 23a8bb913..095fbed4c 100644 --- a/src/serac/physics/tests/contact_finite_diff.cpp +++ b/src/serac/physics/tests/contact_finite_diff.cpp @@ -34,7 +34,7 @@ TEST_P(ContactFiniteDiff, patch) constexpr int p = 1; constexpr int dim = 2; - constexpr double eps = 1.0e-7; + constexpr double eps = 1.0e-4; MPI_Barrier(MPI_COMM_WORLD); @@ -48,7 +48,7 @@ TEST_P(ContactFiniteDiff, patch) double shift = eps * 10; // clang-format off auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.999}).translate({shift, 0.0}).bdrAttribInfo() + shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.9}).translate({shift, 0.0}).bdrAttribInfo() .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5) // shift up height of element diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp index 4a7d80b81..c5aab0ee8 100644 --- a/src/serac/physics/tests/contact_patch_2D.cpp +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -6,7 +6,6 @@ #include "serac/physics/solid_mechanics_contact.hpp" -#include #include #include #include @@ -24,9 +23,6 @@ #include "serac/physics/materials/solid_material.hpp" #include "serac/serac_config.hpp" #include "serac/infrastructure/application_manager.hpp" -#include - - namespace serac { @@ -48,9 +44,9 @@ TEST_P(ContactTest, patch) // Construct the appropriate dimension mesh and give it to the data store auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(1 , 1).translate({0.0, 1.0}).bdrAttribInfo() + shared::MeshBuilder::SquareMesh(64 , 64).translate({0.0, 1.0}).bdrAttribInfo() .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), - shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); + shared::MeshBuilder::SquareMesh(64, 64).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); mfem::VisItDataCollection visit_dc("contact_patch_visit", &mesh->mfemParMesh()); @@ -91,21 +87,21 @@ TEST_P(ContactTest, patch) ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement = std::get<0>(GetParam()), .type = ContactType::Frictionless, - .penalty = 0.2, - .penalty2 = 0.0, + .penalty = 25000, + .penalty2 = 25000, .jacobian = std::get<1>(GetParam())}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, name, mesh); - double K = 10.0; - double G = 0.25; + double K = 1000.0; + double G = 1200; solid_mechanics::NeoHookean mat{1.0, K, G}; solid_solver.setMaterial(mat, mesh->entireBody()); // Define the function for the initial displacement and boundary condition // constexpr int dim = 2; - auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.1}}; }; + auto applied_disp_function = [](tensor, auto) { return tensor{{0, -0.01}}; }; // Define a boundary attribute set and specify initial / boundary conditions solid_solver.setFixedBCs(mesh->domain("x0_faces"), Component::X); @@ -124,44 +120,38 @@ TEST_P(ContactTest, patch) // Perform the quasi-static solve double dt = 1.0; solid_solver.advanceTimestep(dt); - solid_solver.advanceTimestep(dt); - // solid_solver.advanceTimestep(dt); // solid_solver.advanceTimestep(dt); // Output the sidre-based plot files solid_solver.outputStateToDisk(paraview_name); // Check the l2 norm of the displacement dofs - // auto c = 1.0; auto c = (3.0 * K - 2.0 * G) / (3.0 * K + G); mfem::VectorFunctionCoefficient elasticity_sol_coeff(2, [c](const mfem::Vector& x, mfem::Vector& u) { - // u[0] = 0.0; - // u[1] = -0.01 * c * x[1]; - u[0] = c * -0.1 * x[0]; - u[1] = -0.1 * x[1]; + u[0] = 0.005 * c * x[0]; + u[1] = -0.005 * x[1]; // u[2] = -0.5 * 0.01 * x[2]; }); - mfem::ParFiniteElementSpace elasticity_fes(solid_solver.displacement().space()); + mfem::ParFiniteElementSpace elasticity_fes(solid_solver.reactions().space()); mfem::ParGridFunction elasticity_sol(&elasticity_fes); elasticity_sol.ProjectCoefficient(elasticity_sol_coeff); mfem::ParGridFunction approx_error(elasticity_sol); approx_error -= solid_solver.displacement().gridFunction(); auto approx_error_l2 = mfem::ParNormlp(approx_error, 2, MPI_COMM_WORLD); - EXPECT_NEAR(0.0, approx_error_l2, 1.0e-2); + EXPECT_NEAR(0.0, approx_error_l2, 1.0e-3); } INSTANTIATE_TEST_SUITE_P( tribol, ContactTest, testing::Values( - std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ") -)); + std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"))); } // namespace serac int main(int argc, char* argv[]) { -feenableexcept(FE_INVALID | FE_OVERFLOW); + testing::InitGoogleTest(&argc, argv); serac::ApplicationManager applicationManager(argc, argv); diff --git a/src/serac/physics/tests/tribol_finite_diff.cpp b/src/serac/physics/tests/tribol_finite_diff.cpp index c59237fd6..193102c87 100644 --- a/src/serac/physics/tests/tribol_finite_diff.cpp +++ b/src/serac/physics/tests/tribol_finite_diff.cpp @@ -42,10 +42,10 @@ TEST_P(TribolFiniteDiff, patch) // Construct the appropriate dimension mesh and give it to the data store - double shift = eps * 10; + // double shift = eps * 10; // clang-format off auto pmesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.999}).translate({shift, 0.0}).bdrAttribInfo() + shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.999}).translate({0.0, 0.0}).bdrAttribInfo() .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), shared::MeshBuilder::SquareMesh(1, 1).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5) // shift up height of element From 29656e79b40c85deda3bad99d647807596e01dca Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Thu, 14 Aug 2025 09:12:06 -0700 Subject: [PATCH 07/11] patch test --- tribol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tribol b/tribol index 2ac69027c..1379ab29e 160000 --- a/tribol +++ b/tribol @@ -1 +1 @@ -Subproject commit 2ac69027ce48b47c1f30432cf60907b929f485f4 +Subproject commit 1379ab29e9926e4b367bd7a016020efe0c6f2a6b From 8d2daa8a7a4cc4bba85c3ea9c74a954da1fbcce4 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Thu, 14 Aug 2025 17:38:06 -0700 Subject: [PATCH 08/11] smooth mortar fixes --- examples/contact/ironing_2D.cpp | 19 ++++++++++------- .../physics/tests/contact_finite_diff.cpp | 4 ++-- src/serac/physics/tests/contact_patch_2D.cpp | 21 +++++++++++-------- 3 files changed, 25 insertions(+), 19 deletions(-) diff --git a/examples/contact/ironing_2D.cpp b/examples/contact/ironing_2D.cpp index 73541c2ac..a50ae2211 100644 --- a/examples/contact/ironing_2D.cpp +++ b/examples/contact/ironing_2D.cpp @@ -64,8 +64,8 @@ int main(int argc, char* argv[]) // std::shared_ptr mesh = std::make_shared(filename, "ironing_2D_mesh", 2, 0); auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(32, 32).updateBdrAttrib(1, 6).updateBdrAttrib(3, 9).bdrAttribInfo().scale({1.0, 1.0}), - shared::MeshBuilder::SquareMesh(32, 32).scale({0.25, 0.25}).translate({0.0, 1.0}).updateBdrAttrib(3, 5).updateBdrAttrib(1,8).updateAttrib(1, 2)}), "iroing_2D_mesh", 0, 0); + shared::MeshBuilder::SquareMesh(32, 32).updateBdrAttrib(1, 6).updateBdrAttrib(3, 9).bdrAttribInfo().scale({1.0, 0.5}), + shared::MeshBuilder::SquareMesh(8, 8).scale({0.25, 0.25}).translate({0.0, 0.5}).updateBdrAttrib(3, 5).updateBdrAttrib(1,8).updateAttrib(1, 2)}), "iroing_2D_mesh", 0, 0); serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level=0}; @@ -79,7 +79,7 @@ int main(int argc, char* argv[]) return 1; #endif - serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton, + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::TrustRegion, .relative_tol = 1.0e-8, .absolute_tol = 1.0e-10, .max_iterations = 50, @@ -106,7 +106,7 @@ int main(int argc, char* argv[]) serac::FiniteElementState G_field(serac::StateManager::newState(serac::L2<0>{}, "shear_mod", mesh->tag())); - mfem::Vector G_values({0.25, 2.5}); + mfem::Vector G_values({0.1, 2.5}); mfem::PWConstCoefficient G_coeff(G_values); G_field.project(G_coeff); solid_solver.setParameter(1, G_field); @@ -124,12 +124,12 @@ int main(int argc, char* argv[]) serac::tensor u{}; // std::cout << "T ========= " << t << std::endl; if (t <= init_steps + 1.0e-12) { - u[1] = -t * 0.2 / init_steps; + u[1] = -t * 0.05 / init_steps; // std::cout << "In IF statement. u[1] = " << u[1] << " and t = " << t << std::endl; } else { - u[0] = (t - init_steps) * 0.0005; - u[1] = -0.2; + u[0] = (t - init_steps) * 0.01; + u[1] = -0.05; // std::cout << "in ELSE statement. u[1] = " << u[1] << " and u[0] = " << u[0] << " and t = " << t << std::endl; } return u; @@ -156,11 +156,14 @@ int main(int argc, char* argv[]) //Perform the quasi-static solve double dt = 1.0; - for (int i{0}; i < 100; ++i) { + for (int i{0}; i < 150; ++i) { solid_solver.advanceTimestep(dt); visit_dc.SetCycle(i); visit_dc.SetTime((i+1)*dt); visit_dc.Save(); + if (i == 88) { + std::cout << i << std::endl; + } //Output the sidre-based plot files solid_solver.outputStateToDisk(visit_name); diff --git a/src/serac/physics/tests/contact_finite_diff.cpp b/src/serac/physics/tests/contact_finite_diff.cpp index 095fbed4c..d8eef7557 100644 --- a/src/serac/physics/tests/contact_finite_diff.cpp +++ b/src/serac/physics/tests/contact_finite_diff.cpp @@ -34,7 +34,7 @@ TEST_P(ContactFiniteDiff, patch) constexpr int p = 1; constexpr int dim = 2; - constexpr double eps = 1.0e-4; + constexpr double eps = 0.7; MPI_Barrier(MPI_COMM_WORLD); @@ -45,7 +45,7 @@ TEST_P(ContactFiniteDiff, patch) // Construct the appropriate dimension mesh and give it to the data store - double shift = eps * 10; + double shift = 0.0; // clang-format off auto mesh = std::make_shared(shared::MeshBuilder::Unify({ shared::MeshBuilder::SquareMesh(1, 1).translate({0.0, 0.9}).translate({shift, 0.0}).bdrAttribInfo() diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp index c5aab0ee8..11baa22b3 100644 --- a/src/serac/physics/tests/contact_patch_2D.cpp +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include "axom/slic/core/SimpleLogger.hpp" #include @@ -85,17 +86,17 @@ TEST_P(ContactTest, patch) .print_level = 1}; ContactOptions contact_options{.method = ContactMethod::SmoothMortar, - .enforcement = std::get<0>(GetParam()), + .enforcement =serac::ContactEnforcement::Penalty, .type = ContactType::Frictionless, - .penalty = 25000, - .penalty2 = 25000, - .jacobian = std::get<1>(GetParam())}; + .penalty = 100000, + .penalty2 = 0, + .jacobian = serac::ContactJacobian::Exact}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, name, mesh); - double K = 1000.0; - double G = 1200; + double K = 10000.0; + double G = 10; solid_mechanics::NeoHookean mat{1.0, K, G}; solid_solver.setMaterial(mat, mesh->entireBody()); @@ -126,7 +127,8 @@ TEST_P(ContactTest, patch) solid_solver.outputStateToDisk(paraview_name); // Check the l2 norm of the displacement dofs - auto c = (3.0 * K - 2.0 * G) / (3.0 * K + G); + auto c = (3.0 * K - 2.0 * G) / ((3.0 * K + 2*G)); + // auto c = 0.0; mfem::VectorFunctionCoefficient elasticity_sol_coeff(2, [c](const mfem::Vector& x, mfem::Vector& u) { u[0] = 0.005 * c * x[0]; u[1] = -0.005 * x[1]; @@ -138,13 +140,14 @@ TEST_P(ContactTest, patch) mfem::ParGridFunction approx_error(elasticity_sol); approx_error -= solid_solver.displacement().gridFunction(); auto approx_error_l2 = mfem::ParNormlp(approx_error, 2, MPI_COMM_WORLD); - EXPECT_NEAR(0.0, approx_error_l2, 1.0e-3); + EXPECT_NEAR(0.0, approx_error_l2, 1.0e-2); } INSTANTIATE_TEST_SUITE_P( tribol, ContactTest, testing::Values( - std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"))); + std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Approximate, "penalty_approxJ"))); + // std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"))); } // namespace serac From f4c10834f155191644cf50103d5c672df03ddd2c Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Fri, 15 Aug 2025 14:39:27 -0700 Subject: [PATCH 09/11] Added twisted 2D contact test --- examples/contact/CMakeLists.txt | 2 + examples/contact/ironing_2D.cpp | 7 +- examples/contact/twisted_ironing_2D.cpp | 183 ++++++++++++++++++++++++ 3 files changed, 187 insertions(+), 5 deletions(-) create mode 100644 examples/contact/twisted_ironing_2D.cpp diff --git a/examples/contact/CMakeLists.txt b/examples/contact/CMakeLists.txt index 79acabe48..391f944d6 100644 --- a/examples/contact/CMakeLists.txt +++ b/examples/contact/CMakeLists.txt @@ -10,6 +10,8 @@ if(TRIBOL_FOUND AND STRUMPACK_DIR) ironing.cpp sphere.cpp twist.cpp + ironing_2D.cpp + twisted_ironing_2D.cpp ) foreach(filename ${CONTACT_EXAMPLES_SOURCES}) diff --git a/examples/contact/ironing_2D.cpp b/examples/contact/ironing_2D.cpp index a50ae2211..58e3c2232 100644 --- a/examples/contact/ironing_2D.cpp +++ b/examples/contact/ironing_2D.cpp @@ -146,7 +146,7 @@ int main(int argc, char* argv[]) solid_solver.addContactInteraction(contact_interaction_id, surface_1_boundary_attributes, surface_2_boundary_attributes, contact_options); //Finalize the data structures - solid_solver.completeSetup(); + // solid_solver.completeSetup(); std::string visit_name = name + "_visit"; solid_solver.outputStateToDisk(visit_name); @@ -156,14 +156,11 @@ int main(int argc, char* argv[]) //Perform the quasi-static solve double dt = 1.0; - for (int i{0}; i < 150; ++i) { + for (int i{0}; i < 1; ++i) { solid_solver.advanceTimestep(dt); visit_dc.SetCycle(i); visit_dc.SetTime((i+1)*dt); visit_dc.Save(); - if (i == 88) { - std::cout << i << std::endl; - } //Output the sidre-based plot files solid_solver.outputStateToDisk(visit_name); diff --git a/examples/contact/twisted_ironing_2D.cpp b/examples/contact/twisted_ironing_2D.cpp new file mode 100644 index 000000000..ce19e16cc --- /dev/null +++ b/examples/contact/twisted_ironing_2D.cpp @@ -0,0 +1,183 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include +#include + +#include "axom/slic.hpp" + +#include "mfem.hpp" + +#include "serac/numerics/solver_config.hpp" +#include "serac/physics/contact/contact_config.hpp" +#include "serac/physics/materials/parameterized_solid_material.hpp" +#include "serac/physics/solid_mechanics.hpp" +#include "serac/physics/solid_mechanics_contact.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/serac.hpp" +#include + +#include "serac/physics/contact/contact_config.hpp" +#include "serac/physics/solid_mechanics_contact.hpp" + +#include +#include +#include +#include +#include +#include +#include "axom/slic/core/SimpleLogger.hpp" +#include "mfem.hpp" + +#include "shared/mesh/MeshBuilder.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" +#include "serac/physics/boundary_conditions/components.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/mesh.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/serac_config.hpp" +#include "serac/infrastructure/application_manager.hpp" +#include + + + +int main(int argc, char* argv[]) + { + + //Initialize and automatically finalize MPI and other libraries + serac::ApplicationManager applicationManager(argc, argv); + + // NOTE: p must be equal to 1 to work with Tribol's mortar method + constexpr int p = 1; + + //NOTE: dim must be equal to 2 + constexpr int dim = 2; + + //Create DataStore + std::string name = "twisted_contact_ironing_2D_example"; + axom::sidre::DataStore datastore; + serac::StateManager::initialize(datastore, name + "_data"); + + + auto mesh = std::make_shared(shared::MeshBuilder::Unify({ + shared::MeshBuilder::SquareMesh(32, 32).updateBdrAttrib(1, 6).updateBdrAttrib(3, 9).scale({1.0, 0.5}), + shared::MeshBuilder::SquareMesh(8, 8).scale({0.25, 0.25}).translate({0.0, 0.5}).updateBdrAttrib(3, 5).updateBdrAttrib(1, 8).updateBdrAttrib(2, 8).updateBdrAttrib(4, 8).updateAttrib(1, 2)}), "ironing_2D_mesh", 0, 0); + + serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level=0}; + + mfem::VisItDataCollection visit_dc("contact_ironing_twist_vist", &mesh->mfemParMesh()); + + visit_dc.SetPrefixPath("visit_out"); + visit_dc.Save(); + + #ifndef MFEM_USE_STRUMPACK + SLIC_INFO_ROOT("Contact requires MFEM built with strumpack."); + return 1; + #endif + + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::TrustRegion, + .relative_tol = 1.0e-8, + .absolute_tol = 1.0e-10, + .max_iterations = 50, + .print_level = 1}; + + serac::ContactOptions contact_options{.method = serac::ContactMethod::SmoothMortar, + .enforcement = serac::ContactEnforcement::Penalty, + .penalty = 750, + .penalty2 = 0.0, + .jacobian = serac::ContactJacobian::Exact}; + + serac::SolidMechanicsContact, serac::L2<0>>> solid_solver( + nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, mesh, {"bulk_mod", "shear_mod"}); + + + serac::FiniteElementState K_field(serac::StateManager::newState(serac::L2<0>{}, "bulk_mod", mesh->tag())); + + mfem::Vector K_values({10.0, 100.0}); + mfem::PWConstCoefficient K_coeff(K_values); + K_field.project(K_coeff); + solid_solver.setParameter(0, K_field); + + serac::FiniteElementState G_field(serac::StateManager::newState(serac::L2<0>{}, "shear_mod", mesh->tag())); + + mfem::Vector G_values({0.25, 2.5}); + mfem::PWConstCoefficient G_coeff(G_values); + G_field.project(G_coeff); + solid_solver.setParameter(1, G_field); + + serac::solid_mechanics::ParameterizedNeoHookeanSolid mat{1.0, 0.0, 0.0}; + solid_solver.setMaterial(serac::DependsOn<0, 1>{}, mat, mesh->entireBody()); + + //Pass the BC information to the solver object + mesh->addDomainOfBoundaryElements("bottom_of_subtrate", serac::by_attr(6)); + solid_solver.setFixedBCs((mesh->domain("bottom_of_subtrate"))); + + mesh->addDomainOfBoundaryElements("top of indenter", serac::by_attr(5)); + const serac::tensor r0{{0.125, 0.625}}; + auto applied_displacement = [r0](serac::tensor x, double t) { + constexpr double init_steps = 10.0; + constexpr double theta_max = 80.0 * M_PI / 180.0; + serac::tensor u{}; + if (t <= init_steps + 1.0e-12) { + u[1] = -t * 0.03 / init_steps; + } + else { + double hm = (t - init_steps) * 0.01; //horizontal movement + double theta = theta_max * hm; //current rotation angle + double cos_theta = std::cos(theta); + double sin_theta = std::sin(theta); + + //Rotate about r0 + serac::tensor y {{x[0] - r0[0], x[1] - r0[1]}}; + serac::tensor y_rot {{cos_theta*y[0] - sin_theta*y[1], sin_theta*y[0] + cos_theta*y[1]}}; + + u[0] = (y_rot[0] - y[0]) + 0.01 * (t - init_steps); + u[1] = (y_rot[1] - y[1]) - 0.03; + } + return u; + }; + + + solid_solver.setDisplacementBCs(applied_displacement, mesh->domain("top of indenter")); + + //Add contact interaction + + auto contact_interaction_id = 0; + std::set surface_1_boundary_attributes({8}); + std::set surface_2_boundary_attributes({9}); + solid_solver.addContactInteraction(contact_interaction_id, surface_1_boundary_attributes, surface_2_boundary_attributes, contact_options); + + + std::string visit_name = name + "_visit"; + solid_solver.outputStateToDisk(visit_name); + + solid_solver.completeSetup(); + + //Perform the quasi-static solve + double dt = 1.0; + + for (int i{0}; i < 110; ++i) { + solid_solver.advanceTimestep(dt); + visit_dc.SetCycle(i); + visit_dc.SetTime((i+1)*dt); + visit_dc.Save(); + //Output the sidre-based plot files + solid_solver.outputStateToDisk(visit_name); + } + + return 0; + } + + + + + + + + + + + From 1b7463151e222ce4b6cc8465249ea444c647f517 Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Fri, 15 Aug 2025 14:48:50 -0700 Subject: [PATCH 10/11] update tribol --- tribol | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tribol b/tribol index 1379ab29e..5c56ea26e 160000 --- a/tribol +++ b/tribol @@ -1 +1 @@ -Subproject commit 1379ab29e9926e4b367bd7a016020efe0c6f2a6b +Subproject commit 5c56ea26e84bda0595bb935d5ddc292747759a3d From 213230185b770bc0e423d9010f677ba4a209ffdf Mon Sep 17 00:00:00 2001 From: Ryan Lutz Date: Mon, 8 Sep 2025 17:30:48 -0700 Subject: [PATCH 11/11] x and y component wise error --- examples/contact/ironing_2D.cpp | 10 +-- src/serac/physics/tests/contact_patch_2D.cpp | 94 +++++++++++++++++--- 2 files changed, 87 insertions(+), 17 deletions(-) diff --git a/examples/contact/ironing_2D.cpp b/examples/contact/ironing_2D.cpp index 58e3c2232..8b8d113de 100644 --- a/examples/contact/ironing_2D.cpp +++ b/examples/contact/ironing_2D.cpp @@ -64,8 +64,8 @@ int main(int argc, char* argv[]) // std::shared_ptr mesh = std::make_shared(filename, "ironing_2D_mesh", 2, 0); auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(32, 32).updateBdrAttrib(1, 6).updateBdrAttrib(3, 9).bdrAttribInfo().scale({1.0, 0.5}), - shared::MeshBuilder::SquareMesh(8, 8).scale({0.25, 0.25}).translate({0.0, 0.5}).updateBdrAttrib(3, 5).updateBdrAttrib(1,8).updateAttrib(1, 2)}), "iroing_2D_mesh", 0, 0); + shared::MeshBuilder::SquareMesh(4, 4).updateBdrAttrib(1, 6).updateBdrAttrib(3, 9).bdrAttribInfo().scale({1.0, 0.5}), + shared::MeshBuilder::SquareMesh(1, 1).scale({0.25, 0.25}).translate({0.0, 0.5}).updateBdrAttrib(3, 5).updateBdrAttrib(1,8).updateAttrib(1, 2)}), "iroing_2D_mesh", 0, 0); serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level=0}; @@ -124,12 +124,12 @@ int main(int argc, char* argv[]) serac::tensor u{}; // std::cout << "T ========= " << t << std::endl; if (t <= init_steps + 1.0e-12) { - u[1] = -t * 0.05 / init_steps; + u[1] = -t * 0.15 / init_steps; // std::cout << "In IF statement. u[1] = " << u[1] << " and t = " << t << std::endl; } else { u[0] = (t - init_steps) * 0.01; - u[1] = -0.05; + u[1] = -0.15; // std::cout << "in ELSE statement. u[1] = " << u[1] << " and u[0] = " << u[0] << " and t = " << t << std::endl; } return u; @@ -156,7 +156,7 @@ int main(int argc, char* argv[]) //Perform the quasi-static solve double dt = 1.0; - for (int i{0}; i < 1; ++i) { + for (int i{0}; i < 200; ++i) { solid_solver.advanceTimestep(dt); visit_dc.SetCycle(i); visit_dc.SetTime((i+1)*dt); diff --git a/src/serac/physics/tests/contact_patch_2D.cpp b/src/serac/physics/tests/contact_patch_2D.cpp index 11baa22b3..833eb24ee 100644 --- a/src/serac/physics/tests/contact_patch_2D.cpp +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -9,11 +9,17 @@ #include #include #include +#include +#include +#include +#include +#include #include #include #include #include "axom/slic/core/SimpleLogger.hpp" #include +#include #include "mfem.hpp" #include "shared/mesh/MeshBuilder.hpp" @@ -24,6 +30,15 @@ #include "serac/physics/materials/solid_material.hpp" #include "serac/serac_config.hpp" #include "serac/infrastructure/application_manager.hpp" +#include + + +// static void enable_fpe() { +// // trap on invalid ops (NaN), divide-by-zero, and overflow +// feclearexcept(FE_ALL_EXCEPT); +// feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW); + +// } namespace serac { @@ -45,16 +60,16 @@ TEST_P(ContactTest, patch) // Construct the appropriate dimension mesh and give it to the data store auto mesh = std::make_shared(shared::MeshBuilder::Unify({ - shared::MeshBuilder::SquareMesh(64 , 64).translate({0.0, 1.0}).bdrAttribInfo() + shared::MeshBuilder::SquareMesh(100,100 ).translate({0.0, 1.0}).bdrAttribInfo() .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), - shared::MeshBuilder::SquareMesh(64, 64).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); + shared::MeshBuilder::SquareMesh(80, 80).bdrAttribInfo().updateBdrAttrib(4, 7).updateBdrAttrib(1, 8).updateBdrAttrib(3, 5)}), "patch_mesh_2D", 0, 0); mfem::VisItDataCollection visit_dc("contact_patch_visit", &mesh->mfemParMesh()); visit_dc.SetPrefixPath("visit_out"); visit_dc.Save(); - + mesh->addDomainOfBoundaryElements("x0_faces", serac::by_attr(7)); @@ -79,23 +94,24 @@ TEST_P(ContactTest, patch) return; #endif - NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, + NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::NewtonLineSearch, .relative_tol = 1.0e-13, .absolute_tol = 1.0e-13, .max_iterations = 20, + .max_line_search_iterations = 12, .print_level = 1}; ContactOptions contact_options{.method = ContactMethod::SmoothMortar, .enforcement =serac::ContactEnforcement::Penalty, .type = ContactType::Frictionless, - .penalty = 100000, + .penalty = 10000000, .penalty2 = 0, .jacobian = serac::ContactJacobian::Exact}; SolidMechanicsContact solid_solver(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, name, mesh); - double K = 10000.0; + double K = 1000.0; double G = 10; solid_mechanics::NeoHookean mat{1.0, K, G}; solid_solver.setMaterial(mat, mesh->entireBody()); @@ -115,6 +131,7 @@ TEST_P(ContactTest, patch) // Finalize the data structures solid_solver.completeSetup(); + std::string paraview_name = name + "_paraview"; solid_solver.outputStateToDisk(paraview_name); @@ -137,16 +154,69 @@ TEST_P(ContactTest, patch) mfem::ParFiniteElementSpace elasticity_fes(solid_solver.reactions().space()); mfem::ParGridFunction elasticity_sol(&elasticity_fes); elasticity_sol.ProjectCoefficient(elasticity_sol_coeff); - mfem::ParGridFunction approx_error(elasticity_sol); - approx_error -= solid_solver.displacement().gridFunction(); - auto approx_error_l2 = mfem::ParNormlp(approx_error, 2, MPI_COMM_WORLD); - EXPECT_NEAR(0.0, approx_error_l2, 1.0e-2); + // mfem::ParGridFunction approx_error(elasticity_sol); + // approx_error -= solid_solver.displacement().gridFunction(); + // auto approx_error_l2 = mfem::ParNormlp(approx_error, 2, MPI_COMM_WORLD); + // EXPECT_NEAR(0.0, approx_error_l2, 1.0e-2); + + //Set up test to only look at y component of error********* +const mfem::ParFiniteElementSpace& u_space_const = solid_solver.displacement().space(); +auto& u_space = const_cast(u_space_const); +mfem::ParGridFunction U_exact(&u_space); +U_exact.ProjectCoefficient(elasticity_sol_coeff); + +// Numerical displacement +const mfem::ParGridFunction& U_num = solid_solver.displacement().gridFunction(); + +// Overall Error +mfem::ParGridFunction U_err(U_exact); +U_err -= U_num; +const double L2_err_vec = mfem::ParNormlp(U_err, 2, MPI_COMM_WORLD); +std::cout << "L2_err_vec = " << L2_err_vec << std::endl; + +//y-component error +const mfem::FiniteElementCollection* fec = u_space.FEColl(); +mfem::ParFiniteElementSpace y_fes(&mesh->mfemParMesh(), fec, /*vdim=*/1, u_space.GetOrdering()); //builds scalar space on same mesh + +mfem::ParGridFunction uy_ex(&y_fes), uy_num(&y_fes); +const int n = y_fes.GetNDofs(); + +for (int i = 0; i < n; ++i) { + uy_ex(i) = U_exact(n*1 + i); + uy_num(i) = U_num (n*1 + i); +} + +//Same thing for x forces. +mfem::ParGridFunction ux_ex(&y_fes), ux_num(&y_fes); + +for( int i = 0; i < n; ++i) { + ux_ex(i) = U_exact(i); + ux_num(i) = U_num(i); +} + +mfem::ParGridFunction uy_err(uy_ex); +mfem::ParGridFunction ux_err(ux_ex); +uy_err -= uy_num; +ux_err -= ux_num; +const double L2_err_y = mfem::ParNormlp(uy_err, 2, MPI_COMM_WORLD); +const double L2_err_x = mfem::ParNormlp(ux_err, 2, MPI_COMM_WORLD); +std::cout << "L2_err_y = " << L2_err_y << std::endl; +std::cout << "L2_err_x = " << L2_err_x << std::endl; + +EXPECT_NEAR(0.0, L2_err_vec, 1e-2); +EXPECT_NEAR(0.0, L2_err_y, 1e-2); +EXPECT_NEAR(0.0, L2_err_x, 1e-2); + + +std::cout << "check = " + << std::abs(L2_err_vec*L2_err_vec - (L2_err_x*L2_err_x + L2_err_y*L2_err_y)) + << "\n"; } INSTANTIATE_TEST_SUITE_P( tribol, ContactTest, testing::Values( - std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Approximate, "penalty_approxJ"))); + std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_approxJ"))); // std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"))); } // namespace serac @@ -155,7 +225,7 @@ int main(int argc, char* argv[]) { - + // enable_fpe(); testing::InitGoogleTest(&argc, argv); serac::ApplicationManager applicationManager(argc, argv); return RUN_ALL_TESTS();