time_discretized_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"
#include "smith/physics/mesh.hpp"
#include "smith/differentiable_numerics/field_state.hpp"
#include "smith/differentiable_numerics/time_integration_rule.hpp"
 
namespace smith {
 
template <int spatial_dim, typename OutputSpace, typename inputs = Parameters<>>
class TimeDiscretizedWeakForm;
 
template <int spatial_dim, typename OutputSpace, typename... InputSpaces>
class TimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>>
    : public FunctionalWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>> {
 public:
  using WeakFormT = FunctionalWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>>;  
 
  TimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr<Mesh> mesh,
                          const mfem::ParFiniteElementSpace& output_mfem_space,
                          const typename WeakFormT::SpacesT& input_mfem_spaces)
      : WeakFormT(physics_name, mesh, output_mfem_space, input_mfem_spaces)
  {
  }
 
  template <int... active_parameters, typename BodyIntegralType>
  void addBodyIntegral(DependsOn<active_parameters...> depends_on, std::string body_name, BodyIntegralType integrand)
  {
    const double* dt = &this->dt_;
    const size_t* cycle = &this->cycle_;
    WeakFormT::addBodyIntegral(depends_on, body_name, [dt, cycle, integrand](double t, auto X, auto... inputs) {
      TimeInfo time_info(t, *dt, *cycle);
      return integrand(time_info, X, inputs...);
    });
  }
 
  template <typename BodyForceType, int... all_active_parameters>
  void addBodyIntegralImpl(std::string body_name, BodyForceType body_integral,
                           std::integer_sequence<int, all_active_parameters...>)
  {
    addBodyIntegral(DependsOn<all_active_parameters...>{}, body_name, body_integral);
  }
 
  template <typename BodyForceType>
  void addBodyIntegral(std::string body_name, BodyForceType body_integral)
  {
    addBodyIntegralImpl(body_name, body_integral, std::make_integer_sequence<int, sizeof...(InputSpaces)>{});
  }
};
 
class SecondOrderTimeDiscretizedWeakForms {
 public:
  std::shared_ptr<WeakForm> time_discretized_weak_form;  
  std::shared_ptr<WeakForm>
      final_reaction_weak_form;  
};
 
template <int spatial_dim, typename OutputSpace, typename inputs = Parameters<>>
class SecondOrderTimeDiscretizedWeakForm;
 
template <int spatial_dim, typename OutputSpace, typename TrialInputSpace, typename... InputSpaces>
class SecondOrderTimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<TrialInputSpace, InputSpaces...>>
    : public SecondOrderTimeDiscretizedWeakForms {
 public:
  static constexpr int NUM_STATE_VARS = 4;  
 
  using TimeDiscretizedWeakFormT =
      TimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<TrialInputSpace, InputSpaces...>>;        
  using FinalReactionFormT = TimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>>;  
 
  SecondOrderTimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr<Mesh> mesh,
                                     ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule,
                                     const mfem::ParFiniteElementSpace& output_mfem_space,
                                     const typename TimeDiscretizedWeakFormT::SpacesT& input_mfem_spaces)
      : time_rule_(time_rule)
  {
    time_discretized_weak_form_ =
        std::make_shared<TimeDiscretizedWeakFormT>(physics_name, mesh, output_mfem_space, input_mfem_spaces);
    time_discretized_weak_form = time_discretized_weak_form_;
 
    typename TimeDiscretizedWeakFormT::SpacesT input_mfem_spaces_trial_removed(std::next(input_mfem_spaces.begin()),
                                                                               input_mfem_spaces.end());
    final_reaction_weak_form_ =
        std::make_shared<FinalReactionFormT>(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed);
    final_reaction_weak_form = final_reaction_weak_form_;
  }
 
  template <int... active_parameters, typename BodyIntegralType>
  void addBodyIntegral(DependsOn<active_parameters...> /*depends_on*/, std::string body_name,
                       BodyIntegralType integrand)
  {
    auto time_rule = time_rule_;
    time_discretized_weak_form_->addBodyIntegral(
        DependsOn<0, 1, 2, 3, NUM_STATE_VARS + active_parameters...>{}, body_name,
        [integrand, time_rule](const TimeInfo& t, auto X, auto U, auto U_old, auto U_dot_old, auto U_dot_dot_old,
                               auto... inputs) {
          return integrand(t, X, time_rule.value(t, U, U_old, U_dot_old, U_dot_dot_old),
                           time_rule.dot(t, U, U_old, U_dot_old, U_dot_dot_old),
                           time_rule.ddot(t, U, U_old, U_dot_old, U_dot_dot_old), inputs...);
        });
    final_reaction_weak_form_->addBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{},
                                               body_name, integrand);
  }
 
  template <typename BodyForceType>
  void addBodyIntegral(std::string body_name, BodyForceType body_integral)
  {
    addBodyIntegral(DependsOn<>{}, body_name, body_integral);
  }
 
 private:
  std::shared_ptr<TimeDiscretizedWeakFormT>
      time_discretized_weak_form_;  
  std::shared_ptr<FinalReactionFormT>
      final_reaction_weak_form_;  
 
  ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule_;  
};
 
}  // namespace smith