solid_weak_form.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/physics/functional_weak_form.hpp"
 
namespace smith {
 
template <int order, int dim, typename InputSpaces = Parameters<>>
class SolidWeakForm;
 
template <int order, int dim, typename... InputSpaces>
class SolidWeakForm<order, dim, Parameters<InputSpaces...>>
    : public FunctionalWeakForm<dim, H1<order, dim>,
                                Parameters<H1<order, dim>, H1<order, dim>, H1<order, dim>, InputSpaces...>> {
 public:
  using BaseWeakFormT = FunctionalWeakForm<dim, H1<order, dim>,
                                           Parameters<H1<order, dim>, H1<order, dim>, H1<order, dim>, InputSpaces...>>;
 
  template <typename T>
  using qdata_type = std::shared_ptr<QuadratureData<T>>;
 
  static constexpr int NUM_STATE_VARS = 3;
 
  enum STATE
  {
    DISPLACEMENT,
    VELOCITY,
    ACCELERATION,
    NUM_STATES
  };
 
  SolidWeakForm(std::string physics_name, std::shared_ptr<Mesh> mesh, const mfem::ParFiniteElementSpace& test_space,
                std::vector<const mfem::ParFiniteElementSpace*> parameter_fe_spaces = {})
      : BaseWeakFormT(physics_name, mesh, test_space, constructAllSpaces(test_space, parameter_fe_spaces))
  {
  }
 
  template <int... active_parameters, typename MaterialType, typename StateType = Empty>
  void setMaterial(DependsOn<active_parameters...>, std::string body_name, const MaterialType& material,
                   qdata_type<StateType> qdata = EmptyQData)
  {
    static_assert(std::is_same_v<StateType, Empty> || std::is_same_v<StateType, typename MaterialType::State>,
                  "invalid quadrature data provided in setMaterial()");
    MaterialStressFunctor<MaterialType> material_functor(material);
    BaseWeakFormT::weak_form_->AddDomainIntegral(
        Dimension<dim>{}, DependsOn<0, 2, active_parameters + NUM_STATE_VARS...>{}, std::move(material_functor),
        BaseWeakFormT::mesh_->domain(body_name), qdata);
 
    BaseWeakFormT::v_dot_weak_form_residual_->AddDomainIntegral(
        Dimension<dim>{}, DependsOn<0, 1, 3, active_parameters + 1 + NUM_STATE_VARS...>{},
        [material_functor](double t, auto X, auto state, auto V, auto... params) {
          auto flux = material_functor(t, X, state, params...);
          return smith::inner(get<VALUE>(V), get<VALUE>(flux)) +
                 smith::inner(get<DERIVATIVE>(V), get<DERIVATIVE>(flux));
        },
        BaseWeakFormT::mesh_->domain(body_name), qdata);
  }
 
  template <typename MaterialType, typename StateType = Empty>
  void setMaterial(std::string body_name, const MaterialType& material,
                   std::shared_ptr<QuadratureData<StateType>> qdata = EmptyQData)
  {
    setMaterial(DependsOn<>{}, body_name, material, qdata);
  }
 
  template <int... active_parameters, typename MaterialType, typename StateType = Empty>
  void setRateMaterial(DependsOn<active_parameters...>, std::string body_name, const MaterialType& material,
                       qdata_type<StateType> qdata = EmptyQData)
  {
    static_assert(std::is_same_v<StateType, Empty> || std::is_same_v<StateType, typename MaterialType::State>,
                  "invalid quadrature data provided in setMaterial()");
    RateMaterialStressFunctor<MaterialType> material_functor(material, &this->dt_);
    BaseWeakFormT::weak_form_->AddDomainIntegral(
        Dimension<dim>{}, DependsOn<0, 1, 2, active_parameters + NUM_STATE_VARS...>{}, std::move(material_functor),
        BaseWeakFormT::mesh_->domain(body_name), qdata);
 
    BaseWeakFormT::v_dot_weak_form_residual_->AddDomainIntegral(
        Dimension<dim>{}, DependsOn<0, 1, 2, 3, active_parameters + 1 + NUM_STATE_VARS...>{},
        [material_functor, qdata](double t, auto X, auto state, auto V, auto... params) {
          auto flux = material_functor(t, X, state, params...);
          return smith::inner(get<VALUE>(V), get<VALUE>(flux)) +
                 smith::inner(get<DERIVATIVE>(V), get<DERIVATIVE>(flux));
        },
        BaseWeakFormT::mesh_->domain(body_name), qdata);
  }
 
  template <typename MaterialType, typename StateType = Empty>
  void setRateMaterial(std::string body_name, const MaterialType& material,
                       std::shared_ptr<QuadratureData<StateType>> qdata = EmptyQData)
  {
    setRateMaterial(DependsOn<>{}, body_name, material, qdata);
  }
 
  template <int... active_parameters, typename PressureType>
  void addPressure(DependsOn<active_parameters...>, std::string boundary_name, PressureType pressure_function)
  {
    BaseWeakFormT::weak_form_->AddBoundaryIntegral(
        Dimension<dim - 1>{}, DependsOn<0, active_parameters + NUM_STATE_VARS...>{},
        [pressure_function](double t, auto X, auto displacement, auto... params) {
          // Calculate the position and normal in the shape perturbed deformed configuration
          auto x = X + displacement;
          auto n = cross(get<DERIVATIVE>(x));
          // smith::Functional's boundary integrals multiply the q-function output by
          // norm(cross(dX_dxi)) at that quadrature point, but if we impose a shape displacement
          // then that weight needs to be corrected. The new weight should be
          // norm(cross(dX_dxi + du_dxi + dp_dxi)) where u is displacement and p is shape displacement. This implies:
          //
          //   pressure * normalize(normal_new) * w_new
          // = pressure * normalize(normal_new) * (w_new / w_old) * w_old
          // = pressure * normalize(normal_new) * (norm(normal_new) / norm(normal_old)) * w_old
          // = pressure * (normal_new / norm(normal_new)) * (norm(normal_new) / norm(normal_old)) * w_old
          // = pressure * (normal_new / norm(normal_old)) * w_old
 
          // We always query the pressure function in the undeformed configuration
          return pressure_function(t, get<VALUE>(X), params...) * (n / norm(cross(get<DERIVATIVE>(X))));
        },
        BaseWeakFormT::mesh_->domain(boundary_name));
 
    BaseWeakFormT::v_dot_weak_form_residual_->AddBoundaryIntegral(
        Dimension<dim - 1>{}, DependsOn<0, 1, active_parameters + 1 + NUM_STATE_VARS...>{},
        [pressure_function](double t, auto X, auto V, auto displacement, auto... params) {
          auto x = X + displacement;
          auto n = cross(get<DERIVATIVE>(x));
          auto pressure = pressure_function(t, get<VALUE>(X), params...) * (n / norm(cross(get<DERIVATIVE>(X))));
          return inner(get<VALUE>(V), pressure);
        },
        BaseWeakFormT::mesh_->domain(boundary_name));
  }
 
  template <typename PressureType>
  void addPressure(std::string boundary_name, PressureType pressure_function)
  {
    addPressure(DependsOn<>{}, boundary_name, pressure_function);
  }
 
 protected:
  std::vector<const mfem::ParFiniteElementSpace*> constructAllSpaces(
      const mfem::ParFiniteElementSpace& state_space, const std::vector<const mfem::ParFiniteElementSpace*>& spaces)
  {
    std::vector<const mfem::ParFiniteElementSpace*> all_spaces{&state_space, &state_space, &state_space};
    for (auto& s : spaces) {
      all_spaces.push_back(s);
    }
    return all_spaces;
  }
 
  template <typename Material>
  struct MaterialStressFunctor {
    MaterialStressFunctor(Material material) : material_(material) {}
 
    Material material_;
 
    template <typename X, typename State, typename Displacement, typename Acceleration, typename... Params>
    auto SMITH_HOST_DEVICE operator()(double, X, State& state, Displacement displacement, Acceleration acceleration,
                                      Params... params) const
    {
      auto du_dX = get<DERIVATIVE>(displacement);
      auto d2u_dt2 = get<VALUE>(acceleration);
      auto stress = material_.pkStress(state, du_dX, params...);
      return smith::tuple{material_.density(params...) * d2u_dt2, stress};
    }
  };
 
  template <typename Material>
  struct RateMaterialStressFunctor {
    RateMaterialStressFunctor(Material material, const double* dt) : material_(material), dt_(dt) {}
 
    Material material_;
 
    const double* dt_;
 
    template <typename X, typename State, typename Displacement, typename Velocity, typename Acceleration,
              typename... Params>
    auto SMITH_HOST_DEVICE operator()(double /*t*/, X, State& state, Displacement displacement, Velocity velocity,
                                      Acceleration acceleration, Params... params) const
    {
      auto du_dX = get<DERIVATIVE>(displacement);
      auto dv_dX = get<DERIVATIVE>(velocity);
      auto d2u_dt2 = get<VALUE>(acceleration);
      auto stress = material_.pkStress(*dt_, state, du_dX, dv_dX, params...);
      return smith::tuple{material_.density(params...) * d2u_dt2, stress};
    }
  };
};
 
template <int order, int dim, typename... ParameterSpaces>
auto create_solid_weak_form(const std::string& physics_name, std::shared_ptr<smith::Mesh> mesh,
                            const std::vector<smith::FiniteElementState*>& states,  // u, v, a, e
                            const std::vector<smith::FiniteElementState*>& params)
{
  enum FieldNumbering
  {
    DISP,
    VELO,
    ACCEL,
  };
 
  std::vector<const mfem::ParFiniteElementSpace*> parameter_fe_spaces;
 
  if constexpr (sizeof...(ParameterSpaces) > 0) {
    for_constexpr<sizeof...(ParameterSpaces)>([&](auto i) { parameter_fe_spaces.push_back(&params[i]->space()); });
  }
 
  using WeakFormT = SolidWeakForm<order, dim, Parameters<ParameterSpaces...>>;
 
  return std::make_shared<WeakFormT>(physics_name, mesh, states[DISP]->space(), parameter_fe_spaces);
}
 
}  // namespace smith