Template Class ElasticityProblem

Inheritance Relationships

Base Type

  • public Subscriptor

Class Documentation

template<int dim, int spacedim = dim>
class ElasticityProblem : public Subscriptor

Finite element model for (optionally) time-dependent elasticity, coupled to immersed 1D/0D inclusions via reduced Lagrange multipliers.

This class assembles a small set of global sparse matrices and then solves either:

  • a pure displacement problem (no inclusions), or

  • a saddle-point problem for displacement and Lagrange multipliers, or

  • a displacement problem driven by inclusion “pressure” data (when ElasticityProblemParameters::pressure_coupling is enabled).

The intent of this documentation is to state exactly which continuous and discrete equations correspond to the current implementation in source/elasticity.cc (primarily assemble_elasticity_system(), assemble_coupling(), and solve()), and to list the constants and options that select between different model variants.

Unknowns, discrete spaces, and blocks

Let \(\Omega \subset \mathbb{R}^{spacedim}\) be the computational domain.

  • Primary unknown: displacement \(u:\Omega\to\mathbb{R}^{spacedim}\).

  • Optional secondary unknown: reduced Lagrange multipliers \(\lambda \in \mathbb{R}^{N_\lambda}\) attached to immersed inclusions (Inclusions<spacedim>).

In the code, the global solution is stored as a 2-block vector:

  • solution.block(0) \(\equiv u_h\),

  • solution.block(1) \(\equiv \lambda_h\) (may be size 0 if no inclusions).

The displacement is discretized with a continuous \(Q_p\) vector element: FESystem<spacedim>(FE_Q<spacedim>(fe_degree), spacedim).

Assembled matrices and right-hand sides

The following sparse matrices are assembled explicitly:

  • stiffness_matrix \(\equiv A\)

  • mass_matrix \(\equiv C\)

  • damping_matrix \(\equiv D\)

  • coupling_matrix \(\equiv B^T\) (maps \(\lambda\) to displacement dofs)

  • inclusion_matrix \(\equiv M\) (multiplier mass / scaling matrix)

Right-hand side blocks are:

  • system_rhs.block(0) \(\equiv f\) (volumetric + boundary loads and penalty terms),

  • system_rhs.block(1) \(\equiv g\) (inclusion data projected onto the multiplier basis),

  • system_rhs_f.block(0) \(\equiv f_f\) (computed in solve() as \(f_f=B^T\lambda\) in some branches; it is not assembled in assemble_elasticity_system()).

Constants and where they come from

Material constants are taken per-cell from MaterialProperties (see include/material_properties.h):

Model/solver constants are taken from ElasticityProblemParameters (see include/elasticity_problem_parameters.h):

A single length scale \(h\) used in the penalty term is computed as GridTools::minimal_cell_diameter(tria) in assemble_elasticity_system().

Bulk constitutive forms (what goes into A, C, D)

Define the symmetric gradient \(\varepsilon(u)=\frac12(\nabla u+\nabla u^T)\) and divergence \(\nabla\cdot u\).

Let \(\{\varphi_i\}\) be the displacement FE basis functions (vector-valued), and let \(\langle \cdot, \cdot \rangle\) denote the Euclidean inner product on vectors/tensors. The entries assembled in assemble_elasticity_system() are:

Mass matrix (always assembled):

\[ C_{ij} = \int_\Omega \rho\, \varphi_i \cdot \varphi_j\,dx. \]

Stiffness matrix contributions (assembled depending on flags):

Damping matrix contributions:

Body force right-hand side (always assembled):

\[ f_i \mathrel{+}= \int_\Omega \varphi_i \cdot f_{\text{rhs}}(x)\,dx, \]
where f_rhs comes from ElasticityProblemParameters::rhs.

Boundary conditions (as implemented)

Strong Dirichlet boundary conditions:

Normal-flux constraints:

Neumann boundary ids:

  • On boundary ids in ElasticityProblemParameters::neumann_ids, the code adds a scalar traction-like load computed as \(t_n(x)=\frac{1}{spacedim}\,(\texttt{Neumann\_bc}(x)\cdot n(x))\) and assembled as

    \[ f_i \mathrel{+}= -\int_{\Gamma_N} t_n\, \varphi_i\,ds. \]

Weak Dirichlet ids (penalty term):

  • For cells that have at least one boundary face with id in ElasticityProblemParameters::weak_dirichlet_ids, the code adds a penalty term using \(h = \) GridTools::minimal_cell_diameter(tria):

    \[ A_{ij} \mathrel{+}= \frac{\gamma_p}{h}\int_{\Omega_{\text{marked}}} \varphi_i\cdot\varphi_j\,dx,\qquad f_i \mathrel{+}= \frac{\gamma_p}{h}\,a_w\int_{\Omega_{\text{marked}}}\varphi_i\,dx. \]
    with \(\gamma_p=\) penalty_term and \(a_w=\) wave_ampltiude.

Inclusion coupling (B, M, g) and equation

types

