diff --git a/src/serac/physics/state/state_manager.cpp b/src/serac/physics/state/state_manager.cpp index 41a2443cd..538e9f2e4 100644 --- a/src/serac/physics/state/state_manager.cpp +++ b/src/serac/physics/state/state_manager.cpp @@ -53,13 +53,28 @@ double StateManager::newDataCollection(const std::string& name, const std::optio // Functional needs the nodal grid function and neighbor data in the mesh + // Determine if the existing nodal grid function is discontinuous. This + // indicates that the mesh is periodic and the new nodal grid function must also + // be discontinuous. + bool is_discontinuous = false; + auto nodes = mesh(name).GetNodes(); + if (nodes) { + is_discontinuous = nodes->FESpace()->FEColl()->GetContType() == mfem::FiniteElementCollection::DISCONTINUOUS; + SLIC_WARNING_ROOT( + "Periodic mesh detected! This will only work on translational periodic surfaces for vector H1 fields and " + "has not been thoroughly tested. Proceed at your own risk."); + } + // This mfem call ensures the mesh contains an H1 grid function describing nodal // cordinates. The parameters do the following: // 1. Sets the order of the mesh to p = 1 - // 2. Uses a continuous (i.e. H1) finite element space + // 2. Uses the existing continuity of the mesh finite element space (periodic meshes are discontinuous) // 3. Uses the spatial dimension as the mesh dimension (i.e. it is not a lower dimension manifold) // 4. Uses nodal instead of VDIM ordering (i.e. xxxyyyzzz instead of xyzxyzxyz) - mesh(name).SetCurvature(1, false, -1, mfem::Ordering::byNODES); + mesh(name).SetCurvature(1, is_discontinuous, -1, mfem::Ordering::byNODES); + + // Sidre will destruct the nodal grid function instead of the mesh + mesh(name).SetNodesOwner(false); // Generate the face neighbor information in the mesh. This is needed by the face restriction // operators used by Functional @@ -221,6 +236,29 @@ void StateManager::save(const double t, const int cycle, const std::string& mesh mfem::ParMesh* StateManager::setMesh(std::unique_ptr pmesh, const std::string& mesh_tag) { + // Determine if the existing nodal grid function is discontinuous. This + // indicates that the mesh is periodic and the new nodal grid function must also + // be discontinuous. + bool is_discontinuous = false; + auto nodes = pmesh->GetNodes(); + if (nodes) { + is_discontinuous = nodes->FESpace()->FEColl()->GetContType() == mfem::FiniteElementCollection::DISCONTINUOUS; + SLIC_WARNING_ROOT( + "Periodic mesh detected! This will only work on translational periodic surfaces for vector H1 fields and " + "has not been thoroughly tested. Proceed at your own risk."); + } + + // This mfem call ensures the mesh contains an H1 grid function describing nodal + // cordinates. The parameters do the following: + // 1. Sets the order of the mesh to p = 1 + // 2. Uses the existing continuity of the mesh finite element space (periodic meshes are discontinuous) + // 3. Uses the spatial dimension as the mesh dimension (i.e. it is not a lower dimension manifold) + // 4. Uses nodal instead of VDIM ordering (i.e. xxxyyyzzz instead of xyzxyzxyz) + pmesh->SetCurvature(1, is_discontinuous, -1, mfem::Ordering::byNODES); + + // Sidre will destruct the nodal grid function instead of the mesh + pmesh->SetNodesOwner(false); + newDataCollection(mesh_tag); auto& datacoll = datacolls_.at(mesh_tag); datacoll.SetMesh(pmesh.release()); @@ -229,14 +267,6 @@ mfem::ParMesh* StateManager::setMesh(std::unique_ptr pmesh, const // Functional needs the nodal grid function and neighbor data in the mesh auto& new_pmesh = mesh(mesh_tag); - // This mfem call ensures the mesh contains an H1 grid function describing nodal - // cordinates. The parameters do the following: - // 1. Sets the order of the mesh to p = 1 - // 2. Uses a continuous (i.e. H1) finite element space - // 3. Uses the spatial dimension as the mesh dimension (i.e. it is not a lower dimension manifold) - // 4. Uses nodal instead of VDIM ordering (i.e. xxxyyyzzz instead of xyzxyzxyz) - new_pmesh.SetCurvature(1, false, -1, mfem::Ordering::byNODES); - // Generate the face neighbor information in the mesh. This is needed by the face restriction // operators used by Functional new_pmesh.ExchangeFaceNbrData(); diff --git a/src/serac/physics/tests/CMakeLists.txt b/src/serac/physics/tests/CMakeLists.txt index 38ae5f7f9..f551ade7a 100644 --- a/src/serac/physics/tests/CMakeLists.txt +++ b/src/serac/physics/tests/CMakeLists.txt @@ -36,6 +36,7 @@ set(solver_tests solid_legacy_reuse_mesh.cpp solid_legacy_adjoint.cpp solid.cpp + solid_periodic.cpp solid_shape.cpp thermal_legacy.cpp thermal_shape.cpp diff --git a/src/serac/physics/tests/solid_periodic.cpp b/src/serac/physics/tests/solid_periodic.cpp new file mode 100644 index 000000000..bd58bb3f3 --- /dev/null +++ b/src/serac/physics/tests/solid_periodic.cpp @@ -0,0 +1,129 @@ +// Copyright (c) 2019-2023, 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 "axom/slic/core/SimpleLogger.hpp" +#include +#include "mfem.hpp" + +#include "serac/serac_config.hpp" +#include "serac/mesh/mesh_utils.hpp" +#include "serac/physics/solid_mechanics.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/physics/materials/parameterized_solid_material.hpp" +#include "serac/physics/state/state_manager.hpp" + +namespace serac { + +using solid_mechanics::default_static_options; +using solid_mechanics::direct_static_options; + +TEST(SolidMechanics, Periodic) +{ + MPI_Barrier(MPI_COMM_WORLD); + + int serial_refinement = 0; + int parallel_refinement = 0; + + // Create DataStore + axom::sidre::DataStore datastore; + serac::StateManager::initialize(datastore, "solid_periodic"); + + // Construct the appropriate dimension mesh and give it to the data store + int nElem = 2; + double lx = 3.0e-1, ly = 3.0e-1, lz = 0.25e-1; + auto initial_mesh = + mfem::Mesh(mfem::Mesh::MakeCartesian3D(4 * nElem, 4 * nElem, nElem, mfem::Element::HEXAHEDRON, lx, ly, lz)); + + // Create translation vectors defining the periodicity + mfem::Vector x_translation({lx, 0.0, 0.0}); + std::vector translations = {x_translation}; + double tol = 1e-6; + + std::vector periodicMap = initial_mesh.CreatePeriodicVertexMapping(translations, tol); + + // Create the periodic mesh using the vertex mapping defined by the translation vectors + auto periodic_mesh = mfem::Mesh::MakePeriodic(initial_mesh, periodicMap); + auto mesh = mesh::refineAndDistribute(std::move(periodic_mesh), serial_refinement, parallel_refinement); + + serac::StateManager::setMesh(std::move(mesh)); + + constexpr int p = 1; + constexpr int dim = 3; + + // Construct and initialized the user-defined moduli to be used as a differentiable parameter in + // the solid physics module. + FiniteElementState user_defined_shear_modulus(StateManager::newState( + FiniteElementState::Options{.order = 1, .element_type = ElementType::L2, .name = "parameterized_shear"})); + + double shear_modulus_value = 1.0; + + user_defined_shear_modulus = shear_modulus_value; + + FiniteElementState user_defined_bulk_modulus(StateManager::newState( + FiniteElementState::Options{.order = 1, .element_type = ElementType::L2, .name = "parameterized_bulk"})); + + double bulk_modulus_value = 1.0; + + user_defined_bulk_modulus = bulk_modulus_value; + + // Construct a functional-based solid solver + SolidMechanics, L2

>> solid_solver(default_static_options, GeometricNonlinearities::On, + "solid_periodic"); + solid_solver.setParameter(0, user_defined_bulk_modulus); + solid_solver.setParameter(1, user_defined_shear_modulus); + + solid_mechanics::ParameterizedNeoHookeanSolid mat{1.0, 0.0, 0.0}; + solid_solver.setMaterial(DependsOn<0, 1>{}, mat); + + // Boundary conditions: + // Prescribe zero displacement at the supported end of the beam + std::set support = {2}; + auto zero_displacement = [](const mfem::Vector&, mfem::Vector& u) -> void { u = 0.0; }; + solid_solver.setDisplacementBCs(support, zero_displacement); + + double iniDispVal = 5.0e-6; + auto ini_displacement = [iniDispVal](const mfem::Vector&, mfem::Vector& u) -> void { u = iniDispVal; }; + solid_solver.setDisplacement(ini_displacement); + + tensor constant_force; + + constant_force[0] = 0.0; + constant_force[1] = 1.0e-2; + + if (dim == 3) { + constant_force[2] = 0.0; + } + + solid_mechanics::ConstantBodyForce force{constant_force}; + solid_solver.addBodyForce(force); + + // Finalize the data structures + solid_solver.completeSetup(); + + // Perform the quasi-static solve + double dt = 1.0; + solid_solver.advanceTimestep(dt); + + // Output the sidre-based plot files + solid_solver.outputState(); +} + +} // namespace serac + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + MPI_Init(&argc, &argv); + + axom::slic::SimpleLogger logger; + + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + + return result; +}