solid_mechanics_system.hpp

Go to the documentation of this file.

// Copyright (c) Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "smith/differentiable_numerics/field_store.hpp"
#include "smith/differentiable_numerics/nonlinear_block_solver.hpp"
#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp"
#include "smith/differentiable_numerics/state_advancer.hpp"
#include "smith/differentiable_numerics/multiphysics_time_integrator.hpp"
#include "smith/differentiable_numerics/time_integration_rule.hpp"
#include "smith/differentiable_numerics/time_discretized_weak_form.hpp"
#include "smith/differentiable_numerics/differentiable_physics.hpp"
#include "smith/physics/weak_form.hpp"
#include "smith/differentiable_numerics/system_base.hpp"
 
namespace smith {
 
template <int dim, int order, typename DisplacementTimeRule = ImplicitNewmarkSecondOrderTimeIntegrationRule,
          typename... parameter_space>
struct SolidMechanicsSystem : public SystemBase {
  static_assert(DisplacementTimeRule::num_states == 4, "SolidMechanicsSystem requires a 4-state time integration rule");
 
  using SolidWeakFormType =
      TimeDiscretizedWeakForm<dim, H1<order, dim>,
                              TimeRuleParams<DisplacementTimeRule, H1<order, dim>, parameter_space...>>;
 
  using CycleZeroSolidWeakFormType =
      TimeDiscretizedWeakForm<dim, H1<order, dim>,
                              Parameters<H1<order, dim>, H1<order, dim>, H1<order, dim>, parameter_space...>>;
 
  std::shared_ptr<SolidWeakFormType> solid_weak_form;  
  std::shared_ptr<CycleZeroSolidWeakFormType>
      cycle_zero_solid_weak_form;                        
  std::shared_ptr<DirichletBoundaryConditions> disp_bc;  
  std::shared_ptr<DisplacementTimeRule> disp_time_rule;  
 
  std::vector<FieldState> getStateFields() const
  {
    return {field_store->getField(prefix("displacement_solve_state")), field_store->getField(prefix("displacement")),
            field_store->getField(prefix("velocity")), field_store->getField(prefix("acceleration"))};
  }
 
  std::vector<FieldState> getOutputFieldStates() const
  {
    return {field_store->getField(prefix("displacement")), field_store->getField(prefix("velocity")),
            field_store->getField(prefix("acceleration"))};
  }
 
  std::vector<ReactionInfo> getReactionInfos() const
  {
    return {{prefix("reactions"), &field_store->getField(prefix("displacement")).get()->space()}};
  }
 
  std::unique_ptr<DifferentiablePhysics> createDifferentiablePhysics(std::string physics_name)
  {
    return std::make_unique<DifferentiablePhysics>(field_store->getMesh(), field_store->graph(),
                                                   field_store->getShapeDisp(), getStateFields(), getParameterFields(),
                                                   advancer, physics_name, getReactionInfos());
  }
 
  template <typename MaterialType>
  void setMaterial(const MaterialType& material, const std::string& domain_name)
  {
    auto captured_rule = disp_time_rule;
    solid_weak_form->addBodyIntegral(
        domain_name, [=](auto t_info, auto /*X*/, auto u, auto u_old, auto v_old, auto a_old, auto... params) {
          auto [u_current, v_current, a_current] = captured_rule->interpolate(t_info, u, u_old, v_old, a_old);
 
          typename MaterialType::State state;
          auto pk_stress = material(state, get<DERIVATIVE>(u_current), params...);
 
          return smith::tuple{get<VALUE>(a_current) * material.density, pk_stress};
        });
 
    // Add to cycle-zero weak form (at cycle 0, u and v are given, solve for a)
    cycle_zero_solid_weak_form->addBodyIntegral(
        domain_name, [=](auto /*t_info*/, auto /*X*/, auto u, auto /*v_old*/, auto a, auto... params) {
          typename MaterialType::State state;
          auto pk_stress = material(state, get<DERIVATIVE>(u), params...);
 
          return smith::tuple{get<VALUE>(a) * material.density, pk_stress};
        });
  }
 