The immersed inclusions define a reduced multiplier space with basis values (Fourier/harmonic modes) evaluated at quadrature points on each inclusion. In the code, assemble_coupling() builds:

  • Coupling matrix \(B^T\) (coupling_matrix, size \(N_u\times N_\lambda\)),

  • Inclusion matrix \(M\) (inclusion_matrix, size \(N_\lambda\times N_\lambda\), assembled diagonally),

  • Multiplier right-hand side \(g\) (system_rhs.block(1)).

Conceptually, the saddle-point branch (pressure_coupling == false) solves the discrete constrained system

\[\begin{split} \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} u \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ g \end{bmatrix}. \end{split}\]

The pressure_coupling == true branch does not solve the saddle-point problem. Instead it computes \(\lambda \approx M^{-1}g\) and then forms an additional load \(f_f = B^T\lambda\) (stored in system_rhs_f.block(0)).

Multiplier right-hand side \(g\) (as assembled in assemble_coupling()) is a normalized quadrature projection of the inclusion boundary data onto the selected Fourier coefficients. In particular, for each inclusion and each selected mode index \(j\), the code adds terms of the form

\[ g_j \mathrel{+}= \frac{1}{|\Gamma|}\sum_q \psi_j(q)\,d_j(q)\,ds_q, \]
with \(|\Gamma|=\) section_measure, basis value \(\psi_j(q)=\) inclusion_fe_values[j], and \(d_j(q)\) taken either from Inclusions::inclusions_rhs or from Inclusions::inclusions_data (see source/elasticity.cc, assemble_coupling()).

Time integration (Newmark predictor/corrector)

If ElasticityProblemParameters::initial_time != final_time, the solver uses Newmark time integration for the (semi-discrete) second-order system

\[ C\,\ddot u + D\,\dot u + A\,u = \big(f + f_f\big)\,\sin(t), \]
where \(f_f=B^T\lambda\) is present only in the pressure_coupling == true branch (as implemented in solve()).

Given \(u^n, v^n=\dot u^n, a^n=\ddot u^n\), the code computes: Predictor:

\[ u^{pred} = u^n + \Delta t\,v^n + \frac{\Delta t^2}{2}(1-2\beta_N)a^n,\qquad v^{pred} = v^n + \Delta t(1-\gamma_N)a^n. \]

Effective (acceleration) solve:

  • Linear elasticity:

    \[ \big(C + \beta_N\Delta t^2 A\big)a^{n+1} = (f+f_f)\sin(t^n) - A\,u^{pred}. \]

  • With damping (Rayleigh and/or Kelvin–Voigt):

    \[ \big(C + \beta_N\Delta t^2 A + \gamma_N\Delta t\,D\big)a^{n+1} = (f+f_f)\sin(t^n) - D\,v^{pred} - A\,u^{pred}. \]
    Corrector:
    \[ u^{n+1} = u^{pred} + \beta_N\Delta t^2 a^{n+1},\qquad v^{n+1} = v^{pred} + \gamma_N\Delta t\, a^{n+1}. \]

Predictor/corrector matrices

For matrix-based solvers/preconditioners that work on assembled matrices, the key “effective” (acceleration) matrix used in the Newmark solve is:

  • \(A_n = C + \beta_N\Delta t^2 A\) (undamped)

  • \(A_n = C + \beta_N\Delta t^2 A + \gamma_N\Delta t\,D\) (damped)

In the current implementation these combinations are built as LinearOperator sums in solve() (not as explicit sparse matrices).

Note

The bulk constitutive model is selected via ElasticityProblemParameters::elasticity_model (inferred from material parameters). Rayleigh damping is determined by the per-material parameters MaterialProperties::{rayleigh_alpha,rayleigh_beta}.

Note

In the stationary case (initial_time == final_time) with pressure_coupling == true, the current implementation computes \(f_f=B^T\lambda\) but then solves \(A u = f\) (i.e., \(f_f\) is not used in the stationary solve path).

Template Parameters:

  • dim – Topological dimension of the mesh (typically equal to spacedim).

  • spacedim – Embedding space dimension. The displacement has spacedim components.

Unnamed Group

DoFHandler<spacedim> dh

DoF metadata for displacement field.

std::vector<IndexSet> owned_dofs

Locally-owned DoF sets.

std::vector<IndexSet> relevant_dofs

Locally-relevant DoF sets.

Unnamed Group

AffineConstraints<double> constraints

Affine constraints in bulk and coupling spaces.

AffineConstraints<double> inclusion_constraints

Inclusion multiplier constraints.

AffineConstraints<double> mean_value_constraints

Mean-value constraints.

Unnamed Group

LA::MPI::SparseMatrix stiffness_matrix

Assembled sparse operators.

LA::MPI::SparseMatrix newmark_matrix

Effective matrix for Newmark solves.

LA::MPI::SparseMatrix mass_matrix

Bulk mass matrix.

LA::MPI::SparseMatrix coupling_matrix

Coupling matrix B^T.

LA::MPI::SparseMatrix damping_matrix

Viscous/Rayleigh damping matrix.

Unnamed Group

LA::MPI::PreconditionAMG prec_A

AMG preconditioners for displacement, Newmark, coupling, and mass blocks.

LA::MPI::PreconditionAMG prec_newmark

