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 new file mode 100644 index 000000000..8b8d113de --- /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(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}; + + 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::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, + .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.1, 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.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.15; + // 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 < 200; ++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/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; + } + + + + + + + + + + + diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp index 0af69dcff..c451b1b53 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 @@ -86,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) { @@ -192,6 +194,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/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_finite_diff.cpp b/src/serac/physics/tests/contact_finite_diff.cpp index 91fce7e4c..d8eef7557 100644 --- a/src/serac/physics/tests/contact_finite_diff.cpp +++ b/src/serac/physics/tests/contact_finite_diff.cpp @@ -32,9 +32,9 @@ 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; + constexpr double eps = 0.7; MPI_Barrier(MPI_COMM_WORLD); @@ -45,26 +45,26 @@ 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::CubeMesh(1, 1, 1), - shared::MeshBuilder::CubeMesh(1, 1, 1) + 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 - .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 new file mode 100644 index 000000000..833eb24ee --- /dev/null +++ b/src/serac/physics/tests/contact_patch_2D.cpp @@ -0,0 +1,232 @@ +// 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 +#include +#include +#include +#include +#include +#include "axom/slic/core/SimpleLogger.hpp" +#include +#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" +#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 { + +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(100,100 ).translate({0.0, 1.0}).bdrAttribInfo() + .updateBdrAttrib(4, 7).updateBdrAttrib(3, 9).updateBdrAttrib(1, 6), + 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)); + 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::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 = 10000000, + .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 = 10; + 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.01}}; }; + + // 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 + 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]; + // 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-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::Exact, "penalty_approxJ"))); + // std::make_tuple(ContactEnforcement::Penalty, ContactJacobian::Exact, "penalty_exactJ"))); + +} // namespace serac + +int main(int argc, char* argv[]) +{ + + + // enable_fpe(); + testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/serac/physics/tests/tribol_finite_diff.cpp b/src/serac/physics/tests/tribol_finite_diff.cpp index 2d9b41fec..193102c87 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; + + // 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({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 + + // // 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 diff --git a/tribol b/tribol index 8351aa4e1..5c56ea26e 160000 --- a/tribol +++ b/tribol @@ -1 +1 @@ -Subproject commit 8351aa4e104fce00ee25c760a5c3f3c882557d50 +Subproject commit 5c56ea26e84bda0595bb935d5ddc292747759a3d