  template <int... active_parameters, typename BodyForceType>
  void addBodyForce(DependsOn<active_parameters...> depends_on, const std::string& domain_name,
                    BodyForceType force_function)
  {
    auto captured_rule = disp_time_rule;
    solid_weak_form->addBodySource(
        depends_on, domain_name, [=](auto t_info, auto X, auto u, auto u_old, auto v_old, auto a_old, auto... params) {
          auto [u_current, v_current, a_current] = captured_rule->interpolate(t_info, u, u_old, v_old, a_old);
          return force_function(t_info.time(), X, u_current, v_current, a_current, params...);
        });
 
    addCycleZeroBodySourceImpl(
        domain_name,
        [=](auto t_info, auto X, auto u, auto v_old, auto a, auto... params) {
          return force_function(t_info.time(), X, u, v_old, a, params...);
        },
        std::make_index_sequence<3 + sizeof...(parameter_space)>{});
  }
 
  template <typename BodyForceType>
  void addBodyForce(const std::string& domain_name, BodyForceType force_function)
  {
    addBodyForceAllParams(domain_name, force_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{});
  }
 
  template <int... active_parameters, typename TractionType>
  void addTraction(DependsOn<active_parameters...> depends_on, const std::string& domain_name,
                   TractionType traction_function)
  {
    auto captured_rule = disp_time_rule;
    solid_weak_form->addBoundaryFlux(
        depends_on, domain_name,
        [=](auto t_info, auto X, auto n, auto u, auto u_old, auto v_old, auto a_old, auto... params) {
          auto [u_current, v_current, a_current] = captured_rule->interpolate(t_info, u, u_old, v_old, a_old);
          return traction_function(t_info.time(), X, n, u_current, v_current, a_current, params...);
        });
 
    addCycleZeroBoundaryFluxImpl(
        domain_name,
        [=](auto t_info, auto X, auto n, auto u, auto v_old, auto a, auto... params) {
          return traction_function(t_info.time(), X, n, u, v_old, a, params...);
        },
        std::make_index_sequence<3 + sizeof...(parameter_space)>{});
  }
 
  template <typename TractionType>
  void addTraction(const std::string& domain_name, TractionType traction_function)
  {
    addTractionAllParams(domain_name, traction_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{});
  }
 
  template <int... active_parameters, typename PressureType>
  void addPressure(DependsOn<active_parameters...> depends_on, const std::string& domain_name,
                   PressureType pressure_function)
  {
    auto captured_rule = disp_time_rule;
    solid_weak_form->addBoundaryIntegral(
        depends_on, domain_name, [=](auto t_info, auto X, auto u, auto u_old, auto v_old, auto a_old, auto... params) {
          auto u_current = captured_rule->value(t_info, u, u_old, v_old, a_old);
 
          auto x_current = X + u_current;
          auto n_deformed = cross(get<DERIVATIVE>(x_current));
          auto n_shape_norm = norm(cross(get<DERIVATIVE>(X)));
 
          auto pressure = pressure_function(t_info.time(), get<VALUE>(X), get<VALUE>(params)...);
 
          return pressure * n_deformed * (1.0 / n_shape_norm);
        });
 
    addCycleZeroBoundaryIntegralImpl(
        domain_name,
        [=](auto t_info, auto X, auto u, auto /*v_old*/, auto /*a*/, auto... params) {
          auto u_current = u;
 
          auto x_current = X + u_current;
          auto n_deformed = cross(get<DERIVATIVE>(x_current));
          auto n_shape_norm = norm(cross(get<DERIVATIVE>(X)));
 
          auto pressure = pressure_function(t_info.time(), get<VALUE>(X), get<VALUE>(params)...);
 
          return pressure * n_deformed * (1.0 / n_shape_norm);
        },
        std::make_index_sequence<3 + sizeof...(parameter_space)>{});
  }
 
