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/utilities.h>
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/tria.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/grid/tria.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/time_parameters.h>
#include <list>
#include <memory>
#include <set>
#include <string>
#include <utility>
#include <variant>
#include <vector>
namespace ImmersX
{
namespace LA
{
using namespace ImmersXLA;
} // namespace LA
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;
unsigned int initial_refinement = 2;
unsigned int n_refinement_cycles = 1;
std::set<dealii::types::boundary_id> dirichlet_ids{0};
std::set<dealii::types::boundary_id> neumann_ids;
std::string name_of_grid = "hyper_cube";
std::string arguments_for_grid = "-1: 1: false";
std::string triangulation_type = "distributed";
double density = 1.0;
double lame_mu = 1.0;
double lame_lambda = 1.0;
double damping_shear = 0.0;
double damping_bulk = 0.0;
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>>
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;
};
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 VectorType &
body_force_vector() const;
void
body_force_at_time(double time, VectorType &destination) const;
void
update_constraints(double time) 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:
using DistributedTriangulation =
dealii::parallel::distributed::Triangulation<dim, spacedim>;
using FullyDistributedTriangulation =
dealii::parallel::fullydistributed::Triangulation<dim, spacedim>;
using TriangulationVariant =
std::variant<DistributedTriangulation, FullyDistributedTriangulation>;
static TriangulationVariant
make_triangulation_storage(MPI_Comm mpi_communicator);
bool
uses_fully_distributed_triangulation() const;
void
assemble_body_force(double time);
void
add_initial_acceleration_constraint_rhs(VectorType &rhs) const;
void
assemble_backward_euler_system(const VectorType &previous_displacement,
const VectorType &previous_velocity,
double dt);
void
advance_one_timestep(double dt);
void
solve_backward_euler_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);
const ElastodynamicsParameters<dim, spacedim> ∥
MPI_Comm mpi_communicator;
dealii::ConditionalOStream pcout;
mutable dealii::TimerOutput computing_timer;
TriangulationVariant triangulation_storage;
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::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;
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;
MatrixType stiffness_matrix_storage;
MatrixType damping_matrix_storage;
MatrixType system_matrix_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