From 6a55990e09c4acaa944f529055afd02dd13d9f25 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 11 Oct 2024 15:33:49 -0700 Subject: [PATCH 01/17] Fix a missing const in Functional method --- src/serac/numerics/functional/functional.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index 402cab9b3..79311b68c 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -301,7 +301,7 @@ class Functional { /// @overload template - void AddDomainIntegral(Dimension, DependsOn, const lambda& integrand, Domain& domain, + void AddDomainIntegral(Dimension, DependsOn, const lambda& integrand, const Domain& domain, std::shared_ptr> qdata = NoQData) { if (domain.mesh_.GetNE() == 0) return; From e5b101f1e4745c7f46b55abd8042f5adb96de721 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 11 Oct 2024 15:43:50 -0700 Subject: [PATCH 02/17] Implement multi-material solid mechanics, test by eye --- src/serac/physics/solid_mechanics.hpp | 8 +- src/serac/physics/tests/CMakeLists.txt | 1 + .../physics/tests/solid_multi_material.cpp | 121 ++++++++++++++++++ 3 files changed, 126 insertions(+), 4 deletions(-) create mode 100644 src/serac/physics/tests/solid_multi_material.cpp diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 3ba8818ab..562ad571c 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -877,7 +877,7 @@ class SolidMechanics, std::integer_se * @note This method must be called prior to completeSetup() */ template - void setMaterial(DependsOn, const MaterialType& material, + void setMaterial(DependsOn, const MaterialType& material, const Domain& domain, qdata_type qdata = EmptyQData) { static_assert(std::is_same_v || std::is_same_v, @@ -890,14 +890,14 @@ class SolidMechanics, std::integer_se // fact that the displacement, acceleration, and shape // fields are always-on and come first, so the `n`th // parameter will actually be argument `n + NUM_STATE_VARS` - std::move(material_functor), mesh_, qdata); + std::move(material_functor), domain, qdata); } /// @overload template - void setMaterial(const MaterialType& material, std::shared_ptr> qdata = EmptyQData) + void setMaterial(const MaterialType& material, const Domain& domain, std::shared_ptr> qdata = EmptyQData) { - setMaterial(DependsOn<>{}, material, qdata); + setMaterial(DependsOn<>{}, material, domain, qdata); } /** diff --git a/src/serac/physics/tests/CMakeLists.txt b/src/serac/physics/tests/CMakeLists.txt index 698378cd1..fb253d888 100644 --- a/src/serac/physics/tests/CMakeLists.txt +++ b/src/serac/physics/tests/CMakeLists.txt @@ -18,6 +18,7 @@ set(physics_serial_test_sources dynamic_solid_adjoint.cpp quasistatic_solid_adjoint.cpp finite_element_vector_set_over_domain.cpp + solid_multi_material.cpp ) serac_add_tests(SOURCES ${physics_serial_test_sources} diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp new file mode 100644 index 000000000..ef8654b3f --- /dev/null +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -0,0 +1,121 @@ +// Copyright (c) 2019-2024, 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.hpp" + +#include "axom/slic/core/SimpleLogger.hpp" +#include + +#include "serac/mesh/mesh_utils.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/infrastructure/initialize.hpp" +#include "serac/infrastructure/terminator.hpp" + +namespace serac { + +template +tensor average(std::vector >& positions) +{ + tensor total{}; + for (auto x : positions) { + total += x; + } + return total / double(positions.size()); +} + + +TEST(Solid, MultiMaterial) +{ + MPI_Barrier(MPI_COMM_WORLD); + + constexpr int p = 2; + constexpr int dim = 3; + int serial_refinement = 0; + int parallel_refinement = 0; + + // Create DataStore + axom::sidre::DataStore datastore; + serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial"); + + auto mesh = mesh::refineAndDistribute(buildCuboidMesh(8, 1, 1, 8.0, 1.0, 1.0), serial_refinement, parallel_refinement); + + const std::string mesh_tag{"mesh"}; + + auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); + + // _solver_params_start + serac::LinearSolverOptions linear_options{.linear_solver = LinearSolver::SuperLU}; + + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1.0e-12, + .absolute_tol = 1.0e-12, + .max_iterations = 5000, + .print_level = 1}; + + SolidMechanics solid(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, + GeometricNonlinearities::Off, "solid_mechanics", mesh_tag); + // _solver_params_end + + using Material = solid_mechanics::LinearIsotropic; + + constexpr double E_left = 1.0; + constexpr double nu = 0.25; + Material mat_left{.K = E_left/3.0/(1 - 2*nu), .G = 0.5*E_left/(1 + nu)}; + + constexpr double E_right = 2.0; + Material mat_right{.K = E_right/3.0/(1 - 2*nu), .G = 0.5*E_right/(1 + nu)}; +// using Hardening = solid_mechanics::LinearHardening; +// using MaterialB = solid_mechanics::J2SmallStrain; + +// Hardening hardening{.sigma_y = 10.0, .Hi = 0.1, .density = 1.0;}; +// Material mat{ +// .E = 2.0, // Young's modulus +// .nu = 0.25, // Poisson's ratio +// .hardening = hardening, +// .Hk = 0.0, // kinematic hardening constant +// .density = 1.0 // mass density +// }; + +// MaterialB::State initial_state{}; + + //auto qdata = solid_solver.createQuadratureDataBuffer(initial_state); + + auto is_in_left = [](std::vector> coords, int /* attribute */) { + return average(coords)[0] < 4.0; + }; + Domain left = Domain::ofElements(pmesh, is_in_left); + + auto is_in_right = [=](std::vector> coords, int attr){ return !is_in_left(coords, attr); }; + Domain right = Domain::ofElements(pmesh, is_in_right); + + solid.setMaterial(mat_left, left); + solid.setMaterial(mat_right, right); + + solid.setDisplacementBCs({2}, [](auto){ return 0.0; }, 2); + solid.setDisplacementBCs({3}, [](auto){ return 1.0; }, 0); + solid.setDisplacementBCs({5}, [](auto){ return 0.0; }, 0); + solid.setDisplacementBCs({6}, [](auto){ return 0.0; }, 1); + + solid.completeSetup(); + + // Perform the quasi-static solve + solid.advanceTimestep(1.0); + solid.outputStateToDisk("paraview"); +} + +} // namespace serac + +int main(int argc, char* argv[]) +{ + testing::InitGoogleTest(&argc, argv); + + serac::initialize(argc, argv); + + int result = RUN_ALL_TESTS(); + + serac::exitGracefully(result); +} From 4addc24dd9d20d8ba42e6911fd82d038cd576d8a Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 11 Oct 2024 17:04:07 -0700 Subject: [PATCH 03/17] Make test rigorous --- .../physics/tests/solid_multi_material.cpp | 56 ++++++++++++++++--- 1 file changed, 47 insertions(+), 9 deletions(-) diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index ef8654b3f..e5c986d70 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -30,6 +30,15 @@ tensor average(std::vector >& positions) TEST(Solid, MultiMaterial) { + /* + * Checks multi material case with the following uniaxial problem: + * MATERIAL 1 MATERIAL 2 + * E = 1 E = 2 + * u = 0 --------------------|-------------------- stress = 1 + * + * strain = 1 0.5 + * + */ MPI_Barrier(MPI_COMM_WORLD); constexpr int p = 2; @@ -41,7 +50,12 @@ TEST(Solid, MultiMaterial) axom::sidre::DataStore datastore; serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial"); - auto mesh = mesh::refineAndDistribute(buildCuboidMesh(8, 1, 1, 8.0, 1.0, 1.0), serial_refinement, parallel_refinement); + constexpr double L = 8.0; + constexpr double W = 1.0; + constexpr double H = 1.0; + constexpr double VOLUME = L*W*H; + + auto mesh = mesh::refineAndDistribute(buildCuboidMesh(8, 1, 1, L, W, H), serial_refinement, parallel_refinement); const std::string mesh_tag{"mesh"}; @@ -63,11 +77,12 @@ TEST(Solid, MultiMaterial) using Material = solid_mechanics::LinearIsotropic; constexpr double E_left = 1.0; - constexpr double nu = 0.25; - Material mat_left{.K = E_left/3.0/(1 - 2*nu), .G = 0.5*E_left/(1 + nu)}; + constexpr double nu_left = 0.125; + Material mat_left{.K = E_left/3.0/(1 - 2*nu_left), .G = 0.5*E_left/(1 + nu_left)}; - constexpr double E_right = 2.0; - Material mat_right{.K = E_right/3.0/(1 - 2*nu), .G = 0.5*E_right/(1 + nu)}; + constexpr double E_right = 2.0*E_left; + constexpr double nu_right = 2*nu_left; + Material mat_right{.K = E_right/3.0/(1 - 2*nu_right), .G = 0.5*E_right/(1 + nu_right)}; // using Hardening = solid_mechanics::LinearHardening; // using MaterialB = solid_mechanics::J2SmallStrain; @@ -85,7 +100,7 @@ TEST(Solid, MultiMaterial) //auto qdata = solid_solver.createQuadratureDataBuffer(initial_state); auto is_in_left = [](std::vector> coords, int /* attribute */) { - return average(coords)[0] < 4.0; + return average(coords)[0] < 0.5*L; }; Domain left = Domain::ofElements(pmesh, is_in_left); @@ -95,16 +110,39 @@ TEST(Solid, MultiMaterial) solid.setMaterial(mat_left, left); solid.setMaterial(mat_right, right); - solid.setDisplacementBCs({2}, [](auto){ return 0.0; }, 2); - solid.setDisplacementBCs({3}, [](auto){ return 1.0; }, 0); + constexpr double stress = 1.0; + Domain end_face = Domain::ofBoundaryElements(pmesh, by_attr(3)); + solid.setTraction(DependsOn<>{}, [stress](auto, auto n, auto){ return stress*n; }, end_face); + + solid.setDisplacementBCs({2}, [](auto){ return 0.0; }, 1); solid.setDisplacementBCs({5}, [](auto){ return 0.0; }, 0); - solid.setDisplacementBCs({6}, [](auto){ return 0.0; }, 1); + solid.setDisplacementBCs({1}, [](auto){ return 0.0; }, 2); solid.completeSetup(); // Perform the quasi-static solve solid.advanceTimestep(1.0); solid.outputStateToDisk("paraview"); + + // Define output functionals for verification + + constexpr double subdomain_volume = 0.5*VOLUME; + + auto average_strain_integrand = [](auto, auto, auto displacement) { + auto strain = get<1>(displacement); + return strain[0][0]/subdomain_volume; + }; + + Functional)> average_strain_left({&solid.displacement().space()}); + average_strain_left.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, + left); + + Functional)> average_strain_right({&solid.displacement().space()}); + average_strain_right.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, + right); + + EXPECT_NEAR(average_strain_left(solid.time(), solid.displacement()), stress/E_left, 1e-10); + EXPECT_NEAR(average_strain_right(solid.time(), solid.displacement()), stress/E_right, 1e-10); } } // namespace serac From 1381019f20d12219a178a853ad8fe65e3fd3db9a Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sun, 13 Oct 2024 07:34:50 -0700 Subject: [PATCH 04/17] Make domain optional, default to whole mesh --- src/serac/physics/solid_mechanics.hpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 562ad571c..2c9a0adaa 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -869,6 +869,8 @@ class SolidMechanics, std::integer_se * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, * 3>`) * + * @param optional_domain The subdomain which will use this material. Must be a domain of elements. If + * no argument is provided, it defaults to the entire mesh. * @param qdata the buffer of material internal variables at each quadrature point * * @pre MaterialType must have a public member variable `density` @@ -877,9 +879,11 @@ class SolidMechanics, std::integer_se * @note This method must be called prior to completeSetup() */ template - void setMaterial(DependsOn, const MaterialType& material, const Domain& domain, + void setMaterial(DependsOn, const MaterialType& material, + const std::optional& optional_domain = std::nullopt, qdata_type qdata = EmptyQData) { + Domain domain = (optional_domain) ? *optional_domain : EntireDomain(mesh_); static_assert(std::is_same_v || std::is_same_v, "invalid quadrature data provided in setMaterial()"); MaterialStressFunctor material_functor(material, geom_nonlin_); @@ -895,7 +899,9 @@ class SolidMechanics, std::integer_se /// @overload template - void setMaterial(const MaterialType& material, const Domain& domain, std::shared_ptr> qdata = EmptyQData) + void setMaterial(const MaterialType& material, + const std::optional& domain = std::nullopt, + std::shared_ptr> qdata = EmptyQData) { setMaterial(DependsOn<>{}, material, domain, qdata); } From 30d341b3255a393892899bf120a71d13d37a3148 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 11:19:34 -0700 Subject: [PATCH 05/17] Handle internal variables (quadrature data) --- src/serac/numerics/functional/geometry.hpp | 19 +++ src/serac/physics/solid_mechanics.hpp | 5 +- src/serac/physics/state/state_manager.hpp | 24 +++ .../physics/tests/solid_multi_material.cpp | 138 +++++++++++++++++- 4 files changed, 181 insertions(+), 5 deletions(-) diff --git a/src/serac/numerics/functional/geometry.hpp b/src/serac/numerics/functional/geometry.hpp index d56312c22..33bedf76b 100644 --- a/src/serac/numerics/functional/geometry.hpp +++ b/src/serac/numerics/functional/geometry.hpp @@ -2,6 +2,8 @@ #include "mfem.hpp" +#include "serac/numerics/functional/domain.hpp" + namespace serac { /** @@ -76,6 +78,23 @@ inline std::array geometry_counts(cons return counts; } +/** + * @brief count the number of elements of each geometry in a domain + * @param mesh the domain to count + */ +inline std::array geometry_counts(const Domain& domain) +{ + std::array counts{}; + + counts[mfem::Geometry::SEGMENT] = uint32_t(domain.edge_ids_.size()); + counts[mfem::Geometry::TRIANGLE] = uint32_t(domain.tri_ids_.size()); + counts[mfem::Geometry::SQUARE] = uint32_t(domain.quad_ids_.size()); + counts[mfem::Geometry::TETRAHEDRON] = uint32_t(domain.tet_ids_.size()); + counts[mfem::Geometry::CUBE] = uint32_t(domain.hex_ids_.size()); + + return counts; +} + /** * @brief count the number of boundary elements of each geometry in a mesh * @param mesh the mesh to count diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 2c9a0adaa..cd8f00e00 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -418,9 +418,10 @@ class SolidMechanics, std::integer_se * @return std::shared_ptr< QuadratureData > */ template - qdata_type createQuadratureDataBuffer(T initial_state) + qdata_type createQuadratureDataBuffer(T initial_state, const std::optional& optional_domain = std::nullopt) { - return StateManager::newQuadratureDataBuffer(mesh_tag_, order, dim, initial_state); + Domain domain = (optional_domain) ? *optional_domain : EntireDomain(mesh_); + return StateManager::newQuadratureDataBuffer(domain, order, dim, initial_state); } /** diff --git a/src/serac/physics/state/state_manager.hpp b/src/serac/physics/state/state_manager.hpp index acf21d0aa..3fc4c5e78 100644 --- a/src/serac/physics/state/state_manager.hpp +++ b/src/serac/physics/state/state_manager.hpp @@ -130,6 +130,30 @@ class StateManager { return std::make_shared>(elems, qpts_per_elem, initial_state); } + // make qdata by domain + template + static std::shared_ptr> newQuadratureDataBuffer(const Domain& domain, int order, int dim, + T initial_state) + { + int Q = order + 1; + + std::array elems = geometry_counts(domain); + std::array qpts_per_elem{}; + + std::vector geometries; + if (dim == 2) { + geometries = {mfem::Geometry::TRIANGLE, mfem::Geometry::SQUARE}; + } else { + geometries = {mfem::Geometry::TETRAHEDRON, mfem::Geometry::CUBE}; + } + + for (auto geom : geometries) { + qpts_per_elem[size_t(geom)] = uint32_t(num_quadrature_points(geom, Q)); + } + + return std::make_shared>(elems, qpts_per_elem, initial_state); + } + /** * @brief Checks if StateManager has a dual with the given name * @param name A string that uniquely identifies the name diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index e5c986d70..79c55ddde 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -78,11 +78,11 @@ TEST(Solid, MultiMaterial) constexpr double E_left = 1.0; constexpr double nu_left = 0.125; - Material mat_left{.K = E_left/3.0/(1 - 2*nu_left), .G = 0.5*E_left/(1 + nu_left)}; + Material mat_left{.density = 1.0, .K = E_left/3.0/(1 - 2*nu_left), .G = 0.5*E_left/(1 + nu_left)}; constexpr double E_right = 2.0*E_left; constexpr double nu_right = 2*nu_left; - Material mat_right{.K = E_right/3.0/(1 - 2*nu_right), .G = 0.5*E_right/(1 + nu_right)}; + Material mat_right{.density = 1.0, .K = E_right/3.0/(1 - 2*nu_right), .G = 0.5*E_right/(1 + nu_right)}; // using Hardening = solid_mechanics::LinearHardening; // using MaterialB = solid_mechanics::J2SmallStrain; @@ -128,7 +128,7 @@ TEST(Solid, MultiMaterial) constexpr double subdomain_volume = 0.5*VOLUME; - auto average_strain_integrand = [](auto, auto, auto displacement) { + auto average_strain_integrand = [subdomain_volume](auto, auto, auto displacement) { auto strain = get<1>(displacement); return strain[0][0]/subdomain_volume; }; @@ -145,6 +145,138 @@ TEST(Solid, MultiMaterial) EXPECT_NEAR(average_strain_right(solid.time(), solid.displacement()), stress/E_right, 1e-10); } +TEST(Solid, MultiMaterialWithState) +{ + /* + * Checks multi-material case with the following uniaxial problem: + * MATERIAL 1 MATERIAL 2 + * E = 2 E = 2 + * nu_1 nu_2 = 0 + * sigma_y = 0.5*stress + * Hi = E / 3.6 (linear hardening modulus) + * u = 0 --------------------|-------------------- stress = 1 + * + * The solution for the strain is: + * strain = stress/E | (Hi + E)/(Hi*E)*stress - sigma_y/Hi + * + * nu_1 is chosen so that the problem is uniaxial. That is, the lateral contraction in the + * elastic side exactly matches the contraction in the plastic side. For the values given + * above, this corresponds to nu_1 = 0.45. + * + * This problem checks that both materials models are correctly called for their subdomains, and + * that the internal variables are correctly indexed for the multi-material case. + * + */ + MPI_Barrier(MPI_COMM_WORLD); + + constexpr int p = 2; + constexpr int dim = 3; + int serial_refinement = 0; + int parallel_refinement = 0; + + // Create DataStore + axom::sidre::DataStore datastore; + serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial"); + + constexpr double L = 8.0; + constexpr double W = 1.0; + constexpr double H = 1.0; + constexpr double VOLUME = L*W*H; + + constexpr double applied_stress = 1.0; + + auto mesh = mesh::refineAndDistribute(buildCuboidMesh(8, 1, 1, L, W, H), serial_refinement, parallel_refinement); + + const std::string mesh_tag{"mesh"}; + + auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); + + // _solver_params_start + serac::LinearSolverOptions linear_options{.linear_solver = LinearSolver::SuperLU}; + + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1.0e-12, + .absolute_tol = 1.0e-12, + .max_iterations = 5000, + .print_level = 1}; + + SolidMechanics solid(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, + GeometricNonlinearities::Off, "solid_mechanics", mesh_tag); + // _solver_params_end + + auto is_in_left = [](std::vector> coords, int /* attribute */) { + return average(coords)[0] < 0.5*L; + }; + Domain left = Domain::ofElements(pmesh, is_in_left); + Domain right = EntireDomain(pmesh) - left; + + using Hardening = solid_mechanics::LinearHardening; + using MaterialRight = solid_mechanics::J2SmallStrain; + + constexpr double E_right = 2.0; + constexpr double nu_right = 0.0; + constexpr double Hi = E_right/3.6; + constexpr double sigma_y = 0.75*applied_stress; + + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi}; + MaterialRight mat_right{ + .E = E_right, // Young's modulus + .nu = nu_right, // Poisson's ratio + .hardening = hardening, + .Hk = 0.0, // kinematic hardening constant + .density = 1.0 // mass density + }; + + using MaterialLeft = solid_mechanics::LinearIsotropic; + + constexpr double E_left = 2.0; + constexpr double nu_left = 0.5*E_left/Hi*(1 - sigma_y/applied_stress); + MaterialLeft mat_left{.density = 1.0, .K = E_left/3.0/(1 - 2*nu_left), .G = 0.5*E_left/(1 + nu_left)}; + + MaterialRight::State initial_state{}; + auto qdata = solid.createQuadratureDataBuffer(initial_state, right); + + solid.setMaterial(mat_left, left); + solid.setMaterial(mat_right, right, qdata); + + Domain end_face = Domain::ofBoundaryElements(pmesh, by_attr(3)); + solid.setTraction(DependsOn<>{}, [applied_stress](auto, auto n, auto){ return applied_stress*n; }, end_face); + + solid.setDisplacementBCs({2}, [](auto){ return 0.0; }, 1); + solid.setDisplacementBCs({5}, [](auto){ return 0.0; }, 0); + solid.setDisplacementBCs({1}, [](auto){ return 0.0; }, 2); + + solid.completeSetup(); + + std::cout << "setup complete " << std::endl; + + // Perform the quasi-static solve + solid.advanceTimestep(1.0); + solid.outputStateToDisk("paraview"); + + // Define output functionals for verification + + constexpr double subdomain_volume = 0.5*VOLUME; + + auto average_strain_integrand = [subdomain_volume](auto, auto, auto displacement) { + auto strain = get<1>(displacement); + return strain[0][0]/subdomain_volume; + }; + + Functional)> average_strain_left({&solid.displacement().space()}); + average_strain_left.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, + left); + + Functional)> average_strain_right({&solid.displacement().space()}); + average_strain_right.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, + right); + + EXPECT_NEAR(average_strain_left(solid.time(), solid.displacement()), applied_stress/E_left, 1e-10); + + double exact = (E_right + Hi)/(E_right*Hi)*applied_stress - sigma_y/Hi; + EXPECT_NEAR(average_strain_right(solid.time(), solid.displacement()), exact, 1e-10); +} + } // namespace serac int main(int argc, char* argv[]) From e21c0c64d78a598834267ef907177396e1aa8063 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 16:28:58 -0700 Subject: [PATCH 06/17] FIll in missing values for geometry type counts --- src/serac/numerics/functional/geometry.hpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/serac/numerics/functional/geometry.hpp b/src/serac/numerics/functional/geometry.hpp index 33bedf76b..c652b1455 100644 --- a/src/serac/numerics/functional/geometry.hpp +++ b/src/serac/numerics/functional/geometry.hpp @@ -86,12 +86,11 @@ inline std::array geometry_counts(cons { std::array counts{}; - counts[mfem::Geometry::SEGMENT] = uint32_t(domain.edge_ids_.size()); - counts[mfem::Geometry::TRIANGLE] = uint32_t(domain.tri_ids_.size()); - counts[mfem::Geometry::SQUARE] = uint32_t(domain.quad_ids_.size()); - counts[mfem::Geometry::TETRAHEDRON] = uint32_t(domain.tet_ids_.size()); - counts[mfem::Geometry::CUBE] = uint32_t(domain.hex_ids_.size()); - + constexpr std::array geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE, mfem::Geometry::SQUARE, + mfem::Geometry::TETRAHEDRON, mfem::Geometry::CUBE}; + for (auto geom : geometries) { + counts[uint32_t(geom)] = uint32_t(domain.get(geom).size()); + } return counts; } From 9a3e61f99e9ea9034cc9ceb663b2723eff1b44dd Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 16:36:47 -0700 Subject: [PATCH 07/17] Ipmrpove general comments of tests --- src/serac/physics/tests/solid_multi_material.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 79c55ddde..5b51f8f1d 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -154,10 +154,12 @@ TEST(Solid, MultiMaterialWithState) * nu_1 nu_2 = 0 * sigma_y = 0.5*stress * Hi = E / 3.6 (linear hardening modulus) + * x = 0 x = L/2 x = L * u = 0 --------------------|-------------------- stress = 1 * * The solution for the strain is: - * strain = stress/E | (Hi + E)/(Hi*E)*stress - sigma_y/Hi + * strain = | stress/E, x in (0, L/2) + * | (Hi + E)/(Hi*E)*stress - sigma_y/Hi, x in (L/2, L) * * nu_1 is chosen so that the problem is uniaxial. That is, the lateral contraction in the * elastic side exactly matches the contraction in the plastic side. For the values given From 53c7e19c3fa3ce2d38ce8170b39f0713cb3718d5 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 17:54:54 -0700 Subject: [PATCH 08/17] Add another setMaterial overload to avoid breaking changes for anyone --- src/serac/physics/solid_mechanics.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index cd8f00e00..33e01539f 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -907,6 +907,14 @@ class SolidMechanics, std::integer_se setMaterial(DependsOn<>{}, material, domain, qdata); } + /// @overload + template + void setMaterial(const MaterialType& material, + std::shared_ptr> qdata) + { + setMaterial(DependsOn<>{}, material, EntireDomain(mesh_), qdata); + } + /** * @brief Set the underlying finite element state to a prescribed displacement * From ac7aa4bb4965bd6602e22d236936683b593490d6 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 18:02:15 -0700 Subject: [PATCH 09/17] Make style --- src/serac/physics/solid_mechanics.hpp | 8 +- .../physics/tests/solid_multi_material.cpp | 165 +++++++++--------- 2 files changed, 87 insertions(+), 86 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 33e01539f..69dd9cadc 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -882,7 +882,7 @@ class SolidMechanics, std::integer_se template void setMaterial(DependsOn, const MaterialType& material, const std::optional& optional_domain = std::nullopt, - qdata_type qdata = EmptyQData) + qdata_type qdata = EmptyQData) { Domain domain = (optional_domain) ? *optional_domain : EntireDomain(mesh_); static_assert(std::is_same_v || std::is_same_v, @@ -900,8 +900,7 @@ class SolidMechanics, std::integer_se /// @overload template - void setMaterial(const MaterialType& material, - const std::optional& domain = std::nullopt, + void setMaterial(const MaterialType& material, const std::optional& domain = std::nullopt, std::shared_ptr> qdata = EmptyQData) { setMaterial(DependsOn<>{}, material, domain, qdata); @@ -909,8 +908,7 @@ class SolidMechanics, std::integer_se /// @overload template - void setMaterial(const MaterialType& material, - std::shared_ptr> qdata) + void setMaterial(const MaterialType& material, std::shared_ptr> qdata) { setMaterial(DependsOn<>{}, material, EntireDomain(mesh_), qdata); } diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 5b51f8f1d..3c0edbeab 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -18,7 +18,7 @@ namespace serac { template -tensor average(std::vector >& positions) +tensor average(std::vector>& positions) { tensor total{}; for (auto x : positions) { @@ -27,17 +27,16 @@ tensor average(std::vector >& positions) return total / double(positions.size()); } - TEST(Solid, MultiMaterial) { /* * Checks multi material case with the following uniaxial problem: * MATERIAL 1 MATERIAL 2 - * E = 1 E = 2 + * E = 1 E = 2 * u = 0 --------------------|-------------------- stress = 1 - * + * * strain = 1 0.5 - * + * */ MPI_Barrier(MPI_COMM_WORLD); @@ -50,10 +49,10 @@ TEST(Solid, MultiMaterial) axom::sidre::DataStore datastore; serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial"); - constexpr double L = 8.0; - constexpr double W = 1.0; - constexpr double H = 1.0; - constexpr double VOLUME = L*W*H; + constexpr double L = 8.0; + constexpr double W = 1.0; + constexpr double H = 1.0; + constexpr double VOLUME = L * W * H; auto mesh = mesh::refineAndDistribute(buildCuboidMesh(8, 1, 1, L, W, H), serial_refinement, parallel_refinement); @@ -76,47 +75,51 @@ TEST(Solid, MultiMaterial) using Material = solid_mechanics::LinearIsotropic; - constexpr double E_left = 1.0; + constexpr double E_left = 1.0; constexpr double nu_left = 0.125; - Material mat_left{.density = 1.0, .K = E_left/3.0/(1 - 2*nu_left), .G = 0.5*E_left/(1 + nu_left)}; + Material mat_left{.density = 1.0, .K = E_left / 3.0 / (1 - 2 * nu_left), .G = 0.5 * E_left / (1 + nu_left)}; - constexpr double E_right = 2.0*E_left; - constexpr double nu_right = 2*nu_left; - Material mat_right{.density = 1.0, .K = E_right/3.0/(1 - 2*nu_right), .G = 0.5*E_right/(1 + nu_right)}; -// using Hardening = solid_mechanics::LinearHardening; -// using MaterialB = solid_mechanics::J2SmallStrain; + constexpr double E_right = 2.0 * E_left; + constexpr double nu_right = 2 * nu_left; + Material mat_right{.density = 1.0, .K = E_right / 3.0 / (1 - 2 * nu_right), .G = 0.5 * E_right / (1 + nu_right)}; + // using Hardening = solid_mechanics::LinearHardening; + // using MaterialB = solid_mechanics::J2SmallStrain; -// Hardening hardening{.sigma_y = 10.0, .Hi = 0.1, .density = 1.0;}; -// Material mat{ -// .E = 2.0, // Young's modulus -// .nu = 0.25, // Poisson's ratio -// .hardening = hardening, -// .Hk = 0.0, // kinematic hardening constant -// .density = 1.0 // mass density -// }; + // Hardening hardening{.sigma_y = 10.0, .Hi = 0.1, .density = 1.0;}; + // Material mat{ + // .E = 2.0, // Young's modulus + // .nu = 0.25, // Poisson's ratio + // .hardening = hardening, + // .Hk = 0.0, // kinematic hardening constant + // .density = 1.0 // mass density + // }; -// MaterialB::State initial_state{}; + // MaterialB::State initial_state{}; - //auto qdata = solid_solver.createQuadratureDataBuffer(initial_state); + // auto qdata = solid_solver.createQuadratureDataBuffer(initial_state); auto is_in_left = [](std::vector> coords, int /* attribute */) { - return average(coords)[0] < 0.5*L; + return average(coords)[0] < 0.5 * L; }; Domain left = Domain::ofElements(pmesh, is_in_left); - auto is_in_right = [=](std::vector> coords, int attr){ return !is_in_left(coords, attr); }; - Domain right = Domain::ofElements(pmesh, is_in_right); + auto is_in_right = [=](std::vector> coords, int attr) { return !is_in_left(coords, attr); }; + Domain right = Domain::ofElements(pmesh, is_in_right); solid.setMaterial(mat_left, left); solid.setMaterial(mat_right, right); - constexpr double stress = 1.0; - Domain end_face = Domain::ofBoundaryElements(pmesh, by_attr(3)); - solid.setTraction(DependsOn<>{}, [stress](auto, auto n, auto){ return stress*n; }, end_face); + constexpr double stress = 1.0; + Domain end_face = Domain::ofBoundaryElements(pmesh, by_attr(3)); + solid.setTraction( + DependsOn<>{}, [stress](auto, auto n, auto) { return stress * n; }, end_face); - solid.setDisplacementBCs({2}, [](auto){ return 0.0; }, 1); - solid.setDisplacementBCs({5}, [](auto){ return 0.0; }, 0); - solid.setDisplacementBCs({1}, [](auto){ return 0.0; }, 2); + solid.setDisplacementBCs( + {2}, [](auto) { return 0.0; }, 1); + solid.setDisplacementBCs( + {5}, [](auto) { return 0.0; }, 0); + solid.setDisplacementBCs( + {1}, [](auto) { return 0.0; }, 2); solid.completeSetup(); @@ -126,23 +129,21 @@ TEST(Solid, MultiMaterial) // Define output functionals for verification - constexpr double subdomain_volume = 0.5*VOLUME; + constexpr double subdomain_volume = 0.5 * VOLUME; auto average_strain_integrand = [subdomain_volume](auto, auto, auto displacement) { auto strain = get<1>(displacement); - return strain[0][0]/subdomain_volume; + return strain[0][0] / subdomain_volume; }; Functional)> average_strain_left({&solid.displacement().space()}); - average_strain_left.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, - left); + average_strain_left.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, left); Functional)> average_strain_right({&solid.displacement().space()}); - average_strain_right.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, - right); + average_strain_right.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, right); - EXPECT_NEAR(average_strain_left(solid.time(), solid.displacement()), stress/E_left, 1e-10); - EXPECT_NEAR(average_strain_right(solid.time(), solid.displacement()), stress/E_right, 1e-10); + EXPECT_NEAR(average_strain_left(solid.time(), solid.displacement()), stress / E_left, 1e-10); + EXPECT_NEAR(average_strain_right(solid.time(), solid.displacement()), stress / E_right, 1e-10); } TEST(Solid, MultiMaterialWithState) @@ -156,18 +157,18 @@ TEST(Solid, MultiMaterialWithState) * Hi = E / 3.6 (linear hardening modulus) * x = 0 x = L/2 x = L * u = 0 --------------------|-------------------- stress = 1 - * + * * The solution for the strain is: * strain = | stress/E, x in (0, L/2) * | (Hi + E)/(Hi*E)*stress - sigma_y/Hi, x in (L/2, L) - * + * * nu_1 is chosen so that the problem is uniaxial. That is, the lateral contraction in the * elastic side exactly matches the contraction in the plastic side. For the values given * above, this corresponds to nu_1 = 0.45. - * + * * This problem checks that both materials models are correctly called for their subdomains, and * that the internal variables are correctly indexed for the multi-material case. - * + * */ MPI_Barrier(MPI_COMM_WORLD); @@ -180,10 +181,10 @@ TEST(Solid, MultiMaterialWithState) axom::sidre::DataStore datastore; serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial"); - constexpr double L = 8.0; - constexpr double W = 1.0; - constexpr double H = 1.0; - constexpr double VOLUME = L*W*H; + constexpr double L = 8.0; + constexpr double W = 1.0; + constexpr double H = 1.0; + constexpr double VOLUME = L * W * H; constexpr double applied_stress = 1.0; @@ -207,46 +208,50 @@ TEST(Solid, MultiMaterialWithState) // _solver_params_end auto is_in_left = [](std::vector> coords, int /* attribute */) { - return average(coords)[0] < 0.5*L; + return average(coords)[0] < 0.5 * L; }; - Domain left = Domain::ofElements(pmesh, is_in_left); + Domain left = Domain::ofElements(pmesh, is_in_left); Domain right = EntireDomain(pmesh) - left; - using Hardening = solid_mechanics::LinearHardening; - using MaterialRight = solid_mechanics::J2SmallStrain; + using Hardening = solid_mechanics::LinearHardening; + using MaterialRight = solid_mechanics::J2SmallStrain; - constexpr double E_right = 2.0; + constexpr double E_right = 2.0; constexpr double nu_right = 0.0; - constexpr double Hi = E_right/3.6; - constexpr double sigma_y = 0.75*applied_stress; + constexpr double Hi = E_right / 3.6; + constexpr double sigma_y = 0.75 * applied_stress; - Hardening hardening{.sigma_y = sigma_y, .Hi = Hi}; + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi}; MaterialRight mat_right{ - .E = E_right, // Young's modulus - .nu = nu_right, // Poisson's ratio - .hardening = hardening, - .Hk = 0.0, // kinematic hardening constant - .density = 1.0 // mass density + .E = E_right, // Young's modulus + .nu = nu_right, // Poisson's ratio + .hardening = hardening, + .Hk = 0.0, // kinematic hardening constant + .density = 1.0 // mass density }; using MaterialLeft = solid_mechanics::LinearIsotropic; - constexpr double E_left = 2.0; - constexpr double nu_left = 0.5*E_left/Hi*(1 - sigma_y/applied_stress); - MaterialLeft mat_left{.density = 1.0, .K = E_left/3.0/(1 - 2*nu_left), .G = 0.5*E_left/(1 + nu_left)}; + constexpr double E_left = 2.0; + constexpr double nu_left = 0.5 * E_left / Hi * (1 - sigma_y / applied_stress); + MaterialLeft mat_left{.density = 1.0, .K = E_left / 3.0 / (1 - 2 * nu_left), .G = 0.5 * E_left / (1 + nu_left)}; MaterialRight::State initial_state{}; - auto qdata = solid.createQuadratureDataBuffer(initial_state, right); + auto qdata = solid.createQuadratureDataBuffer(initial_state, right); solid.setMaterial(mat_left, left); solid.setMaterial(mat_right, right, qdata); Domain end_face = Domain::ofBoundaryElements(pmesh, by_attr(3)); - solid.setTraction(DependsOn<>{}, [applied_stress](auto, auto n, auto){ return applied_stress*n; }, end_face); + solid.setTraction( + DependsOn<>{}, [applied_stress](auto, auto n, auto) { return applied_stress * n; }, end_face); - solid.setDisplacementBCs({2}, [](auto){ return 0.0; }, 1); - solid.setDisplacementBCs({5}, [](auto){ return 0.0; }, 0); - solid.setDisplacementBCs({1}, [](auto){ return 0.0; }, 2); + solid.setDisplacementBCs( + {2}, [](auto) { return 0.0; }, 1); + solid.setDisplacementBCs( + {5}, [](auto) { return 0.0; }, 0); + solid.setDisplacementBCs( + {1}, [](auto) { return 0.0; }, 2); solid.completeSetup(); @@ -258,28 +263,26 @@ TEST(Solid, MultiMaterialWithState) // Define output functionals for verification - constexpr double subdomain_volume = 0.5*VOLUME; + constexpr double subdomain_volume = 0.5 * VOLUME; auto average_strain_integrand = [subdomain_volume](auto, auto, auto displacement) { auto strain = get<1>(displacement); - return strain[0][0]/subdomain_volume; + return strain[0][0] / subdomain_volume; }; Functional)> average_strain_left({&solid.displacement().space()}); - average_strain_left.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, - left); + average_strain_left.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, left); Functional)> average_strain_right({&solid.displacement().space()}); - average_strain_right.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, - right); + average_strain_right.AddDomainIntegral(Dimension{}, DependsOn<0>{}, average_strain_integrand, right); - EXPECT_NEAR(average_strain_left(solid.time(), solid.displacement()), applied_stress/E_left, 1e-10); + EXPECT_NEAR(average_strain_left(solid.time(), solid.displacement()), applied_stress / E_left, 1e-10); - double exact = (E_right + Hi)/(E_right*Hi)*applied_stress - sigma_y/Hi; + double exact = (E_right + Hi) / (E_right * Hi) * applied_stress - sigma_y / Hi; EXPECT_NEAR(average_strain_right(solid.time(), solid.displacement()), exact, 1e-10); } -} // namespace serac +} // namespace serac int main(int argc, char* argv[]) { From 9297ff9fa631b2af4fb2000f7ba71a2eaea1a24b Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 18:07:18 -0700 Subject: [PATCH 10/17] Remove dead code --- .../physics/tests/solid_multi_material.cpp | 24 ++++--------------- 1 file changed, 4 insertions(+), 20 deletions(-) diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 3c0edbeab..906b4aff3 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -35,6 +35,7 @@ TEST(Solid, MultiMaterial) * E = 1 E = 2 * u = 0 --------------------|-------------------- stress = 1 * + * Solution: * strain = 1 0.5 * */ @@ -82,29 +83,12 @@ TEST(Solid, MultiMaterial) constexpr double E_right = 2.0 * E_left; constexpr double nu_right = 2 * nu_left; Material mat_right{.density = 1.0, .K = E_right / 3.0 / (1 - 2 * nu_right), .G = 0.5 * E_right / (1 + nu_right)}; - // using Hardening = solid_mechanics::LinearHardening; - // using MaterialB = solid_mechanics::J2SmallStrain; - - // Hardening hardening{.sigma_y = 10.0, .Hi = 0.1, .density = 1.0;}; - // Material mat{ - // .E = 2.0, // Young's modulus - // .nu = 0.25, // Poisson's ratio - // .hardening = hardening, - // .Hk = 0.0, // kinematic hardening constant - // .density = 1.0 // mass density - // }; - - // MaterialB::State initial_state{}; - - // auto qdata = solid_solver.createQuadratureDataBuffer(initial_state); auto is_in_left = [](std::vector> coords, int /* attribute */) { return average(coords)[0] < 0.5 * L; }; - Domain left = Domain::ofElements(pmesh, is_in_left); - - auto is_in_right = [=](std::vector> coords, int attr) { return !is_in_left(coords, attr); }; - Domain right = Domain::ofElements(pmesh, is_in_right); + Domain left = Domain::ofElements(pmesh, is_in_left); + Domain right = EntireDomain(pmesh) - left; solid.setMaterial(mat_left, left); solid.setMaterial(mat_right, right); @@ -179,7 +163,7 @@ TEST(Solid, MultiMaterialWithState) // Create DataStore axom::sidre::DataStore datastore; - serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial"); + serac::StateManager::initialize(datastore, "solid_mechanics_multimaterial_with_state"); constexpr double L = 8.0; constexpr double W = 1.0; From 9e01a2d4cd0bdbf3d0c7f0e68bb7174e3bc43d83 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 18:09:06 -0700 Subject: [PATCH 11/17] Fix broken doxygen comment --- src/serac/numerics/functional/geometry.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/numerics/functional/geometry.hpp b/src/serac/numerics/functional/geometry.hpp index c652b1455..0f98f4aed 100644 --- a/src/serac/numerics/functional/geometry.hpp +++ b/src/serac/numerics/functional/geometry.hpp @@ -80,7 +80,7 @@ inline std::array geometry_counts(cons /** * @brief count the number of elements of each geometry in a domain - * @param mesh the domain to count + * @param domain the domain to count */ inline std::array geometry_counts(const Domain& domain) { From 502b15faa96ea40a996225d18f972e8837c9df45 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 19:02:09 -0700 Subject: [PATCH 12/17] Add missing doxygen --- src/serac/physics/solid_mechanics.hpp | 1 + src/serac/physics/state/state_manager.hpp | 11 ++++++++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 69dd9cadc..789294bc1 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -415,6 +415,7 @@ class SolidMechanics, std::integer_se * * @tparam T the type to be created at each quadrature point * @param initial_state the value to be broadcast to each quadrature point + * @param optional_domain restricts the quadrature buffer to the given domain * @return std::shared_ptr< QuadratureData > */ template diff --git a/src/serac/physics/state/state_manager.hpp b/src/serac/physics/state/state_manager.hpp index 3fc4c5e78..35a3f4f75 100644 --- a/src/serac/physics/state/state_manager.hpp +++ b/src/serac/physics/state/state_manager.hpp @@ -130,7 +130,16 @@ class StateManager { return std::make_shared>(elems, qpts_per_elem, initial_state); } - // make qdata by domain + /** + * @overload + * + * @tparam T the type to be created at each quadrature point + * @param domain The domain to which this buffer will pertain + * @param order The order of the discretization of the displacement and velocity fields + * @param dim The spatial dimension of the mesh + * @param initial_state the value to be broadcast to each quadrature point + * @return shared pointer to quadrature data buffer + */ template static std::shared_ptr> newQuadratureDataBuffer(const Domain& domain, int order, int dim, T initial_state) From 0bdb257a3ba9c9e76536c6767b63363344775666 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 20:38:49 -0700 Subject: [PATCH 13/17] Add yet another overload of setMaterial to avoid breaking changes --- src/serac/physics/solid_mechanics.hpp | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 789294bc1..d48e458f7 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -871,8 +871,7 @@ class SolidMechanics, std::integer_se * values will change to `dual` numbers rather than `double`. (e.g. `tensor` becomes `tensor, * 3>`) * - * @param optional_domain The subdomain which will use this material. Must be a domain of elements. If - * no argument is provided, it defaults to the entire mesh. + * @param domain The subdomain which will use this material. Must be a domain of elements. * @param qdata the buffer of material internal variables at each quadrature point * * @pre MaterialType must have a public member variable `density` @@ -882,10 +881,9 @@ class SolidMechanics, std::integer_se */ template void setMaterial(DependsOn, const MaterialType& material, - const std::optional& optional_domain = std::nullopt, - qdata_type qdata = EmptyQData) + const Domain& domain, + qdata_type qdata = EmptyQData) { - Domain domain = (optional_domain) ? *optional_domain : EntireDomain(mesh_); static_assert(std::is_same_v || std::is_same_v, "invalid quadrature data provided in setMaterial()"); MaterialStressFunctor material_functor(material, geom_nonlin_); @@ -899,17 +897,26 @@ class SolidMechanics, std::integer_se std::move(material_functor), domain, qdata); } + /// @overload + template + void setMaterial(DependsOn, const MaterialType& material, + qdata_type qdata = EmptyQData) + { + setMaterial(DependsOn{}, material, EntireDomain(mesh_), qdata); + } + /// @overload template - void setMaterial(const MaterialType& material, const std::optional& domain = std::nullopt, + void setMaterial(const MaterialType& material, const Domain& domain, std::shared_ptr> qdata = EmptyQData) { setMaterial(DependsOn<>{}, material, domain, qdata); } /// @overload - template - void setMaterial(const MaterialType& material, std::shared_ptr> qdata) + template + void setMaterial(const MaterialType& material, + std::shared_ptr> qdata = EmptyQData) { setMaterial(DependsOn<>{}, material, EntireDomain(mesh_), qdata); } From e6d34c5b22b5b56c67cb9b74145eea0a6d9c877d Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 18 Oct 2024 20:40:10 -0700 Subject: [PATCH 14/17] Make style --- src/serac/physics/solid_mechanics.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index d48e458f7..d0974e3de 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -880,8 +880,7 @@ class SolidMechanics, std::integer_se * @note This method must be called prior to completeSetup() */ template - void setMaterial(DependsOn, const MaterialType& material, - const Domain& domain, + void setMaterial(DependsOn, const MaterialType& material, const Domain& domain, qdata_type qdata = EmptyQData) { static_assert(std::is_same_v || std::is_same_v, @@ -914,9 +913,8 @@ class SolidMechanics, std::integer_se } /// @overload - template - void setMaterial(const MaterialType& material, - std::shared_ptr> qdata = EmptyQData) + template + void setMaterial(const MaterialType& material, std::shared_ptr> qdata = EmptyQData) { setMaterial(DependsOn<>{}, material, EntireDomain(mesh_), qdata); } From f647b0d7983ed8e5acf1ff7ee8d719a370d345bc Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 21 Oct 2024 17:00:53 -0700 Subject: [PATCH 15/17] Add explicit template arguments to fix errors on ancient clang --- src/serac/numerics/functional/geometry.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/serac/numerics/functional/geometry.hpp b/src/serac/numerics/functional/geometry.hpp index 0f98f4aed..eb4da0b2a 100644 --- a/src/serac/numerics/functional/geometry.hpp +++ b/src/serac/numerics/functional/geometry.hpp @@ -86,8 +86,9 @@ inline std::array geometry_counts(cons { std::array counts{}; - constexpr std::array geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE, mfem::Geometry::SQUARE, - mfem::Geometry::TETRAHEDRON, mfem::Geometry::CUBE}; + constexpr std::array geometries = {mfem::Geometry::SEGMENT, mfem::Geometry::TRIANGLE, + mfem::Geometry::SQUARE, mfem::Geometry::TETRAHEDRON, + mfem::Geometry::CUBE}; for (auto geom : geometries) { counts[uint32_t(geom)] = uint32_t(domain.get(geom).size()); } From d62448a4ee0464447f892089ebe1f1ff17a0f7ad Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 25 Oct 2024 15:18:01 -0700 Subject: [PATCH 16/17] Remove redundant allocation function for quadrature data I think it's better to force an explicit domain. We don't need an overloaded versino for the whole mesh, since we have the EntireDomain(mesh) function. --- src/serac/physics/state/state_manager.hpp | 39 ++--------------------- 1 file changed, 2 insertions(+), 37 deletions(-) diff --git a/src/serac/physics/state/state_manager.hpp b/src/serac/physics/state/state_manager.hpp index 35a3f4f75..34207a8b7 100644 --- a/src/serac/physics/state/state_manager.hpp +++ b/src/serac/physics/state/state_manager.hpp @@ -99,43 +99,8 @@ class StateManager { * @brief Create a shared ptr to a quadrature data buffer for the given material type * * @tparam T the type to be created at each quadrature point - * @param mesh_tag The tag for the stored mesh used to construct the finite element state - * @param order The order of the discretization of the displacement and velocity fields - * @param dim The spatial dimension of the mesh - * @param initial_state the value to be broadcast to each quadrature point - * @return shared pointer to quadrature data buffer - */ - template - static std::shared_ptr> newQuadratureDataBuffer(const std::string& mesh_tag, int order, int dim, - T initial_state) - { - SLIC_ERROR_ROOT_IF(!hasMesh(mesh_tag), axom::fmt::format("Mesh tag '{}' not found in the data store", mesh_tag)); - - int Q = order + 1; - - std::array elems = geometry_counts(mesh(mesh_tag)); - std::array qpts_per_elem{}; - - std::vector geometries; - if (dim == 2) { - geometries = {mfem::Geometry::TRIANGLE, mfem::Geometry::SQUARE}; - } else { - geometries = {mfem::Geometry::TETRAHEDRON, mfem::Geometry::CUBE}; - } - - for (auto geom : geometries) { - qpts_per_elem[size_t(geom)] = uint32_t(num_quadrature_points(geom, Q)); - } - - return std::make_shared>(elems, qpts_per_elem, initial_state); - } - - /** - * @overload - * - * @tparam T the type to be created at each quadrature point - * @param domain The domain to which this buffer will pertain - * @param order The order of the discretization of the displacement and velocity fields + * @param domain The spatial domain over which to allocate the quadrature data + * @param order The order of the discretization of the primal fields * @param dim The spatial dimension of the mesh * @param initial_state the value to be broadcast to each quadrature point * @return shared pointer to quadrature data buffer From d6c8b36bc88591f763b14f02995d625130b7066d Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 25 Oct 2024 15:28:21 -0700 Subject: [PATCH 17/17] Address Chris's other comments --- src/serac/physics/tests/solid_multi_material.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 906b4aff3..3debb2c54 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -6,7 +6,6 @@ #include "serac/physics/solid_mechanics.hpp" -#include "axom/slic/core/SimpleLogger.hpp" #include #include "serac/mesh/mesh_utils.hpp" @@ -178,7 +177,6 @@ TEST(Solid, MultiMaterialWithState) auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); - // _solver_params_start serac::LinearSolverOptions linear_options{.linear_solver = LinearSolver::SuperLU}; serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, @@ -189,7 +187,6 @@ TEST(Solid, MultiMaterialWithState) SolidMechanics solid(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options, GeometricNonlinearities::Off, "solid_mechanics", mesh_tag); - // _solver_params_end auto is_in_left = [](std::vector> coords, int /* attribute */) { return average(coords)[0] < 0.5 * L;