  template <typename PressureType>
  void addPressure(const std::string& domain_name, PressureType pressure_function)
  {
    addPressureAllParams(domain_name, pressure_function, std::make_index_sequence<4 + sizeof...(parameter_space)>{});
  }
 
 private:
  template <typename BodyForceType, std::size_t... Is>
  void addBodyForceAllParams(const std::string& domain_name, BodyForceType force_function, std::index_sequence<Is...>)
  {
    addBodyForce(DependsOn<static_cast<int>(Is)...>{}, domain_name, force_function);
  }
 
  template <typename TractionType, std::size_t... Is>
  void addTractionAllParams(const std::string& domain_name, TractionType traction_function, std::index_sequence<Is...>)
  {
    addTraction(DependsOn<static_cast<int>(Is)...>{}, domain_name, traction_function);
  }
 
  template <typename PressureType, std::size_t... Is>
  void addPressureAllParams(const std::string& domain_name, PressureType pressure_function, std::index_sequence<Is...>)
  {
    addPressure(DependsOn<static_cast<int>(Is)...>{}, domain_name, pressure_function);
  }
 
  // Cycle-zero helpers: always use all-params DependsOn with the 3-state cycle-zero form count
  template <typename IntegrandType, std::size_t... Is>
  void addCycleZeroBodySourceImpl(const std::string& name, IntegrandType f, std::index_sequence<Is...>)
  {
    cycle_zero_solid_weak_form->addBodySource(DependsOn<static_cast<int>(Is)...>{}, name, f);
  }
 
  template <typename IntegrandType, std::size_t... Is>
  void addCycleZeroBoundaryFluxImpl(const std::string& name, IntegrandType f, std::index_sequence<Is...>)
  {
    cycle_zero_solid_weak_form->addBoundaryFlux(DependsOn<static_cast<int>(Is)...>{}, name, f);
  }
 