AMG for Newmark matrix.

LA::MPI::PreconditionAMG prec_C

AMG for coupling-related solves.

LA::MPI::PreconditionAMG prec_M

AMG for inclusion mass matrix.

Unnamed Group

LA::MPI::SparseMatrix inclusion_matrix

Coupling mass matrix and block vectors for state and rhs quantities.

LA::MPI::BlockVector solution

Current solution.

LA::MPI::BlockVector velocity

Current velocity.

LA::MPI::BlockVector acceleration

Current acceleration.

LA::MPI::BlockVector predictor

Newmark predictor.

LA::MPI::BlockVector corrector

Newmark corrector.

LA::MPI::BlockVector locally_relevant_solution

Ghosted solution.

LA::MPI::BlockVector force_rhs

Volumetric forcing rhs.

LA::MPI::BlockVector bc_rhs

Dirichlet-penalty rhs.

LA::MPI::BlockVector neumann_bc_rhs

Neumann rhs contribution.

LA::MPI::BlockVector system_rhs

Total rhs.

Unnamed Group

std::map<types::boundary_id, Tensor<1, spacedim>> forces

Boundary postprocessing accumulators and stress output vectors.

std::map<types::boundary_id, Tensor<1, spacedim>> average_displacements

Mean displacement per boundary id.

std::map<types::boundary_id, Tensor<1, spacedim>> average_normals

Mean normal per boundary id.

std::map<types::boundary_id, double> areas

Integrated boundary area/measure.

TrilinosWrappers::MPI::Vector sigma_n

Normal stress samples.

Public Functions

ElasticityProblem(const ElasticityProblemParameters<dim, spacedim> &par)

Build the elasticity problem from parsed parameters.

void make_grid()

Create or import computational mesh and manifolds.

void setup_fe()

Initialize finite element and quadrature objects.

void setup_dofs()

Distribute DoFs and initialize matrices/vectors.

void assemble_elasticity_system()

Assemble bulk operators (stiffness/mass/damping) and basic right-hand side.

void assemble_forcing_terms()

Assemble external forcing terms only.

void compute_system_rhs()

Rebuild the system right-hand side from component vectors.

void assemble_coupling()

Assemble immersed coupling operators and multiplier data.

void run()

Execute the run mode selected from parsed parameters.

set up, assemble and run the problem

void run_static()

Internal run implementations selected by time_mode.

void run_quasistatic()

Execute quasi-static time stepping.

void run_newmark()

Execute dynamic Newmark time stepping.

IndexSet assemble_coupling_sparsity(DynamicSparsityPattern &dsp)

Builds coupling sparsity, and returns locally relevant inclusion dofs.

void solve()

Solve the linear system for the active time mode.

void solve_static()

Internal solver implementations selected by time_mode.

void solve_quasistatic()

Solve one quasi-static step.

void solve_newmark()

Solve one dynamic Newmark step.

void refine_and_transfer()

Mark cells for refinement according to deal.ii routines fixed_number, fixed_fraction or global.

void refine_and_transfer_around_inclusions()

Mark cells around inclusions and transfer state vectors.

void execute_actual_refine_and_transfer()

Execute the pending refinement/coarsening and transfer on the cells marked by either “refine_and_transfer()” or “refine_and_transfer_around_inclusions()”.

std::string output_solution() const

Return output base filename for current cycle/time step.

void output_results() const

Write current solution and diagnostics to output files.

void print_parameters() const

Print selected runtime and model parameters.

void compute_internal_and_boundary_stress(bool) const

compute stresses on boundaries (2D and 3D) and internal (2D) this function makes use of boundary id, so make sure that the ifd starts from 0 and are sequential when importing a mesh, this is automatically taken care of for meshes generated with GridTools output is a txt file containing stresses at each cycle/time

void output_lambda() const

Write multiplier field values to disk.

std::string output_stresses() const

Write stress quantities and return produced filename stem.

Public Members

const ElasticityProblemParameters<dim, spacedim> &par

Parsed parameter object.

MPI_Comm mpi_communicator

MPI communicator used by distributed data structures.

ConditionalOStream pcout

Stream that prints only on rank zero.

mutable TimerOutput computing_timer

Timer collecting wall times by section.

parallel::distributed::Triangulation<spacedim> tria

Distributed bulk triangulation.

std::unique_ptr<FiniteElement<spacedim>> fe

Finite element used for displacement.

Inclusions<spacedim> inclusions

Inclusion geometry and reduced basis data.

std::unique_ptr<Quadrature<spacedim>> quadrature

Cell quadrature for bulk assembly.

std::unique_ptr<Quadrature<spacedim - 1>> face_quadrature_formula

Face quadrature for boundary terms.

std::vector<std::vector<BoundingBox<spacedim>>> global_bounding_boxes

Process-local coverings of the background mesh for particle insertion.

unsigned int cycle = 0

Current adaptive-refinement cycle index.

unsigned int time_step = 0

Current time-step index.

FEValuesExtractors::Vector displacement

Extractor for displacement components in FEValues.

double current_time = 0.0

Current physical time for transient simulations.