Contents
Linear Elasticity
Introduction
This tutorial provides a step-by-step guide for implementing a linear elasticity problem using the Finite Element Method (FEM) in the deal.II library. Linear elasticity governs the behavior of solid materials under small deformations, where stress is linearly proportional to strain. We will focus on setting up the weak form of the elasticity equations, assembling the finite element system, applying appropriate boundary conditions, and solving the resulting linear system. This tutorial assumes familiarity with basic FEM concepts and provides a practical foundation for more advanced developments such as nonlinear elasticity or dynamic simulations within the deal.II framework.
Theory
Linear elasticity is a branch of solid mechanics that deals with the behavior of solid materials under small deformations, where stress is proportional to strain. It is governed by Hooke's Law, which relates the stress and strain through material properties such as Young's modulus and Poisson's ratio. In this section, we will discuss the foundational equations of linear elasticity, including the stress-strain relationship, the equilibrium equations, and boundary conditions used in finite element analysis.
Stress-Strain Relationship:
$$ \sigma = \mathbf{C} \epsilon $$
Where:
- σ is the stress tensor
- C is the stiffness matrix (material property)
- ε is the strain tensor
Equilibrium Equation:
$$ \nabla \cdot \sigma + \mathbf{f} = 0 $$
These equations are essential for solving the linear elasticity problem using the Finite Element Method. The solution involves discretizing the domain into finite elements and solving the corresponding system of equations to obtain the displacement field.
Implementation
Include Files
In this section, we import all the necessary deal.II library components as well as standard C++ headers required to set up and solve the linear elasticity problem. These include modules for finite element definitions, mesh handling, quadrature formulas, linear solvers, and other utilities. Including these files ensures that we have access to all the classes and functions needed throughout the program.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/quadrature_point_data.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <fstream>
#include <iostream>
#include <random>
Namespace definition
We start by defining the namespace elas, which encapsulates all components related to the linear elasticity simulation. Namespaces in C++ help organize the code and prevent name conflicts, keeping related classes and functions grouped logically.
namespace elas
{
using namespace dealii;
Elasticity class
The Elasticity class is templated on the problem dimension and manages the overall workflow.
Public functions:
Elasticity(): Constructor to initialize the class.run(): Main driver function to execute the simulation.
Private functions:
create_mesh(): Generate the initial mesh.refine_mesh(): Refine the mesh after load increments.setup_constraints(): Apply boundary conditions and hanging node constraints.setup(): Set up system structures (matrices, vectors, DoFs).assemble(): Assemble the system matrix and RHS.solve_linear_problem(): Solve the elasticity system.output_results()andoutput_stress(): Export displacement and stress fields.reaction(): Compute reaction forces at the boundaries.
template <int dim>
class Elasticity
{
public:
Elasticity();
void run();
private:
// Mesh generation and refinement
void create_mesh();
void refine_mesh(const unsigned int load_step);
// Elasticity problem functions
void setup_constraints (const unsigned int load_step);
void setup(const unsigned int load_step);
void assemble();
void solve_linear_problem();
// Save the results - Disp, stress, reaction, energies
void output_results(const unsigned int load_step) const;
void output_stress(const unsigned int load_step);
void reaction(Tensor<1, dim> &reaction_stress, const types::boundary_id &boundary_id);
Parallel and MPI setup
We define these objects in order to store all the relevant information that allow the program to operate efficiently across distributed memory systems, specifically:
mpi_communicator: MPI communicator for all processes.n_mpi_processes: Number of active MPI processes.this_mpi_process: ID of the current process.
// Parallel setup and print on screen
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
parallel::distributed::Triangulation<dim> triangulation;
Core finite element objects
These objects define the mathematical and computational framework for the simulation
triangulation: Distributed mesh representation.fe: Finite element system for elasticity (displacement field).dof_handler: Manages degrees of freedom.locally_owned_dofs: Contains information of all the dofs of the mesh.locally_relevant_dofs: Stores only the relevant dofs for the current process.constraints: Handles boundary conditions, constraint equations and hanging nodes.
// Objects for elasticity
const FESystem<dim> fe; // FE System
DoFHandler<dim> dof_handler; // DoF Handler
IndexSet locally_owned_dofs; // IndexSet - Locally Owned
IndexSet locally_relevant_dofs; // IndexSet - Locally relevant
AffineConstraints<double> constraints; // Affine Constraint
const QGauss<dim> quadrature_formula; // Quadrature formula
System matrices and solution vectors:
system_matrix: The global elasticity stiffness matrix.system_rhs: The right-hand side vector.local_solution: The solution vector split among MPI processes. Each process owns only part of the global solution corresponding to its locally owned degrees of freedom.distributed_solution: The complete distributed solution vector where all processes share a consistent view across the entire computational domain.distributed_solution_old: A backup of the previous load step's distributed solution, primarily used during adaptive mesh refinement or solution continuation between steps.reaction_stress: A tensor storing computed reaction stresses, typically used for evaluating boundary tractions and reaction forces at constrained surfaces.
PETScWrappers::MPI::SparseMatrix system_matrix; // Elastcitiy Matrix
PETScWrappers::MPI::Vector system_rhs; // RHS
PETScWrappers::MPI::Vector local_solution; // MPI Split solution
PETScWrappers::MPI::Vector distributed_solution; // Full Solution
PETScWrappers::MPI::Vector distributed_solution_old; // This is for the refinement
Tensor<1, dim> reaction_stress;
Material Properties and Simulation Parameters
This section defines the physical and numerical setup for the elasticity simulation. The domain dimensions are specified by L, H, and W, which represent the length, height, and width, respectively. The initial subdivisions of the domain along the x, y, and z axes are controlled by the nx, ny, and nz variables.
The material properties are defined using E for Young's modulus and nu for Poisson's ratio. These values are essential for simulating the material's response to stress. The simulation load is managed with the load_increment, which dictates the amount of load applied at each step, and the applied_load, which tracks the cumulative load over time. The num_load_steps variable defines the total number of load increments to be applied during the simulation.
Two auxiliary tools are utilized in this setup: pcout, which handles parallel printing, allowing output to be displayed during parallel computations, and computing_timer, which measures and reports the computation time for each simulation step.
This configuration forms the foundation for the elasticity simulation, specifying the material properties, domain dimensions, and numerical controls required to run the simulation efficiently and accurately.
// Domain dimensions
double L = 1;
double H = 0.05;
double W = 0.05;
const unsigned int nx = 25;
const unsigned int ny = 2;
const unsigned int nz = 2;
// Material properties
double E = 100;
double nu = 0.3;
double load_increment = 0.01; // Load increment applied
double applied_load = 0; // Current applied_load += load_increment
double num_load_steps = 5; // Total load increment
ConditionalOStream pcout;
TimerOutput computing_timer;
}; // End of Elast class
Constructor
The constructor initializes all essential components required for the elasticity problem, including the core finite element objects and the MPI communicator-related objects. Each object associated with the communicator is initialized in the same order in which it is declared within the class.
// Constructor
template <int dim>
Elasticity<dim>::Elasticity()
: mpi_communicator(MPI_COMM_WORLD)
, n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
, this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
, triangulation (mpi_communicator)
, fe (FE_Q<dim> (1), dim)
, dof_handler (triangulation)
, quadrature_formula (fe.degree + 1)
, pcout (std::cout, this_mpi_process == 0)
, computing_timer (mpi_communicator, pcout,
TimerOutput::never,
TimerOutput::wall_times)
{}
Evaluation of Elastic Constants:
These two inline functions compute the Lamé constants, lambda and mu, based on the material's Young's modulus (E) and Poisson's ratio (nu). These constants are fundamental for solving elasticity problems in solid mechanics.
lambda(E, nu):$$\lambda = \frac{E \, \nu}{(1 + \nu) \, (1 - 2\nu)}$$
-
mu(E, nu):$$ \mu = \frac{E}{(2 (1 + \nu))} $$
// Evaluation of the elastic constant
inline double lambda (const float E, const float nu)
{
return (E * nu) / ((1 + nu) * (1 - 2 * nu));
}
inline double mu (const float E, const float nu)
{
return E / (2 * (1 + nu));
}
Create mesh
This function generates the computational mesh for the elasticity problem, providing flexibility for both 2D and 3D cases.
It starts by defining a coarse subdivision based on the given domain dimensions, which is stored in a standard C++ vector.
The next step is to define the bottom-left point p1 and the top-right point p2 of the rectangular domain as Point objects.
To ensure the solution is dimension-independent, a check is performed to determine the problem's dimensionality, and 2D or 3D points are created accordingly.
Finally, the mesh is created using the GridGenerator::subdivided_hyper_rectangle function.
template <int dim>
void Elasticity<dim>::create_mesh()
{
// Define the initial coarse subdivision
const std::vector<unsigned int> repetitions = {nx,ny,nz};
// Create Mesh
Point<dim> p1, p2;
if constexpr (dim == 2)
{
p1 = Point<dim>(0, 0);
p2 = Point<dim>(L, H);
}
else if constexpr (dim == 3)
{
p1 = Point<dim>(0, 0, 0);
p2 = Point<dim>(L, H, W);
}
// create coarse mesh
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, p1, p2);
Setting Boundary Face IDs:
The code iterates over all active cells in the triangulation object, and for each cell, it iterates over its faces using cell->face_iterators(), it checks if the face is at the boundary with face->at_boundary() and lastly, if the face is at the boundary, the center of the face is calculated using face->center().
Next, the code checks the location of the face's center to determine its position relative to the boundaries of the domain. It first checks if the face is located at the left, right, bottom, or top, assigning the IDs 0, 1, 2, and 3, respectively.
To remain dimension independent, the code uses if constexpr (dim == 3) to ensure that the check for the third coordinate is performed only in 3D cases. If the face is located at z = 0, it assigns a boundary ID of 4, and if the face is at z = W, it assigns a boundary ID of 5.
// Set Face IDs
for (const auto &cell : triangulation.active_cell_iterators ())
{
for (const auto &face : cell->face_iterators ())
{
if (face->at_boundary ())
{
const auto center = face->center ();
if (std::fabs (center(0) - 0) < 1e-12)
face->set_boundary_id (0);
else if (std::fabs (center(0) - L) < 1e-12)
face->set_boundary_id (1);
else if (std::fabs (center(1) - 0) < 1e-12)
face->set_boundary_id (2);
else if (std::fabs (center(1) - H) < 1e-12)
face->set_boundary_id (3);
// Check if we are in 3D before accessing center(2)
if constexpr (dim == 3)
{
if (std::fabs(center(2) - 0) < 1e-12)
face->set_boundary_id(4);
else if (std::fabs(center(2) - W) < 1e-12)
face->set_boundary_id(5);
}
}
}
}
Initialize DoF Handler and Solution Setup:
This section of the code initializes the degrees of freedom (DoFs) and sets up the solution for the elasticity problem in a parallel execution environment. The first operation performed is the distribution of DoFs across the mesh using dof_handler.distribute_dofs(fe). This method ensures that the degrees of freedom are allocated according to the finite element method configuration, which is crucial for efficient parallel execution, where each process is responsible for a part of the problem.
Next, the code retrieves the locally owned degrees of freedom using locally_owned_dofs = dof_handler.locally_owned_dofs(). This step ensures that each MPI process knows which part of the problem it is responsible for. Additionally, it extracts the locally relevant DoFs, which may include DoFs owned by neighboring processes. This is done using DoFTools::extract_locally_relevant_dofs(dof_handler), enabling the necessary inter-process communication for solving the problem.
To maintain continuity between load steps, especially when using adaptive mesh refinement techniques, the solution vector from the previous load step is reinitialized. The vector is called distributed_solution_old, and it is reinitialized using distributed_solution_old.reinit(locally_owned_dofs, mpi_communicator). This ensures that the previous solution is preserved for refinement purposes.
The code also prints some key information regarding the triangulation and the distribution of the mesh. Using the pcout stream, the number of global levels in the triangulation is displayed with the command. This gives an overview of the mesh refinement hierarchy. Similarly, it prints the number of locally owned cells, which informs the user about the number of cells assigned to the current MPI process. The total number of cells across the global domain is also printed.
// Initialize dof handler and extact Owned/Relevant DoFs
dof_handler.distribute_dofs (fe);
locally_owned_dofs = dof_handler.locally_owned_dofs ();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs (dof_handler);
// Initialize the old solution for refinement purposes
distributed_solution_old.reinit (locally_owned_dofs, mpi_communicator);
pcout << "No. of levels in triangulation: "
<< triangulation.n_global_levels () << std::endl;
pcout << " Number of locally owned cells on the process: "
<< triangulation.n_locally_owned_active_cells () << std::endl;
pcout << "Number of global cells:" << triangulation.n_global_active_cells () << std::endl;
pcout << " Total Number of globally active cells: "
<< triangulation.n_global_active_cells () << std::endl
<< " Number of degrees of freedom for elasticity: "
<< dof_handler.n_dofs () << std::endl;
}
Mesh Refinement
This function prepares the computational mesh for adaptive refinement based on an error estimation from the current solution.
It is typically called at the end of each load step to improve the solution's accuracy by refining the mesh where needed.
The first step in the mesh refinement process is the creation of a vector, Vector<float> estimated_error_per_cell,
which stores the estimated error for each locally owned active cell.
The vector is initialized to the number of cells owned by the current MPI process, ensuring that each process has its own set of error estimates.
The next part of the process involves the use of the KellyErrorEstimator<dim>::estimate function.
This function estimates the error for each cell in the mesh using the Kelly error estimator,
which is a standard technique for error estimation in finite element methods.
The estimator requires several inputs: the dof_handler (which contains the degrees of freedom),
a quadrature rule QGauss<dim-1>(fe.degree + 1) that specifies the degree of quadrature for integration,
a map of boundary IDs, a function representing boundary conditions, and the current solution local_solution.
The estimated error for each cell is then stored in the estimated_error_per_cell vector.
Finally, the code initializes a SolutionTransfer object.
The parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector> is responsible for transferring
the solution from the old mesh to the refined mesh.
This ensures that the solution is correctly mapped after the adaptive refinement, especially in a parallel execution environment.
The SolutionTransfer object, named soltrans, is initialized with the dof_handler
to enable the correct transfer of solution data.
template<int dim>
void Elasticity<dim>::refine_mesh(const unsigned int load_step)
{
Vector<float> estimated_error_per_cell (triangulation.n_locally_owned_active_cells ());
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1> (fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
local_solution, estimated_error_per_cell);
// Initialize SolutionTransfer object
parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector>
soltrans (dof_handler);
Mesh Refinement and Solution Transfer:
This section of the code performs adaptive mesh refinement based on error estimation, followed by transferring the solution across the updated mesh.
The parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction function refines the top 30% of cells with the highest error and coarsens the bottom 1% to improve accuracy.
If the mesh is refined to level 4 or higher, further refinement is prevented at level 3 to avoid unnecessary computational overhead. This is controlled by checking the number of global levels with if (triangulation.n_global_levels() >= 4).
The triangulation is then prepared for coarsening and refinement using triangulation.prepare_coarsening_and_refinement(), and the SolutionTransfer object is initialized with soltrans.prepare_for_coarsening_and_refinement(local_solution) to ensure the current solution is correctly mapped onto the refined mesh.
After the mesh is refined, the degrees of freedom are redistributed using dof_handler.distribute_dofs(fe), ensuring the solution is consistent with the new mesh configuration.
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (
triangulation, estimated_error_per_cell, 0.30, 0.01);
if (triangulation.n_global_levels () >= 4)
{
for (const auto &cell : triangulation.active_cell_iterators_on_level (3))
if (cell->is_locally_owned ())
cell->clear_refine_flag ();
}
// Prepare the triangulation
triangulation.prepare_coarsening_and_refinement ();
// Prepare the SolutionTransfer object for coarsening and refinement
// and give the solution vector that we intend to interpolate later,
soltrans.prepare_for_coarsening_and_refinement (local_solution);
// Do the refinement
triangulation.execute_coarsening_and_refinement ();
// redistribute dofs,
dof_handler.distribute_dofs (fe);
Reconstruction of Degrees of Freedom and Solution Interpolation:
This part of the code manages the redistribution of degrees of freedom (DoFs) and interpolates the old solution onto the newly refined mesh. Constraints are then reapplied to ensure the solution remains valid and consistent with the updated mesh structure.
The locally_owned_dofs = dof_handler.locally_owned_dofs() updates the degrees of freedom owned by the current MPI process after mesh refinement. Similarly, locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler) extracts the locally relevant degrees of freedom, including those needed for MPI communication.
The distributed_solution_old.reinit(locally_owned_dofs, mpi_communicator) reinitializes the distributed solution vector to align with the new DoF layout, followed by interpolation of the old solution onto the refined mesh with soltrans.interpolate(distributed_solution_old).
Next, the setup_constraints(load_step) function re-applies the boundary conditions and other constraints, and constraints.distribute(distributed_solution_old) ensures the constraints are applied to the interpolated solution. Finally, the new solution is copied into local_solution to maintain continuity in the simulation.
// Recreate locally_owned_dofs and locally_relevant_dofs index sets
locally_owned_dofs = dof_handler.locally_owned_dofs ();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs (dof_handler);
distributed_solution_old.reinit (locally_owned_dofs, mpi_communicator);
soltrans.interpolate (distributed_solution_old);
// Apply constraints on the interpolated solution
setup_constraints(load_step);
constraints.distribute (distributed_solution_old);
// Copy distributed_solution_old to local_solution
local_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
local_solution = distributed_solution_old;
}
Setup Constraints
This function configures the boundary constraints needed to enforce boundary conditions on the finite element solution.
It starts by clearing and reinitializing the constraint system based on the currently relevant degrees of freedom.
Hanging node constraints, which ensure continuity across mesh refinements, are automatically handled through DoFTools::make_hanging_node_constraints.
Specific displacement conditions are applied next.
Displacement along the x-axis is extracted using u_x and its corresponding mask.
The displacement values on the left and right boundaries are computed based on the current load step,
introducing symmetrical loading: negative on the left side and positive on the right side.
Similarly, displacement along the y-axis is extracted using u_y, with a fixed (zero) displacement applied
on the relevant boundary to prevent vertical motion.
If the simulation is in three dimensions, displacement along the z-axis (u_z) is also constrained to zero
on a specified boundary to avoid out-of-plane motion.
Finally, the constraint system is closed to finalize all imposed conditions before solving the system.
template<int dim>
void Elasticity<dim>::setup_constraints (const unsigned int load_step)
{
constraints.clear ();
constraints.reinit (locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
const FEValuesExtractors::Scalar u_x (0);
const FEValuesExtractors::Scalar u_y (1);
const ComponentMask u_x_mask = fe.component_mask (u_x);
const ComponentMask u_y_mask = fe.component_mask (u_y);
applied_load = load_increment*load_step;
const double u_x_values_right = 0.5*applied_load;
const double u_x_values_left = -0.5*applied_load;
const double u_fix = 0.0;
VectorTools::interpolate_boundary_values (dof_handler, 0,
Functions::ConstantFunction<dim> (u_x_values_left, dim),
constraints,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler, 1,
Functions::ConstantFunction<dim> (u_x_values_right, dim),
constraints,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler, 2,
Functions::ConstantFunction<dim> (u_fix, dim),
constraints,
u_y_mask);
if constexpr (dim == 3)
{
const FEValuesExtractors::Scalar u_z (2);
const ComponentMask u_z_mask = fe.component_mask (u_z);
VectorTools::interpolate_boundary_values (dof_handler, 4,
Functions::ConstantFunction<dim> (u_fix, dim),
constraints,
u_z_mask);
}
constraints.close ();
}
System setup
This function initializes all essential structures required for solving the
finite element system at the current load step.
It starts by creating and extracting the degrees of freedom that are locally
owned and those that are locally relevant for communication across MPI processes.
Solution vectors such as local_solution, system_rhs,
and distributed_solution are reinitialized accordingly to match
the new DoF layout.
Boundary conditions and hanging node constraints are enforced through a call
to setup_constraints(load_step).
A dynamic sparsity pattern is then generated,
defining where non-zero entries are expected in the system matrix.
This sparsity pattern is distributed across MPI processes to ensure
efficient parallel assembly.
Finally, system_matrix is reinitialized
based on the new sparsity pattern, completing the setup required before
assembling and solving the elasticity problem for the given load step.
template<int dim>
void Elasticity<dim>::setup(const unsigned int load_step)
{
TimerOutput::Scope ts(computing_timer, "setup_system");
// Initialize the dofs contained locally as well as those of the near cells
locally_owned_dofs = dof_handler.locally_owned_dofs ();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs (dof_handler);
local_solution.reinit (locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
// Initialize the contraints objects for BC and Hanging Nodes during refinement
system_rhs.reinit (locally_owned_dofs, mpi_communicator);
distributed_solution.reinit (locally_owned_dofs, mpi_communicator);
// Initialize hanging nodes and constrains
setup_constraints (load_step);
// Sparsity pattern initialization (where matrix has nonzero entries)
DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern (dof_handler,
dsp,
constraints,
false);
SparsityTools::distribute_sparsity_pattern (dsp,
dof_handler.locally_owned_dofs (),
mpi_communicator,
locally_relevant_dofs);
// Elasticity Matrix
system_matrix.reinit (locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
Assembly
This function performs the assembly process of the system matrix and rhs.
First it clears the previous system_rhs and system_matrix.
It then initializes the FEValues object,
which manages shape function evaluations, gradients, quadrature points,
and Jacobian-weighted measures for each cell during integration.
The number of degrees of freedom per cell and the number of
quadrature points are determined.
Local objects, including the element cell_matrix and cell_rhs,
are created and sized appropriately. Additionally, a vector
local_dof_indices is initialized to map local degrees of freedom
to the global system, setting the stage for the element-wise assembly process.
template<int dim>
void Elasticity<dim>::assemble()
{
TimerOutput::Scope ts (computing_timer, "assembly");
// Clean the current matrices
system_rhs = 0;
system_matrix = 0;
// Initialize the FEValues Objects
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
// Number of DoFs and Gauss points per cell
const unsigned int dofs_per_cell = fe.n_dofs_per_cell ();
// Local Matrices and Vector
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
// Create vector to store the local dof indices
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
The function loops over all locally owned active cells.
For each cell, the local stiffness matrix cell_matrix
and local right-hand side vector cell_rhs are reset to zero.
The finite element information is updated for the current cell using
fe_values.reinit(cell), preparing the shape functions and
their gradients for integration.
At each Gauss quadrature point, the code computes contributions to the local
matrix by assembling the linear elasticity terms based on Lame parameters
lambda and mu. It also calculates contributions
to the local right-hand side vector, though in this case the body force is zero.
To apply boundary forces, you need to loop over the faces of the cells and check if each face belongs to a boundary.
If the face is part of the boundary, you can then add the corresponding contribution
to the system. This requires to setup and use a FEFaceValues object in addition
to the FEValues already in place, which evaluates
the shape functions and gradients at the domain boundaries.
By using FEFaceValues, you can compute the necessary values
(such as shape functions, gradients, and normal vectors) on the boundary faces,
allowing you to correctly apply the boundary forces to the right-hand side vector (system_rhs) of the linear system.
Afterward, the local degrees of freedom are mapped to their global indices, and the local contributions are inserted into the global system matrix and right-hand side using constraints.distribute_local_to_global, which also enforces boundary and hanging node conditions.
Finally, after all cells have been processed, the system matrix and right-hand side vector are compressed to finalize assembly across all MPI processes, ensuring a consistent global system ready for solving.
// loop over the cells
for (const auto &cell : dof_handler.active_cell_iterators ())
{
if (cell->is_locally_owned ())
{
// Reset the local matrix/rhs of the cell
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
// loop over the gauss point of the cell
for (const unsigned int q_point : fe_values.quadrature_point_indices ())
{
for (const unsigned int i : fe_values.dof_indices ())
{
const unsigned int component_i = fe.system_to_component_index (i).first;
for (const unsigned int j : fe_values.dof_indices ())
{
const unsigned int component_j = fe.system_to_component_index (j).first;
cell_matrix (i, j) += ((fe_values.shape_grad (i, q_point)[component_i]
* lambda (E, nu) *
fe_values.shape_grad (j, q_point)[component_j])
+
(fe_values.shape_grad (i, q_point)[component_j]
* mu (E, nu) *
fe_values.shape_grad (j, q_point)[component_i])
+
((component_i == component_j) ?
(fe_values.shape_grad (i, q_point) * mu (E, nu)
* fe_values.shape_grad (j, q_point))
: 0)
) * fe_values.JxW (q_point); //
}
// Cell rhs set to zero for the current example
// Modify this entry to account for body forces
cell_rhs (i) += fe_values.shape_value (i, q_point) * 0
* fe_values.JxW (q_point);
// To apply boundary forces you need to loop over the exterior boundaries
// to populate the corresponding entries of the rhs
}
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
system_matrix.compress (VectorOperation::add);
system_rhs.compress (VectorOperation::add);
}
Solve system
The function begins by measuring the computation time for solving the linear system. A Conjugate Gradient (CG) solver is initialized with a maximum iteration count and a stopping tolerance proportional to the norm of the right-hand side vector. The system is solved using the CG method with a block Jacobi preconditioner built from the system matrix to accelerate convergence.
After solving, the number of iterations taken to converge is printed using the parallel output stream. Finally, the solution vector is adjusted by applying the constraints, ensuring the solution respects boundary conditions and hanging node continuity requirements.
template<int dim>
void Elasticity<dim>::solve_linear_problem()
{
TimerOutput::Scope ts (computing_timer, "solve_linear_system");
// CG
SolverControl solver_control (100000, 1e-12* system_rhs.l2_norm());
SolverCG<PETScWrappers::MPI::Vector> solver (solver_control);
PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
solver.solve (system_matrix,
distributed_solution,
system_rhs,
preconditioner);
pcout << " Solved in " << solver_control.last_step () << " iterations." << std::endl;
constraints.distribute (distributed_solution);
}
Output results
This function handles writing the displacement solution and subdomain information after each load step. Displacement components are labeled appropriately and interpreted as vector fields for visualization. Each processor’s subdomain is also recorded to distinguish parallel ownership of cells.
The data is organized into patches and written into VTU format files, with a PVTU master record created for parallel output. On the root MPI process, a PVD file is updated to index all the time steps, enabling easy loading of the full simulation history in visualization tools like ParaView.
template <int dim>
void Elasticity<dim>::output_results (const unsigned int load_step) const
{
std::vector<std::string> displacement_names (dim, "displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
displacement_component_interpretation (dim,
DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> data_out_phasefield;
data_out_phasefield.add_data_vector (dof_handler,
local_solution,
displacement_names,
displacement_component_interpretation);
Vector<double> subdomain (triangulation.n_active_cells ());
for (unsigned int i = 0; i < subdomain.size (); ++i)
{
subdomain (i) = triangulation.locally_owned_subdomain ();
}
data_out_phasefield.add_data_vector (subdomain, "subdomain");
data_out_phasefield.build_patches ();
const std::string pvtu_filename = data_out_phasefield.write_vtu_with_pvtu_record ("./",
"solution",
load_step, mpi_communicator,
2,
0);
if (this_mpi_process == 0)
{
static std::vector<std::pair<double, std::string>> times_and_names;
times_and_names.push_back(std::pair<double, std::string>(load_step, pvtu_filename));
std::ofstream pvd_output("./solution_u.pvd");
DataOutBase::write_pvd_record(pvd_output, times_and_names);
}
}
Output stresses
This function calculates the stress tensor at the quadrature points
by post-processing the displacement solution.
The stress components are averaged per cell and stored in vectors.
Separate stress components such as s_xx, s_xy, and s_yy
(and s_zz, s_xz, s_yz in 3D)
are created and attached to the mesh for visualization.
Each stress component is exported into VTU files using DataOut.
A parallel VTU (PVTU) record is generated for distributed visualization,
and a PVD file is updated on the root process to maintain the simulation timeline.
template <int dim>
void Elasticity<dim>::output_stress(const unsigned int load_step)
{
// Vettori per cella in Vector<float>
Vector<float> s_xx(triangulation.n_active_cells());
Vector<float> s_xy(triangulation.n_active_cells());
Vector<float> s_yy(triangulation.n_active_cells());
Vector<float> s_zz(triangulation.n_active_cells());
Vector<float> s_xz(triangulation.n_active_cells());
Vector<float> s_yz(triangulation.n_active_cells());
FEValues<dim> fe_values(fe, quadrature_formula,
update_gradients | update_JxW_values);
const unsigned int n_q_points = quadrature_formula.size();
const FEValuesExtractors::Vector displacements(0);
std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
fe_values[displacements].get_function_symmetric_gradients(local_solution,
strain_tensor);
double c_s_xx = 0, c_s_xy = 0, c_s_yy = 0;
double c_s_zz = 0, c_s_xz = 0, c_s_yz = 0;
for (unsigned int q = 0; q < n_q_points; ++q)
{
const auto &eps_u = strain_tensor[q];
const double tr = trace(eps_u);
const auto stress = 2.0 * mu(E, nu) * eps_u + lambda(E, nu)
* tr * unit_symmetric_tensor<dim>();
c_s_xx += stress[0][0];
c_s_xy += stress[0][1];
c_s_yy += stress[1][1];
if constexpr (dim == 3)
{
c_s_zz += stress[2][2];
c_s_xz += stress[0][2];
c_s_yz += stress[1][2];
}
}
const unsigned int cell_index = cell->active_cell_index();
s_xx[cell_index] = c_s_xx / n_q_points;
s_xy[cell_index] = c_s_xy / n_q_points;
s_yy[cell_index] = c_s_yy / n_q_points;
if constexpr (dim == 3)
{
s_zz[cell_index] = c_s_zz / n_q_points;
s_xz[cell_index] = c_s_xz / n_q_points;
s_yz[cell_index] = c_s_yz / n_q_points;
}
}
}
DataOut<dim> data_out;
data_out.attach_triangulation(triangulation);
data_out.add_data_vector(s_xx, "s_xx", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_xy, "s_xy", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_yy, "s_yy", DataOut<dim>::type_cell_data);
if constexpr (dim == 3)
{
data_out.add_data_vector(s_zz, "s_zz", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_xz, "s_xz", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_yz, "s_yz", DataOut<dim>::type_cell_data);
}
data_out.build_patches();
const std::string pvtu_filename = data_out.write_vtu_with_pvtu_record(
"./", "stresses", load_step, mpi_communicator, 4);
if (this_mpi_process == 0)
{
static std::vector<std::pair<double, std::string>> times_and_names;
times_and_names.emplace_back(load_step, pvtu_filename);
std::ofstream pvd_output("./stresses.pvd");
DataOutBase::write_pvd_record(pvd_output, times_and_names);
}
}
Compute reactions
This function computes the total reaction force on a specified boundary by integrating the stress tensor over the boundary faces. The stress is derived from the strain, which is computed from the displacement solution at quadrature points on the faces.
Only faces matching the given boundary_id are considered. Local contributions are accumulated per MPI rank and summed globally. The resulting reaction vector components are printed using parallel output utilities.
// Reaction
template <int dim>
void Elasticity<dim>::reaction(Tensor<1, dim> &reaction_stress,
const types::boundary_id &boundary_id)
{
// Face quadrature formula
QGauss<dim - 1> face_quadrature_formula(5);
// Set up FE face values for vector-valued and scalar-valued solution
FEFaceValues<dim> fe_face_values(fe, face_quadrature_formula,
UpdateFlags(update_values |
update_gradients |
update_quadrature_points |
update_normal_vectors |
update_JxW_values));
const unsigned int n_face_q_points = face_quadrature_formula.size();
// Extract displacements for gradient calculation
const FEValuesExtractors::Vector displacements(0);
std::vector<SymmetricTensor<2, dim>> strain_tensor(n_face_q_points);
// Parallel DoFHandler iterator
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
Tensor<1, dim> reaction_mpi;
// Loop over all cells in the triangulation
for (; cell != endc; ++cell)
{
if (cell->is_locally_owned())
{
// Loop over all faces of the current cell
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
{
// Check if the current face belongs to the specified boundary
if (cell->face(face)->boundary_id() == boundary_id)
{
// Reinitialize FE face values for the current face
fe_face_values.reinit(cell, face);
fe_face_values[displacements].get_function_symmetric_gradients(local_solution, strain_tensor);
// Loop over quadrature points on the face
for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
{
const auto &eps_u = strain_tensor[q_point];
const Tensor<1, dim> &N = fe_face_values.normal_vector(q_point);
const double JxW_f = fe_face_values.JxW(q_point);
// Compute the stress tensor using the strain
const double tr = trace(eps_u);
const auto stress = 2.0 * mu(E, nu) * eps_u + lambda(E, nu)
* tr * unit_symmetric_tensor<dim>();
// Accumulate the contribution to the reaction stress
reaction_mpi += stress * N * JxW_f;
}
}
}
}
}
// Sum the reaction stresses from all MPI ranks and print it
reaction_stress = Utilities::MPI::sum(reaction_mpi, mpi_communicator);
for (unsigned int i = 0; i < dim; ++i)
{
pcout << " R_" << (char)('x' + i) << ": " << reaction_stress[i] << std::endl;
}
}
Run
This is the main driver function. It starts by creating the mesh and then iterates over each load step. For each step, it sets up the system, assembles the linear system, solves it, updates the solution, refines the mesh, and outputs displacement, stress, and reaction results. Performance information is printed at each load step, and the overall runtime is recorded at the end.
template <int dim>
void Elasticity<dim>::run()
{
Timer timer;
timer.start ();
pcout << "Running on " << Utilities::MPI::n_mpi_processes (mpi_communicator)
<< " MPI rank(s)..." << std::endl;
create_mesh();
// Loop over load steps
for (unsigned int load_step = 0; load_step <= num_load_steps; load_step++)
{
pcout << " --- Load increment number : " << load_step << std::endl;
setup (load_step);
assemble ();
solve_linear_problem ();
// Update ghost values after the solution
local_solution = distributed_solution;
local_solution.update_ghost_values();
// Prepare for refinement
distributed_solution_old = local_solution;
// Refine
refine_mesh(load_step);
pcout << " Total New Number of globally active cells: "
<< triangulation.n_global_active_cells () << std::endl
<< " New Number of degrees of freedom for elasticity: "
<< dof_handler.n_dofs () << std::endl;
// Save results
output_results (load_step);
output_stress (load_step);
reaction(reaction_stress, 1);
// Print the table of cpmputing times
computing_timer.print_summary ();
computing_timer.reset ();
pcout << std::endl;
}
timer.stop ();
pcout << "Total run time: " << timer.wall_time () << " seconds." << std::endl;
}
} // End of namespace Elasticity
Main
The entry point of the program. It initializes the MPI environment and creates
an instance of the Elasticity<dim> class (with <dim> to be set equal to 2 for the 2D problem or 3 for 3D one).
It then calls the run() method to execute the simulation.
Exceptions are caught and printed clearly to standard error to ensure that any
runtime issues are reported properly before the program aborts.
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace elas;
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
Elasticity<2> elasticity; // For the 2D problem
// Elasticity<3> elasticity; // For the 3D problem
elasticity.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what ()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
return 1;
}
return 0;
}
Results
Sezione dei risultati in cui andare a mettere:
- Plot degli spostamenti ad ogni passo temporale
- Plot dei dof per ciclo di raffittimento
- Plot delle tensioni
- Grafico carico vs spostamento
Complete code
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/quadrature_point_data.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <fstream>
#include <iostream>
#include <random>
namespace elas
{
using namespace dealii;
template <int dim>
class Elasticity
{
public:
Elasticity();
void run();
private:
// Mesh generation and refinement
void create_mesh();
void refine_mesh(const unsigned int load_step);
// Elasticity problem functions
void setup_constraints (const unsigned int load_step);
void setup(const unsigned int load_step);
void assemble();
void solve_linear_problem();
// Save the results - Disp, stress, reaction
void output_results(const unsigned int load_step) const;
void output_stress(const unsigned int load_step);
void reaction(Tensor<1, dim> &reaction_stress, const types::boundary_id &boundary_id);
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
parallel::distributed::Triangulation<dim> triangulation;
// Elasticity objects
const FESystem<dim> fe; // FE System
DoFHandler<dim> dof_handler; // DoF Handler
IndexSet locally_owned_dofs; // IndexSet - Locally Owned
IndexSet locally_relevant_dofs; // IndexSet - Locally relevant
AffineConstraints<double> constraints; // Affine Constraint
PETScWrappers::MPI::SparseMatrix system_matrix; // Elastcitiy Matrix
PETScWrappers::MPI::Vector system_rhs; // RHS
PETScWrappers::MPI::Vector local_solution; // MPI Split solution
PETScWrappers::MPI::Vector distributed_solution; // Full Solution
PETScWrappers::MPI::Vector distributed_solution_old; // This is for the refinement
const QGauss<dim> quadrature_formula; // Quadrature formula
Tensor<1, dim> reaction_stress;
// Domain dimensions
double L = 1;
double H = 0.05;
double W = 0.05;
const unsigned int nx = 25;
const unsigned int ny = 2;
const unsigned int nz = 2;
// Material properties
double E = 100;
double nu = 0.3;
double load_increment = 0.01; // Load increment applied
double applied_load = 0; // Current applied_load += load_increment
double num_load_steps = 5; // Total load increment
ConditionalOStream pcout;
TimerOutput computing_timer;
}; // End of Elast class
// Constructor
template <int dim>
Elasticity<dim>::Elasticity()
: mpi_communicator(MPI_COMM_WORLD)
, n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
, this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
, triangulation (mpi_communicator)
, fe (FE_Q<dim> (1), dim)
, dof_handler (triangulation)
, quadrature_formula (fe.degree + 1)
, pcout (std::cout, this_mpi_process == 0)
, computing_timer (mpi_communicator, pcout, TimerOutput::never, TimerOutput::wall_times)
{}
// Evaluation of the elastic constant
inline double lambda (const float E, const float nu)
{
return (E * nu) / ((1 + nu) * (1 - 2 * nu));
}
inline double mu (const float E, const float nu)
{
return E / (2 * (1 + nu));
}
// -----------------------------
// ----------- MESH ------------
// -----------------------------
template <int dim>
void Elasticity<dim>::create_mesh()
{
// Define the initial coarse subdivision
const std::vector<unsigned int> repetitions = {nx,ny,nz};
// Create Mesh
Point<dim> p1, p2;
if constexpr (dim == 2)
{
p1 = Point<dim>(0, 0);
p2 = Point<dim>(L, H);
}
else if constexpr (dim == 3)
{
p1 = Point<dim>(0, 0, 0);
p2 = Point<dim>(L, H, W);
}
// create coarse mesh
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, p1, p2);
// Set Fce IDs
for (const auto &cell : triangulation.active_cell_iterators ())
{
for (const auto &face : cell->face_iterators ())
{
if (face->at_boundary ())
{
const auto center = face->center ();
if (std::fabs (center(0) - 0) < 1e-12)
face->set_boundary_id (0);
else if (std::fabs (center(0) - L) < 1e-12)
face->set_boundary_id (1);
else if (std::fabs (center(1) - 0) < 1e-12)
face->set_boundary_id (2);
else if (std::fabs (center(1) - H) < 1e-12)
face->set_boundary_id (3);
// Check if we are in 3D before accessing center(2)
if constexpr (dim == 3)
{
if (std::fabs(center(2) - 0) < 1e-12)
face->set_boundary_id(4);
else if (std::fabs(center(2) - W) < 1e-12)
face->set_boundary_id(5);
}
}
}
}
// Initialize dof handler and extact Owned/Relevant DoFs
dof_handler.distribute_dofs (fe);
locally_owned_dofs = dof_handler.locally_owned_dofs ();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs (dof_handler);
// Initialize the old solution for refinement purposes
distributed_solution_old.reinit (locally_owned_dofs, mpi_communicator);
pcout << "No. of levels in triangulation: "
<< triangulation.n_global_levels () << std::endl;
pcout << " Number of locally owned cells on the process: "
<< triangulation.n_locally_owned_active_cells () << std::endl;
pcout << "Number of global cells:" << triangulation.n_global_active_cells () << std::endl;
pcout << " Total Number of globally active cells: "
<< triangulation.n_global_active_cells () << std::endl
<< " Number of degrees of freedom for elasticity: "
<< dof_handler.n_dofs () << std::endl;
}
template<int dim>
void Elasticity<dim>::refine_mesh(const unsigned int load_step)
{
Vector<float> estimated_error_per_cell (triangulation.n_locally_owned_active_cells ());
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1> (fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
local_solution, estimated_error_per_cell);
// Initialize SolutionTransfer object
parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector>
soltrans (dof_handler);
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (
triangulation, estimated_error_per_cell, 0.30, 0.01);
if (triangulation.n_global_levels () >= 4)
{
for (const auto &cell : triangulation.active_cell_iterators_on_level (3))
if (cell->is_locally_owned ())
cell->clear_refine_flag ();
}
// Prepare the triangulation
triangulation.prepare_coarsening_and_refinement ();
// Prepare the SolutionTransfer object for coarsening and refinement
// and give the solution vector that we intend to interpolate later,
soltrans.prepare_for_coarsening_and_refinement (local_solution);
// Do the refinement
triangulation.execute_coarsening_and_refinement ();
// redistribute dofs,
dof_handler.distribute_dofs (fe);
// Recreate locally_owned_dofs and locally_relevant_dofs index sets - ELASTICITY
locally_owned_dofs = dof_handler.locally_owned_dofs ();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs (dof_handler);
distributed_solution_old.reinit (locally_owned_dofs, mpi_communicator);
soltrans.interpolate (distributed_solution_old);
// Apply constraints on the interpolated solution
setup_constraints(load_step);
constraints.distribute (distributed_solution_old);
// Copy distributed_solution_old to local_solution
local_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
local_solution = distributed_solution_old;
}
// -----------------------------
// -------- ELASTICITY ---------
// -----------------------------
template<int dim>
void Elasticity<dim>::setup_constraints (const unsigned int load_step)
{
constraints.clear ();
constraints.reinit (locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
const FEValuesExtractors::Scalar u_x (0);
const FEValuesExtractors::Scalar u_y (1);
const ComponentMask u_x_mask = fe.component_mask (u_x);
const ComponentMask u_y_mask = fe.component_mask (u_y);
applied_load = load_increment*load_step;
const double u_x_values_right = 0.5*applied_load;
const double u_x_values_left = -0.5*applied_load;
const double u_fix = 0.0;
VectorTools::interpolate_boundary_values (dof_handler, 0,
Functions::ConstantFunction<dim> (u_x_values_left, dim),
constraints,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler, 1,
Functions::ConstantFunction<dim> (u_x_values_right, dim),
constraints,
u_x_mask);
VectorTools::interpolate_boundary_values (dof_handler, 2,
Functions::ConstantFunction<dim> (u_fix, dim),
constraints,
u_y_mask);
if constexpr (dim == 3)
{
const FEValuesExtractors::Scalar u_z (2);
const ComponentMask u_z_mask = fe.component_mask (u_z);
VectorTools::interpolate_boundary_values (dof_handler, 4,
Functions::ConstantFunction<dim> (u_fix, dim),
constraints,
u_z_mask);
}
constraints.close ();
}
template<int dim>
void Elasticity<dim>::setup(const unsigned int load_step)
{
TimerOutput::Scope ts(computing_timer, "setup_system");
// Initialize the dofs contained locally as well as those of the near cells
locally_owned_dofs = dof_handler.locally_owned_dofs ();
locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs (dof_handler);
local_solution.reinit (locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
// Initialize the contraints objects for BC and Hanging Nodes during refinement
system_rhs.reinit (locally_owned_dofs, mpi_communicator);
distributed_solution.reinit (locally_owned_dofs, mpi_communicator);
// Initialize hanging nodes and constrains
setup_constraints (load_step);
// Sparsity pattern initialization (where matrix has nonzero entries)
DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern (dof_handler,
dsp,
constraints,
false);
SparsityTools::distribute_sparsity_pattern (dsp,
dof_handler.locally_owned_dofs (),
mpi_communicator,
locally_relevant_dofs);
// Elasticity Matrix
system_matrix.reinit (locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
template<int dim>
void Elasticity<dim>::assemble()
{
TimerOutput::Scope ts (computing_timer, "assembly");
// Clean the current matrices
system_rhs = 0;
system_matrix = 0;
// Initialize the FEValues Objects
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
// Number of DoFs and Gauss points per cell
const unsigned int dofs_per_cell = fe.n_dofs_per_cell ();
// Local Matrices and Vector
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
// Create vector to store the local dof indices
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
// loop over the cells
for (const auto &cell : dof_handler.active_cell_iterators ())
{
if (cell->is_locally_owned ())
{
// Reset the local matrix/rhs of the cell
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
// loop over the gauss point of the cell
for (const unsigned int q_point : fe_values.quadrature_point_indices ())
{
for (const unsigned int i : fe_values.dof_indices ())
{
const unsigned int component_i = fe.system_to_component_index (i).first;
for (const unsigned int j : fe_values.dof_indices ())
{
const unsigned int component_j = fe.system_to_component_index (j).first;
cell_matrix (i, j) += ((fe_values.shape_grad (i, q_point)[component_i]
* lambda (E, nu) *
fe_values.shape_grad (j, q_point)[component_j])
+
(fe_values.shape_grad (i, q_point)[component_j]
* mu (E, nu) *
fe_values.shape_grad (j, q_point)[component_i])
+
((component_i == component_j) ?
(fe_values.shape_grad (i, q_point)
* mu (E, nu) *
fe_values.shape_grad (j, q_point))
: 0)
) * fe_values.JxW (q_point); //
}
// Cell rhs set to zero for the current example
// Modify this entry to account for body forces
cell_rhs (i) += fe_values.shape_value (i, q_point) * 0
* fe_values.JxW (q_point);
// To apply boundary forces you need to loop over the exterior boundaries
// to populate the corresponding entries of the rhs
}
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
system_matrix.compress (VectorOperation::add);
system_rhs.compress (VectorOperation::add);
}
template<int dim>
void Elasticity<dim>::solve_linear_problem()
{
TimerOutput::Scope ts (computing_timer, "solve_linear_system");
// CG
SolverControl solver_control (100000, 1e-12* system_rhs.l2_norm());
SolverCG<PETScWrappers::MPI::Vector> solver (solver_control);
PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
solver.solve (system_matrix,
distributed_solution,
system_rhs,
preconditioner);
pcout << " Solved in " << solver_control.last_step () << " iterations." << std::endl;
constraints.distribute (distributed_solution);
}
// -----------------------------
// ---------- OUTPUT -----------
// -----------------------------
// Solution
template <int dim>
void Elasticity<dim>::output_results (const unsigned int load_step) const
{
std::vector<std::string> displacement_names (dim, "displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
displacement_component_interpretation (dim,
DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> data_out_phasefield;
data_out_phasefield.add_data_vector (dof_handler,
local_solution,
displacement_names,
displacement_component_interpretation);
Vector<double> subdomain (triangulation.n_active_cells ());
for (unsigned int i = 0; i < subdomain.size (); ++i)
{
subdomain (i) = triangulation.locally_owned_subdomain ();
}
data_out_phasefield.add_data_vector (subdomain, "subdomain");
data_out_phasefield.build_patches ();
const std::string pvtu_filename = data_out_phasefield.write_vtu_with_pvtu_record ("./",
"solution",
load_step, mpi_communicator,
2,
0);
if (this_mpi_process == 0)
{
static std::vector<std::pair<double, std::string>> times_and_names;
times_and_names.push_back(std::pair<double, std::string>(load_step, pvtu_filename));
std::ofstream pvd_output("./solution_u.pvd");
DataOutBase::write_pvd_record(pvd_output, times_and_names);
}
}
// Stresses
template <int dim>
void Elasticity<dim>::output_stress(const unsigned int load_step)
{
// Vettori per cella in Vector<float>
Vector<float> s_xx(triangulation.n_active_cells());
Vector<float> s_xy(triangulation.n_active_cells());
Vector<float> s_yy(triangulation.n_active_cells());
Vector<float> s_zz(triangulation.n_active_cells());
Vector<float> s_xz(triangulation.n_active_cells());
Vector<float> s_yz(triangulation.n_active_cells());
FEValues<dim> fe_values(fe, quadrature_formula,
update_gradients | update_JxW_values);
const unsigned int n_q_points = quadrature_formula.size();
const FEValuesExtractors::Vector displacements(0);
std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
fe_values[displacements].get_function_symmetric_gradients(local_solution,
strain_tensor);
double c_s_xx = 0, c_s_xy = 0, c_s_yy = 0;
double c_s_zz = 0, c_s_xz = 0, c_s_yz = 0;
for (unsigned int q = 0; q < n_q_points; ++q)
{
const auto &eps_u = strain_tensor[q];
const double tr = trace(eps_u);
const auto stress = 2.0 * mu(E, nu) * eps_u + lambda(E, nu)
* tr * unit_symmetric_tensor<dim>();
c_s_xx += stress[0][0];
c_s_xy += stress[0][1];
c_s_yy += stress[1][1];
if constexpr (dim == 3)
{
c_s_zz += stress[2][2];
c_s_xz += stress[0][2];
c_s_yz += stress[1][2];
}
}
const unsigned int cell_index = cell->active_cell_index();
s_xx[cell_index] = c_s_xx / n_q_points;
s_xy[cell_index] = c_s_xy / n_q_points;
s_yy[cell_index] = c_s_yy / n_q_points;
if constexpr (dim == 3)
{
s_zz[cell_index] = c_s_zz / n_q_points;
s_xz[cell_index] = c_s_xz / n_q_points;
s_yz[cell_index] = c_s_yz / n_q_points;
}
}
}
DataOut<dim> data_out;
data_out.attach_triangulation(triangulation);
data_out.add_data_vector(s_xx, "s_xx", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_xy, "s_xy", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_yy, "s_yy", DataOut<dim>::type_cell_data);
if constexpr (dim == 3)
{
data_out.add_data_vector(s_zz, "s_zz", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_xz, "s_xz", DataOut<dim>::type_cell_data);
data_out.add_data_vector(s_yz, "s_yz", DataOut<dim>::type_cell_data);
}
data_out.build_patches();
const std::string pvtu_filename = data_out.write_vtu_with_pvtu_record(
"./", "stresses", load_step, mpi_communicator, 4);
if (this_mpi_process == 0)
{
static std::vector<std::pair<double, std::string>> times_and_names;
times_and_names.emplace_back(load_step, pvtu_filename);
std::ofstream pvd_output("./stresses.pvd");
DataOutBase::write_pvd_record(pvd_output, times_and_names);
}
}
// Reaction
template <int dim>
void Elasticity<dim>::reaction(Tensor<1, dim> &reaction_stress,
const types::boundary_id &boundary_id)
{
// Face quadrature formula
QGauss<dim - 1> face_quadrature_formula(5);
// Set up FE face values for vector-valued and scalar-valued solution
FEFaceValues<dim> fe_face_values(fe, face_quadrature_formula,
UpdateFlags(update_values |
update_gradients |
update_quadrature_points |
update_normal_vectors |
update_JxW_values));
const unsigned int n_face_q_points = face_quadrature_formula.size();
// Extract displacements for gradient calculation
const FEValuesExtractors::Vector displacements(0);
std::vector<SymmetricTensor<2, dim>> strain_tensor(n_face_q_points);
// Parallel DoFHandler iterator
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
Tensor<1, dim> reaction_mpi;
// Loop over all cells in the triangulation
for (; cell != endc; ++cell)
{
if (cell->is_locally_owned())
{
// Loop over all faces of the current cell
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
{
// Check if the current face belongs to the specified boundary
if (cell->face(face)->boundary_id() == boundary_id)
{
// Reinitialize FE face values for the current face
fe_face_values.reinit(cell, face);
fe_face_values[displacements].get_function_symmetric_gradients(local_solution, strain_tensor);
// Loop over quadrature points on the face
for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
{
const auto &eps_u = strain_tensor[q_point];
const Tensor<1, dim> &N = fe_face_values.normal_vector(q_point);
const double JxW_f = fe_face_values.JxW(q_point);
// Compute the stress tensor using the strain
const double tr = trace(eps_u);
const auto stress = 2.0 * mu(E, nu) * eps_u + lambda(E, nu)
* tr * unit_symmetric_tensor<dim>();
// Accumulate the contribution to the reaction stress
reaction_mpi += stress * N * JxW_f;
}
}
}
}
}
// Sum the reaction stresses from all MPI ranks and print it
reaction_stress = Utilities::MPI::sum(reaction_mpi, mpi_communicator);
for (unsigned int i = 0; i < dim; ++i)
{
pcout << " R_" << (char)('x' + i) << ": " << reaction_stress[i] << std::endl;
}
}
// -----------------------------
// ----------- RUN -------------
// -----------------------------
template <int dim>
void Elasticity<dim>::run()
{
Timer timer;
timer.start ();
pcout << "Running on " << Utilities::MPI::n_mpi_processes (mpi_communicator)
<< " MPI rank(s)..." << std::endl;
create_mesh();
// Loop over load steps
for (unsigned int load_step = 0; load_step <= num_load_steps; load_step++)
{
pcout << " --- Load increment number : " << load_step << std::endl;
setup (load_step);
assemble ();
solve_linear_problem ();
// Update ghost values after the solution
local_solution = distributed_solution;
local_solution.update_ghost_values();
// Prepare for refinement
distributed_solution_old = local_solution;
// Refine
refine_mesh(load_step);
pcout << " Total New Number of globally active cells: "
<< triangulation.n_global_active_cells () << std::endl
<< " New Number of degrees of freedom for elasticity: "
<< dof_handler.n_dofs () << std::endl;
// Save results
output_results (load_step);
output_stress (load_step);
reaction(reaction_stress, 1);
// Print the table of cpmputing times
computing_timer.print_summary ();
computing_timer.reset ();
pcout << std::endl;
}
timer.stop ();
pcout << "Total run time: " << timer.wall_time () << " seconds." << std::endl;
}
} // End of namespace Elasticity
// -----------------------------
// ---------- MAIIN ------------
// -----------------------------
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace elas;
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
Elasticity<3> elasticity;
elasticity.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what ()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
return 1;
}
return 0;
}