  template <typename IntegrandType, std::size_t... Is>
  void addCycleZeroBoundaryIntegralImpl(const std::string& name, IntegrandType f, std::index_sequence<Is...>)
  {
    cycle_zero_solid_weak_form->addBoundaryIntegral(DependsOn<static_cast<int>(Is)...>{}, name, f);
  }
};
 
template <int dim, int order, typename DisplacementTimeRule, typename... parameter_space>
SolidMechanicsSystem<dim, order, DisplacementTimeRule, parameter_space...> buildSolidMechanicsSystem(
    std::shared_ptr<Mesh> mesh, std::shared_ptr<CoupledSystemSolver> solver, DisplacementTimeRule disp_time_rule,
    std::string prepend_name = "", std::shared_ptr<CoupledSystemSolver> cycle_zero_solver = nullptr,
    FieldType<parameter_space>... parameter_types)
{
  auto field_store = std::make_shared<FieldStore>(mesh, 100);
 
  auto prefix = [&](const std::string& name) {
    if (prepend_name.empty()) {
      return name;
    }
    return prepend_name + "_" + name;
  };
 
  // Add shape displacement
  FieldType<H1<1, dim>> shape_disp_type(prefix("shape_displacement"));
  field_store->addShapeDisp(shape_disp_type);
 
  // Add displacement as independent (unknown) with time integration rule
  auto disp_time_rule_ptr = std::make_shared<DisplacementTimeRule>(disp_time_rule);
  FieldType<H1<order, dim>> disp_type(prefix("displacement_solve_state"));
  auto disp_bc = field_store->addIndependent(disp_type, disp_time_rule_ptr);
 
  // Add dependent fields for time integration history
  auto disp_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::VAL, prefix("displacement"));
  auto velo_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::DOT, prefix("velocity"));
  auto accel_old_type = field_store->addDependent(disp_type, FieldStore::TimeDerivative::DDOT, prefix("acceleration"));
 
  // Add parameters
  std::vector<FieldState> parameter_fields;
  (field_store->addParameter(FieldType<parameter_space>(prefix("param_" + parameter_types.name))), ...);
  (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...);
 
  using SystemType = SolidMechanicsSystem<dim, order, DisplacementTimeRule, parameter_space...>;
 
  // Create solid mechanics weak form (u, u_old, v_old, a_old)
  std::string force_name = prefix("solid_force");
  auto solid_weak_form = std::make_shared<typename SystemType::SolidWeakFormType>(
      force_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(),
      field_store->createSpaces(force_name, disp_type.name, disp_type, disp_old_type, velo_old_type, accel_old_type,
                                FieldType<parameter_space>(prefix("param_" + parameter_types.name))...));
 
  // Create cycle-zero weak form (u, v, a) for initial acceleration solve — 3 states, no u_old
  std::string cycle_zero_name = prefix("solid_reaction");
  auto cycle_zero_solid_weak_form = std::make_shared<typename SystemType::CycleZeroSolidWeakFormType>(
      cycle_zero_name, field_store->getMesh(), field_store->getField(accel_old_type.name).get()->space(),
      field_store->createSpaces(cycle_zero_name, accel_old_type.name, disp_type, velo_old_type, accel_old_type,
                                FieldType<parameter_space>(prefix("param_" + parameter_types.name))...));
 
  if (cycle_zero_solver == nullptr) {
    cycle_zero_solver = solver->singleBlockSolver(0);
  }
  SLIC_ERROR_IF(cycle_zero_solver == nullptr,
                "Could not derive a cycle-zero solver for block 0 from the provided solid mechanics solver.");
 
  // Build advancer
  std::vector<std::shared_ptr<WeakForm>> weak_forms{solid_weak_form};
  auto advancer = std::make_shared<MultiphysicsTimeIntegrator>(field_store, weak_forms, solver,
                                                               cycle_zero_solid_weak_form, cycle_zero_solver);
 
  return SystemType{{field_store, solver, advancer, parameter_fields, prepend_name},
                    solid_weak_form,
                    cycle_zero_solid_weak_form,
                    disp_bc,
                    disp_time_rule_ptr};
}
 
template <int dim, int order, typename DisplacementTimeRule, typename... parameter_space>
SolidMechanicsSystem<dim, order, DisplacementTimeRule, parameter_space...> buildSolidMechanicsSystem(
    std::shared_ptr<Mesh> mesh, std::shared_ptr<CoupledSystemSolver> solver, DisplacementTimeRule disp_time_rule,
    std::string prepend_name, FieldType<parameter_space>... parameter_types)
{
  return buildSolidMechanicsSystem<dim, order>(mesh, solver, disp_time_rule, std::move(prepend_name), nullptr,
                                               parameter_types...);
}
 
template <int dim, int order, typename DisplacementTimeRule, typename... parameter_space>
SolidMechanicsSystem<dim, order, DisplacementTimeRule, parameter_space...> buildSolidMechanicsSystem(
    std::shared_ptr<Mesh> mesh, std::shared_ptr<CoupledSystemSolver> solver, DisplacementTimeRule disp_time_rule,
    std::shared_ptr<CoupledSystemSolver> cycle_zero_solver, FieldType<parameter_space>... parameter_types)
{
  return buildSolidMechanicsSystem<dim, order>(mesh, solver, disp_time_rule, "", cycle_zero_solver, parameter_types...);
}
 
template <int dim, int order, typename DisplacementTimeRule, typename... parameter_space>
SolidMechanicsSystem<dim, order, DisplacementTimeRule, parameter_space...> buildSolidMechanicsSystem(
    std::shared_ptr<Mesh> mesh, std::shared_ptr<CoupledSystemSolver> solver, DisplacementTimeRule disp_time_rule,
    FieldType<parameter_space>... parameter_types)
{
  return buildSolidMechanicsSystem<dim, order>(mesh, solver, disp_time_rule, "", nullptr, parameter_types...);
}
 
}  // namespace smith