Program Listing for File elasticity.h¶
↰ Return to documentation for file (include/immersx/physics/elasticity.h)
// ---------------------------------------------------------------------
//
// Copyright (C) 2024 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 3.0 of the License, or (at your option) any later version. The
// full text of the license can be found in the file LICENSE.md at the top level
// of the ImmersX distribution.
//
// ---------------------------------------------------------------------
/* ---------------------------------------------------------------------
*/
#ifndef dealii_distributed_lagrange_multiplier_elasticity_h
#define dealii_distributed_lagrange_multiplier_elasticity_h
#include <deal.II/base/function.h>
#include <deal.II/base/parsed_convergence_table.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/block_linear_operator.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/linear_operator_tools.h>
#include <variant>
#define FORCE_USE_OF_TRILINOS
#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_function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.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_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_fe_field.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_tools_cache.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/physics/elasticity/standard_tensors.h>
// #include <deal.II/trilinos/parameter_acceptor.h>
#include <deal.II/lac/vector.h>
#ifdef Handle
# pragma push_macro("Handle")
# undef Handle
# define IMMERSX_RESTORE_OPENCASCADE_HANDLE
#endif
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/loop.h>
#include <deal.II/meshworker/scratch_data.h>
#include <deal.II/meshworker/simple.h>
#ifdef IMMERSX_RESTORE_OPENCASCADE_HANDLE
# pragma pop_macro("Handle")
# undef IMMERSX_RESTORE_OPENCASCADE_HANDLE
#endif
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/data_out_faces.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/vector_tools.h>
// #include <deal.II/numerics/matrix_tools.h>
#include <deal.II/opencascade/manifold_lib.h>
#include <deal.II/opencascade/utilities.h>
#include <immersx/coupling/inclusions.h>
#ifdef DEAL_II_WITH_VTK
# include <immersx/coupling/reduced_coupling.h>
#endif
#ifdef DEAL_II_WITH_OPENCASCADE
# include <TopoDS.hxx>
#endif
#include <deal.II/base/hdf5.h>
#include <immersx/physics/elasticity_problem_parameters.h>
#include <immersx/physics/material_properties.h>
#include <map>
#include <memory>
#include <string>
namespace ImmersX
{
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
template <int dim>
class RigidBodyMotion : public Function<dim>
{
public:
RigidBodyMotion(const unsigned int type_);
virtual double
value(const Point<dim> &p, const unsigned int component) const override;
private:
const unsigned int type;
};
template <int dim, int spacedim = dim>
class ElasticityProblem : public EnableObserverPointer
{
public:
ElasticityProblem(const ElasticityProblemParameters<dim, spacedim> &par);
void
make_grid();
void
setup_fe();
void
setup_dofs();
void
setup_constraints();
void
assemble_elasticity_system();
void
add_dynamic_constraint_rhs(
const AffineConstraints<double> &displacement_predictor_constraints,
const AffineConstraints<double> &velocity_predictor_constraints,
const AffineConstraints<double> &acceleration_constraints,
LA::MPI::Vector &rhs) const;
void
assemble_forcing_terms();
void
assemble_forcing_terms(const double evaluation_time);
void
compute_system_rhs();
void
compute_system_rhs(const double evaluation_time);
void
make_newmark_value_constraints(
const LA::MPI::Vector &values,
const AffineConstraints<double> &reference_constraints,
AffineConstraints<double> &value_constraints) const;
void
make_newmark_acceleration_constraints(
const LA::MPI::Vector &displacement_predictor,
AffineConstraints<double> &acceleration_constraints) const;
void
assemble_coupling();
void
run();
void
run_static();
void
run_quasistatic();
void
run_newmark();
IndexSet
assemble_coupling_sparsity(DynamicSparsityPattern &dsp);
void
solve();
void
solve_static();
void
solve_quasistatic();
void
solve_newmark();
void
refine_and_transfer();
void
refine_and_transfer_around_inclusions();
void
execute_actual_refine_and_transfer();
std::string
output_solution() const;
void
output_results() const;
void
print_parameters() const;
void
compute_internal_and_boundary_stress(bool) const;
bool
uses_fully_distributed_triangulation() const;
void
output_lambda() const;
std::string
output_stresses() const;
// void
// compute_face_stress();
const ElasticityProblemParameters<dim, spacedim> ∥
MPI_Comm mpi_communicator;
ConditionalOStream pcout;
mutable TimerOutput computing_timer;
using DistributedTriangulation =
parallel::distributed::Triangulation<spacedim>;
using FullyDistributedTriangulation =
parallel::fullydistributed::Triangulation<spacedim>;
using TriangulationVariant =
std::variant<DistributedTriangulation, FullyDistributedTriangulation>;
TriangulationVariant triangulation_storage;
parallel::TriangulationBase<spacedim> *tria;
std::unique_ptr<FiniteElement<spacedim>> fe;
Inclusions<spacedim> inclusions;
#ifdef DEAL_II_WITH_VTK
std::unique_ptr<ReducedCoupling<1, 2, spacedim, spacedim>>
tensor_product_coupling;
#endif
std::unique_ptr<Quadrature<spacedim>> quadrature;
std::unique_ptr<Quadrature<spacedim - 1>> face_quadrature_formula;
DoFHandler<spacedim> dh;
std::vector<IndexSet> owned_dofs;
std::vector<IndexSet> relevant_dofs;
AffineConstraints<double> constraints;
AffineConstraints<double>
inclusion_constraints;
AffineConstraints<double>
mean_value_constraints;
LA::MPI::SparseMatrix stiffness_matrix;
LA::MPI::SparseMatrix
newmark_matrix;
LA::MPI::SparseMatrix mass_matrix;
LA::MPI::SparseMatrix coupling_matrix;
LA::MPI::SparseMatrix damping_matrix;
LA::MPI::PreconditionAMG prec_A;
LA::MPI::PreconditionAMG prec_newmark;
LA::MPI::PreconditionAMG prec_C;
LA::MPI::PreconditionAMG prec_M;
LA::MPI::SparseMatrix inclusion_matrix;
LA::MPI::BlockVector solution;
LA::MPI::BlockVector velocity;
LA::MPI::BlockVector acceleration;
LA::MPI::BlockVector predictor;
LA::MPI::BlockVector corrector;
LA::MPI::BlockVector locally_relevant_solution;
LA::MPI::BlockVector force_rhs;
LA::MPI::BlockVector bc_rhs;
LA::MPI::BlockVector neumann_bc_rhs;
LA::MPI::BlockVector system_rhs;
std::vector<std::vector<BoundingBox<spacedim>>> global_bounding_boxes;
unsigned int cycle = 0;
unsigned int time_step = 0;
FEValuesExtractors::Vector displacement;
std::map<types::boundary_id, Tensor<1, spacedim>> forces;
std::map<types::boundary_id, Tensor<1, spacedim>>
average_displacements;
std::map<types::boundary_id, Tensor<1, spacedim>>
average_normals;
std::map<types::boundary_id, double>
areas;
TrilinosWrappers::MPI::Vector sigma_n;
double current_time = 0.0;
// Coupling representation queries.
bool
uses_tensor_product_coupling() const;
types::global_dof_index
n_multiplier_dofs() const;
unsigned int
n_multiplier_components_per_coupling_dof() const;
// Coupling-specific setup and assembly helpers.
void
setup_point_coupling_dofs();
void
setup_tensor_product_coupling_dofs();
void
assemble_point_coupling(const bool rebuild_matrices);
void
assemble_tensor_product_coupling();
bool coupling_matrices_assembled = false;
// Solver and output helpers
bool
has_immersed_coupling() const;
void
distribute_multiplier_solution(LA::MPI::Vector &lambda) const;
void
output_immersed_particles(const std::string &filename) const;
class Postprocessor;
};
} // namespace ImmersX
#endif