Program Listing for File elastodynamics.h¶
↰ Return to documentation for file (include/immersx/physics/elastodynamics.h)
// ---------------------------------------------------------------------
//
// Copyright (C) 2026 by Luca Heltai
//
// This file is part of the ImmersX application, based on
// the deal.II library.
//
// The ImmersX application is free software; you can use
// it, redistribute it, and/or modify it under the terms of the Apache-2.0
// License WITH LLVM-exception as published by the Free Software Foundation;
// either version 2.0 of the License, or (at your option) any later version.
// The full text of the license can be found in the LICENSE.md file at the top
// level of the ImmersX distribution.
//
// ---------------------------------------------------------------------
#ifndef immersx_elastodynamics_h
#define immersx_elastodynamics_h
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/parsed_convergence_table.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/tria_base.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/numerics/data_out.h>
#include <immersx/algebra/linear_algebra.h>
#include <immersx/core/domain.h>
#include <immersx/core/robin_boundary.h>
#include <immersx/core/second_order_preconditioner.h>
#include <immersx/core/time_parameters.h>
#include <list>
#include <map>
#include <memory>
#include <optional>
#include <set>
#include <string>
#include <utility>
#include <vector>
namespace ImmersX
{
namespace LA
{
using namespace ImmersXLA;
} // namespace LA
struct ElastodynamicsMaterialCoefficients
{
double density = 1.0;
double lame_mu = 1.0;
double lame_lambda = 1.0;
double rayleigh_alpha = 0.0;
double rayleigh_beta = 0.0;
double viscosity_eta = 0.0;
};
class ElastodynamicsMaterialParameters : public dealii::ParameterAcceptor
{
public:
explicit ElastodynamicsMaterialParameters(const std::string &subsection);
ElastodynamicsMaterialCoefficients coefficients;
};
template <int dim, int spacedim = dim>
class ElastodynamicsParameters : public dealii::ParameterAcceptor
{
private:
std::unique_ptr<TimeIntervalParameters> owned_time_parameters;
std::unique_ptr<FixedStepParameters> owned_fixed_step_parameters;
std::unique_ptr<IDAParameters> owned_ida_parameters;
public:
explicit ElastodynamicsParameters(
const std::string &subsection = "/Elastodynamics/",
TimeIntervalParameters *shared_time_parameters = nullptr,
FixedStepParameters *shared_fixed_step_parameters = nullptr,
IDAParameters *shared_ida_parameters = nullptr);
std::string output_directory = ".";
std::string output_name = "elastodynamics";
unsigned int fe_degree = 1;
DomainParameters<dim, spacedim> domain_parameters;
unsigned int n_refinement_cycles = 1;
std::set<dealii::types::boundary_id> dirichlet_ids{0};
std::set<dealii::types::boundary_id> neumann_ids;
std::set<dealii::types::boundary_id> robin_ids;
ElastodynamicsMaterialParameters default_material;
std::map<dealii::types::material_id, std::string>
material_tags_by_material_id;
std::map<std::string, std::unique_ptr<ElastodynamicsMaterialParameters>>
material_parameters_by_tag;
ElastodynamicsMaterialCoefficients
get_material_properties(dealii::types::material_id material_id) const;
TimeIntervalParameters &time_parameters;
FixedStepParameters &fixed_step_parameters;
IDAParameters &ida_parameters;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
body_force;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
displacement_boundary;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
neumann_boundary;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
robin_coefficient;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
robin_rhs;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
velocity_boundary;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
initial_displacement;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
initial_velocity;
mutable dealii::ParameterAcceptorProxy<
dealii::Functions::ParsedFunction<spacedim>>
exact_solution;
mutable dealii::ParameterAcceptorProxy<dealii::ReductionControl>
solver_control;
mutable dealii::ParsedConvergenceTable convergence_table;
private:
template <int, int>
friend class ElastodynamicsSolver;
bool robin_coefficient_depends_on_time = true;
bool robin_rhs_depends_on_time = true;
bool body_force_depends_on_time = true;
bool neumann_boundary_depends_on_time = true;
bool body_force_is_zero = false;
bool neumann_boundary_is_zero = false;
bool displacement_boundary_depends_on_time = true;
bool velocity_boundary_depends_on_time = true;
};
template <int dim, int spacedim = dim>
class ElastodynamicsSolver : public dealii::EnableObserverPointer
{
static_assert(dim >= 1 && dim <= spacedim && spacedim <= 3,
"ElastodynamicsSolver requires 1 <= dim <= spacedim <= 3.");
public:
using VectorType = LA::MPI::Vector;
using MatrixType = LA::MPI::SparseMatrix;
explicit ElastodynamicsSolver(
const ElastodynamicsParameters<dim, spacedim> &par);
void
make_grid();
void
setup_fe();
void
setup_system();
void
assemble_operators();
void
set_initial_conditions();
void
refine_global();
void
initial_acceleration(VectorType &acceleration) const;
void
advance_one_timestep();
void
solve();
void
output_results() const;
void
compute_error() const;
void
run();
void
set_displacement(const VectorType &new_displacement);
void
set_velocity(const VectorType &new_velocity);
void
accept_state(const VectorType &new_displacement,
const VectorType &new_velocity,
double time,
unsigned int step_number);
dealii::types::global_dof_index
n_dofs() const;
bool
state_is_finite() const;
const dealii::parallel::TriangulationBase<dim, spacedim> &
triangulation() const;
const dealii::FiniteElement<dim, spacedim> &
fe() const;
const dealii::Mapping<dim, spacedim> &
mapping() const;
const dealii::DoFHandler<dim, spacedim> &
dof_handler() const;
const dealii::AffineConstraints<double> &
constraints() const;
const dealii::AffineConstraints<double> &
velocity_constraints() const;
const dealii::IndexSet &
locally_owned_dofs() const;
const dealii::IndexSet &
locally_relevant_dofs() const;
const MatrixType &
mass_matrix() const;
const MatrixType &
stiffness_matrix() const;
const MatrixType &
damping_matrix() const;
const MatrixType &
homogeneous_mass_matrix() const;
const MatrixType &
homogeneous_stiffness_matrix() const;
const MatrixType &
homogeneous_damping_matrix() const;
const VectorType &
body_force_vector() const;
void
body_force_at_time(double time, VectorType &destination) const;
void
update_constraints(double time) const;
void
update_boundary_operators(double time) const;
struct AssemblyStatistics
{
unsigned long spatial_operators = 0;
unsigned long robin_operators = 0;
unsigned long body_force = 0;
unsigned long constraints = 0;
unsigned long time_system_matrix = 0;
};
const AssemblyStatistics &
assembly_statistics() const;
const MatrixType &
system_matrix() const;
const VectorType &
system_rhs() const;
const VectorType &
displacement() const;
const VectorType &
velocity() const;
double
current_time() const;
double
time_step() const;
unsigned int
time_step_number() const;
private:
void
assemble_body_force(double time);
void
update_constraints(double time, bool update_velocity) const;
void
prepare_time_integration_matrix(double dt, double theta);
void
assemble_time_integration_rhs(const VectorType &previous_displacement,
const VectorType &previous_velocity,
const VectorType &previous_body_force,
double dt,
double theta);
void
assemble_backward_euler_system(const VectorType &previous_displacement,
const VectorType &previous_velocity,
double dt);
void
assemble_trapezoidal_system(const VectorType &previous_displacement,
const VectorType &previous_velocity,
const VectorType &previous_body_force,
double dt);
void
advance_one_timestep(double dt);
void
solve_time_integration_system();
void
update_locally_relevant_state();
void
run_time_integration();
static void
copy_constraints(const dealii::AffineConstraints<double> &source,
dealii::AffineConstraints<double> &target,
dealii::types::global_dof_index shift);
void
rebuild_combined_constraints() const;
void
update_discrete_velocity_constraints(
const dealii::AffineConstraints<double> &previous_displacement,
double previous_time,
double current_time) const;
void
update_trapezoidal_velocity_constraints(
const dealii::AffineConstraints<double> &previous_displacement,
const dealii::AffineConstraints<double> &previous_velocity,
double previous_time,
double current_time) const;
const ElastodynamicsParameters<dim, spacedim> ∥
MPI_Comm mpi_communicator;
dealii::ConditionalOStream pcout;
mutable dealii::TimerOutput computing_timer;
Domain<dim, spacedim> domain;
dealii::parallel::TriangulationBase<dim, spacedim> *tria;
std::unique_ptr<dealii::FiniteElement<dim, spacedim>> fe_storage;
std::unique_ptr<dealii::Quadrature<dim>> quadrature;
std::unique_ptr<dealii::Quadrature<dim - 1>> face_quadrature;
std::unique_ptr<dealii::Mapping<dim, spacedim>> mapping_storage;
dealii::DoFHandler<dim, spacedim> dh;
dealii::IndexSet owned_dofs;
dealii::IndexSet relevant_dofs;
dealii::IndexSet combined_owned_dofs;
dealii::IndexSet combined_relevant_dofs;
dealii::AffineConstraints<double> hanging_node_constraints_storage;
mutable dealii::AffineConstraints<double> displacement_constraints_storage;
mutable dealii::AffineConstraints<double> velocity_constraints_storage;
mutable dealii::AffineConstraints<double> combined_constraints_storage;
MatrixType mass_matrix_storage;
mutable MatrixType stiffness_matrix_storage;
MatrixType damping_matrix_storage;
MatrixType system_matrix_storage;
MatrixType homogeneous_mass_matrix_storage;
mutable MatrixType homogeneous_stiffness_matrix_storage;
mutable std::optional<double> robin_operator_time;
mutable VectorType cached_robin_force;
mutable std::optional<double> robin_force_time;
VectorType previous_stiffness_action;
MatrixType homogeneous_damping_matrix_storage;
mutable std::unique_ptr<
SecondOrderBlockPreconditioner<VectorType, MatrixType>>
time_preconditioner;
bool spatial_operators_ready = false;
mutable std::optional<double> system_matrix_time_step;
double system_matrix_theta = 0.;
mutable std::optional<double> displacement_constraints_time;
mutable std::optional<double> velocity_constraints_time;
mutable bool velocity_constraints_are_boundary_data = false;
mutable bool combined_constraints_initialized = false;
mutable VectorType cached_body_force;
mutable VectorType cached_volume_force;
mutable VectorType cached_traction_force;
mutable std::optional<double> volume_force_cache_time;
mutable std::optional<double> traction_force_cache_time;
mutable AssemblyStatistics assembly_statistics_storage;
VectorType body_force_storage;
VectorType displacement_storage;
VectorType velocity_storage;
VectorType locally_relevant_displacement;
VectorType locally_relevant_velocity;
VectorType system_rhs_storage;
mutable std::vector<std::pair<double, std::string>> cycles_and_solutions;
double current_time_storage = 0.0;
double current_time_step = 0.0;
unsigned int time_step_number_storage = 0;
unsigned int refinement_cycle_storage = 0;
};
} // namespace ImmersX
#endif // immersx_elastodynamics_h