This program was contributed by Jean-Paul Pelteret <jppelteret@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
Pictures from this code gallery program
Annotated version of README.md
Overview
Many rubber-like materials are not only near-incompressible in nature, but also exhibit a viscoelastic response (within the tested load and time scales). In this example, we extend the near-incompressible rate-independent constitutive used in step-44 (which implements three-field quasi-static quasi-incompressible finite elasticity) to one that displays rate-dependent behavior. It may seem that there is a contradiction of terms here, so lets clarify by saying the the problem remains "quasi-static" in the sense that inertial terms remain insignificant, even though the material response itself is rate-dependent. This implies that, for these fictitious material parameters, it is assumed that the timescale for material relaxation is much longer than that of elastic wave propagation.
We've also taken the opportunity to extend the code first shown in step-44 to parallel (the primary reason for this contribution), using Metis
as a grid partitioner and Trilinos
for linear algebra. As a motivation as to why one might still choose to use Metis
(also associated with parallel::shared::Triangulation in step-18, although this triangulation is not used in this instance) over p4est
(also associated with parallel::distributed::Triangulation) as a grid partitioner, at this point in time it is not possible to use the hp
finite-element in conjunction with the distributed grid, meaning that this code could not, for example, be readily extended to the application shown in
J-P. V. Pelteret, D. Davydov, A. McBride, D. Vu, P. Steinmann, Computational electro-elasticity and magneto-elasticity for quasi-incompressible media immersed in free space. International Journal for Numerical Methods in Engineering, 2016, 108, 1307-1342. DOI: 10.1002/nme.5254
The discerning reader will observe that we've chosen to employ deal.II
's built in solvers as opposed to using Trilinos
solvers. This is because the system matrices K_Jp
and K_pJ
, although block diagonal and well conditioned, and for some reason (perhaps pertaining to the negative definite nature of these blocks, or that the entries are very small in magnitude) Trilinos
solvers are not sufficiently robust to compute inverse matrix-vector multiplication with. We do stress, however, that to date no great attempt has been made by the author to overcome this issue other than by making an entirely different choice of solver.
Finite deformation of a thin strip with a hole.
Various permutations of the problem of an elastomeric strip with a centered cut-out can be found in the literature for solid mechanics, in particular (but not limited to) that pertaining to
incompressible elasticity elasto-plasticity electro-elasticity thermo-elasticity.
Here we implement another permutation (one that is not necessarily benchmarked elsewhere), simply for demonstration purposes. The basic problem configuration is summarized in the following image.

A thin strip of material with a circular hole is (in 3d
) constrained in the Z
direction and loaded in the direction of its long edge. In our implementation, this load is applied to the +Y
surface and may either be displacement control (a Dirichlet
condition) or a traction load (a Neumann
boundary condition). Due to the symmetry of both the geometry and load, the problem can be simplified by modeling only an eighth of the geometry in 3d
or a quarter in 2d
. By doing so, it it necessary to then implement symmetry conditions on the surfaces coincident with the X-Z
and Y-Z
planes (and the X-Y
plane in 3d
). The +X
surface, and that of the hole itself, remain traction-free.
In three dimensions, the geometry (and a potential partitioning over 8 processors) looks as follows:

Note that, for this particular formulation, the two-dimensional case corresponds to neither plane-strain nor plane-stress conditions.
Requirements
Version 8.5.0
or greater of deal.II
C++11
and MPI
must be enabled The following packages must also be enabled: Metis
Trilinos
Compiling and running
Similar to the example programs, run
cmake -DDEAL_II_DIR=/path/to/deal.II .
in this directory to configure the problem.
You can switch between debug and release mode by calling either
or
The problem may then be run in serial mode with
and in parallel (in this case, on 4
processors) with
mpirun -np 4 ./viscoelastic_strip_with_hole
This program can be run in 2d
or 3d
; this choice can be made by making the appropriate changes in the main()
function.
Reference for this work
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation: J-P. V. Pelteret, The deal.II code gallery: Quasi-static Finite-Strain Quasi-incompressible Visco-elasticity, 2017. DOI: 10.5281/zenodo.437604 
Acknowledgements
The support of this work by the European Research Council (ERC) through the Advanced Grant 289049 MOCOPOLY is gratefully acknowledged by the author.
Recommended Literature
C. Miehe (1994), Aspects of the formulation and finite element implementation of large strain isotropic elasticity. International Journal for Numerical Methods in Engineering 37 , 12, 1981-2004. DOI: 10.1002/nme.1620371202; G.A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons. ISBN: 978-0-471-82319-3; P. Wriggers (2008), Nonlinear finite element methods, Springer. DOI: 10.1007/978-3-540-71001-1; T.J.R. Hughes (2000), The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover. ISBN: 978-0486411811
The derivation of the finite-element problem, namely the definition and linearization of the residual and their subsequent discretization are quite lengthy and involved. However, this is already detailed in step-44, some of the aforementioned literature, as well as J-P. V. Pelteret, A computational neuromuscular model of the human upper airway with application to the study of obstructive sleep apnoea. PhD Thesis, University of Cape Town, 2013. http://hdl.handle.net/11427/9519
and need not be repeated here. As for the viscoelastic constitutive law (which satisfies the dissipation inequality through the definition of an internal variable that we denote as Q
in the code), this is derived and presented in detail in C. Linder, M. Tkachuk & C. Miehe, A micromechanically motivated diffusion-based transient network model and its incorporation into finite rubber viscoelasticity. Journal of the Mechanics and Physics of Solids, 2011, 59, 2134-2156. DOI: 10.1016/j.jmps.2011.05.005
In particular, the relevant equations from Linder et al.'s work that are implemented in this work are equations 47, 54 and 56. Note that the time discretization for the rate-dependent internal variable for this single dissipative mechanism is only first-order accurate.
Results
To begin, here is a comparison of the initial grid for the 2d
version of the problem (mirrored about the x-axis and rotated 90
degrees anti-clockwise)

and of the final, displaced grid after the load has been applied and the material is in a (near-completely) relaxed state.

These results, as well as those that follow, were produced using the following material properties: The Poisson ratio is 0.4995
The elastic shear modulus is 80MPa
The viscoelastic shear modulus is 160MPa
The viscoelastic relaxation time is 2.5s
while the boundary conditions were configured in the following manner: The "driving" boundary condition on the short-side (+Y
face) is of the Neumann
variety The applied hydrostatic pressure is -150Pa
(i.e. a tensile load) The load ramp time is 1s
The following chart highlights the displacement of several vertices and clearly illustrates the viscoelastic nature of the material.

During the initial phase, the load is applied over a period of time much shorter than the material's characteristic relaxation time. The material therefore exhibits a very stiff response, and relaxes as the load remains constant for the remainder of the simulation. This deformation that occurs under constant load is commonly known as creep.
We've been lazy and stopped the simulation slightly prematurely, but it is clear that the material's displacement is moving asymptotically towards a equilibrium solution. You can also check what the true resting load state is by removing the dissipative mechanism (setting its shear modulus to zero) and rerunning the simulation with the material then exhibiting rate-independent behavior. Note that in this case, the number of time step over which the load is applied must likely be increased in order to ensure stability of the nonlinear problem.
Animations
Just in case you found the presented chart a little dry and undigestible, below are a couple of animations that demonstrate the viscoelastic nature of the material in a more visually appealing manner.
Animation of the evolution of the displacement field** 
Animation of the evolution of the pressure field** 
{
"Driver boundary condition for the problem");
"Positive stretch applied length-ways to the strip");
"Hydrostatic pressure applied (in the referential configuration) to the interior surface of the hole");
"Total time over which the stretch/pressure is ramped up");
{
{
"Displacement system polynomial order");
"Gauss quadrature order");
{
{
"Total sample thickness");
"A geometric factor affecting the discretisation near the hole");
prm.
declare_entry(
"Number of subdivisions in cross-section",
"2",
"A factor defining the number of initial grid subdivisions in the cross-section");
"A factor defining the number of initial grid subdivisions through the thickness");
"Global refinement level");
"Global grid scaling factor");
{
{
"Elastic shear modulus");
"Viscous shear modulus");
"Viscous relaxation time");
{
{
"Type of solver used to solve the linear system");
"Linear solver residual (scaled by residual norm)");
"Linear solver iterations (multiples of the system matrix size)");
{
{
"Number of Newton-Raphson iterations allowed");
"Displacement error tolerance");
"Force residual tolerance");
{
{
{
{
{
BoundaryConditions::declare_parameters(prm);
FESystem::declare_parameters(prm);
Geometry::declare_parameters(prm);
Materials::declare_parameters(prm);
LinearSolver::declare_parameters(prm);
NonlinearSolver::declare_parameters(prm);
Time::declare_parameters(prm);
{
BoundaryConditions::parse_parameters(prm);
FESystem::parse_parameters(prm);
Geometry::parse_parameters(prm);
Materials::parse_parameters(prm);
LinearSolver::parse_parameters(prm);
NonlinearSolver::parse_parameters(prm);
Time::parse_parameters(prm);
:
:
{
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
long int get_integer(const std::string &entry_string) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
std::string get(const std::string &entry_string) const
double get_double(const std::string &entry_name) const
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Linder2011 eq 54 Assumes first-oder backward Euler time discretisation
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Elastic Neo-Hookean + Linder2011 eq 47
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
Elastic Neo-Hookean + Linder2011 eq 56
{
parameters.mu_e, parameters.nu_e,
parameters.mu_v, parameters.tau_v,
{
std::shared_ptr< Material_Compressible_Three_Field_Linear_Viscoelastic<dim> >
material;
std::pair<unsigned int, double>
const double current_time)
const;
typename ActiveSelector::active_cell_iterator active_cell_iterator
std::vector< index_type > data
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Parallel communication
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
const unsigned int degree;
const unsigned int dofs_per_cell;
static const unsigned int n_blocks = 3;
static const unsigned int n_components = dim + 2;
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Block data
DoF index data
const unsigned int n_q_points;
norm(1.0),
u(1.0), p(1.0), J(1.0)
{
std::pair<double, std::pair<double,double> >
:
time(parameters.end_time, parameters.delta_t),
dofs_per_cell (fe.dofs_per_cell),
Assert(dim==2 || dim==3, ExcMessage(
"This problem only works in 2 or 3 space dimensions."));
numbers::NumberTraits< Number >::real_type norm() const
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Some points for post-processing
p[1] = parameters.length/2.0;
p[1] = parameters.hole_diameter/2.0;
p[0] = parameters.hole_diameter/2.0;
p[0] = parameters.width/2.0;
while (time.current() < time.end()+0.01*time.get_delta_t())
pcout <<
"\n\n*** Spatial position history for tracked vertices ***" << std::endl;
for (
unsigned int t=0; t<
real_time.size(); ++t)
for (
unsigned int d=0;
d<dim; ++
d)
pcout <<
"Point " << p <<
" [" <<
d <<
"]";
pcout << std::setprecision(6);
ExcMessage(
"Vertex not tracked at each timestep"));
for (
unsigned int d=0;
d<dim; ++
d)
std::vector<types::global_dof_index> local_dof_indices;
:
local_dof_indices(dofs_per_cell)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Integration helper
Quadrature point solution
Shape function values and gradients
std::vector<std::vector<double> >
Nx;
std::vector<std::vector<Tensor<2, dim> > >
grad_Nx;
std::vector<std::vector<SymmetricTensor<2, dim> > >
symm_grad_Nx;
:
const unsigned int n_dofs_per_cell =
Nx[0].size();
for (
unsigned int k = 0;
k < n_dofs_per_cell; ++
k)
const double tol = 1
e-12;
parameters.hole_diameter/2.0,
parameters.n_repetitions_xy,
parameters.hole_division_fraction);
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Clear boundary ID's
if (cell->face(face)->at_boundary())
cell->face(face)->set_all_boundary_ids(0);
Set boundary IDs and and manifolds
if (cell->face(face)->at_boundary())
Set boundary IDs
if (
std::abs(cell->face(face)->center()[0] - 0.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
else if (
std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
else if (
std::abs(cell->face(face)->center()[1] - 0.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
else if (
std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
if (
std::abs(cell->vertex(vertex).distance(
centre) - parameters.hole_diameter/2.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Set manifold IDs
if (
std::abs(cell->vertex(vertex).distance(
centre) - parameters.hole_diameter/2.0) < tol)
cell->face(face)->set_manifold_id(parameters.manifold_id_hole);
const double tol = 1
e-12;
parameters.hole_diameter/2.0,
parameters.n_repetitions_xy,
parameters.hole_division_fraction);
parameters.n_repetitions_z+1,
parameters.thickness/2.0,
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
Clear boundary ID's
if (cell->face(face)->at_boundary())
cell->face(face)->set_all_boundary_ids(0);
Set boundary IDs and and manifolds
if (cell->face(face)->at_boundary())
Set boundary IDs
if (
std::abs(cell->face(face)->center()[0] - 0.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X);
else if (
std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X);
else if (
std::abs(cell->face(face)->center()[1] - 0.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y);
else if (
std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y);
else if (
std::abs(cell->face(face)->center()[2] - 0.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Z);
else if (
std::abs(cell->face(face)->center()[2] - parameters.thickness/2.0) < tol)
cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Z);
Project the cell vertex to the XY plane and test the distance from the cylinder axis
cell->face(face)->set_boundary_id(parameters.boundary_id_hole);
Set manifold IDs
Project the cell vertex to the XY plane and test the distance from the cylinder axis
Set manifold ID on face and edges
cell->face(face)->set_all_manifold_ids(parameters.manifold_id_hole);
{
std::set<typename Triangulation<2>::active_cell_iterator >
cells_to_remove;
void hyper_cube_with_cylindrical_hole(Triangulation< dim, spacedim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
Remove all cells that are not in the first quadrant
if (cell->center()[0] < 0.0 || cell->center()[1] < 0.0)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
Subdivide the plate so that we're left one cell to remove (we'll replace this with the plate with the hole) and then make the rest of the subdivisions so that we're left with cells with a decent aspect ratio
for (
unsigned int s=0; s<
n_subs; ++s)
for (
unsigned int s=0; s<
n_subs; ++s)
std::set<typename Triangulation<2>::active_cell_iterator >
cells_to_remove;
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Remove all cells that are in the first quadrant
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
Attach a manifold to the curved boundary and refine Note: We can only guarantee that the vertices sit on the curve, so we must test with their position instead of the cell centre.
if (cell->face(face)->at_boundary())
cell->face(face)->set_manifold_id(10);
{
timer.enter_subsection(
"Setup system");
pcout <<
"Setting up linear system..." << std::endl;
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
Partition triangulation
dof_handler.distribute_dofs(fe);
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
Count DoFs in each block
locally_owned_dofs.
clear();
for (
unsigned int b=0; b<n_blocks; ++b)
for (
unsigned int p=0; p<n_mpi_processes; ++p)
pcout <<
")" << std::endl;
<<
" Number of degrees of freedom: " << dof_handler.n_dofs()
for (
unsigned int p=0; p<n_mpi_processes; ++p)
<< (
DoFTools::count_dofs_with_subdomain_association (dof_handler,p));
pcout <<
")" << std::endl;
<<
" Number of degrees of freedom per block: "
for (
unsigned int ii = 0;
ii < n_components; ++
ii)
for (
unsigned int jj = 0;
jj < n_components; ++
jj)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
We then set up storage vectors
timer.leave_subsection();
for (
unsigned int k = 0;
k < fe.dofs_per_cell; ++
k)
const unsigned int k_group = fe.system_to_base_index(
k).first.first;
pcout <<
"Setting up quadrature point data..." << std::endl;
dof_handler.begin_active()),
for (; cell!=
endc; ++cell)
const std::vector<std::shared_ptr<PointHistory<dim> > >
lqph =
{
<<
"Timestep " << time.get_timestep() <<
" @ "
<< time.current() <<
"s of "
<< time.end() <<
"s" << std::endl;
pcout <<
" CONVERGED! " << std::endl;
const std::pair<unsigned int, double>
pcout <<
" | " << std::fixed << std::setprecision(3) << std::setw(7)
pcout << std::string(132,
'_') << std::endl;
<<
" | LIN_IT LIN_RES RES_NORM "
<<
" RES_U RES_P RES_J NU_NORM "
<<
" NU_U NU_P NU_J " << std::endl;
pcout << std::string(132,
'_') << std::endl;
{
pcout << std::string(132,
'_') << std::endl;
pcout <<
"Relative errors:" << std::endl
<<
"Dilatation:\t" <<
error_dil.first << std::endl
std::pair<double,std::pair<double,double> >
dof_handler.begin_active()),
for (; cell !=
endc; ++cell)
const std::vector<std::shared_ptr<const PointHistory<dim> > >
lqph =
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
Sum across all processors
{
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Construct a residual vector that has the values for all of its constrained DoFs set to zero.
Construct a update vector that has the values for all of its constrained DoFs set to zero.
Cell interpolation -> Ghosted vector
{
timer.enter_subsection(
"Assemble system");
pcout <<
" ASM_SYS " << std::flush;
dof_handler.begin_active()),
for (; cell !=
endc; ++cell)
timer.leave_subsection();
{
constraints.distribute_local_to_global(
data.cell_matrix,
data.cell_rhs,
Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
scratch.fe_values_ref.reinit(cell);
cell->get_dof_indices(
data.local_dof_indices);
const std::vector<std::shared_ptr<const PointHistory<dim> > >
lqph =
Assert(
lqph.size() == n_q_points, ExcInternalError());
@ update_normal_vectors
Normal vectors.
Update quadrature point solution
scratch.fe_values_ref[
u_fe].get_function_gradients(scratch.solution_total,
scratch.solution_grads_u_total);
scratch.fe_values_ref[
p_fe].get_function_values(scratch.solution_total,
scratch.solution_values_p_total);
scratch.fe_values_ref[
J_fe].get_function_values(scratch.solution_total,
scratch.solution_values_J_total);
Update shape functions and their gradients (push-forward)
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
const unsigned int k_group = fe.system_to_base_index(
k).first.first;
const std::vector<double> &
Nx = scratch.Nx[
q_point];
const std::vector<Tensor<2, dim> > &
grad_Nx = scratch.grad_Nx[
q_point];
const double JxW = scratch.fe_values_ref.JxW(
q_point);
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
const unsigned int component_i = fe.system_to_component_index(i).first;
const unsigned int i_group = fe.system_to_base_index(i).first.first;
for (
unsigned int j = 0;
j <= i; ++
j)
const unsigned int component_j = fe.system_to_component_index(
j).first;
const unsigned int j_group = fe.system_to_base_index(
j).first.first;
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
for (
unsigned int j = i + 1;
j < dofs_per_cell; ++
j)
if (parameters.driver ==
"Neumann")
if (cell->face(face)->at_boundary() ==
true
&& cell->face(face)->boundary_id() == parameters.boundary_id_plus_Y)
scratch.fe_face_values_ref.reinit(cell, face);
scratch.fe_face_values_ref.normal_vector(
f_q_point);
/ (parameters.scale * parameters.scale);
const double time_ramp = (time.current() < parameters.load_time ?
time.current() / parameters.load_time : 1.0);
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
fe.system_to_base_index(i).first.first;
fe.system_to_component_index(i).first;
scratch.fe_face_values_ref.shape_value(i,
const double JxW = scratch.fe_face_values_ref.JxW(
{
pcout <<
" CST " << std::flush;
const int boundary_id = parameters.boundary_id_minus_X;
const int boundary_id = parameters.boundary_id_minus_Y;
const int boundary_id = parameters.boundary_id_minus_Z;
if (parameters.driver ==
"Dirichlet")
if (time.current() < parameters.load_time+0.01*time.get_delta_t())
const double delta_length = parameters.length*(parameters.stretch - 1.0)*parameters.scale;
const unsigned int n_stretch_steps =
static_cast<unsigned int>(parameters.load_time/time.get_delta_t());
std::pair<unsigned int, double>
{
timer.enter_subsection(
"Linear solver");
pcout <<
" SLV " << std::flush;
LA::PreconditionJacobi::AdditionalData());
* parameters.max_iterations_lin),
LA::PreconditionAMG::AdditionalData(
(parameters.poly_degree > 1 )) );
* parameters.max_iterations_lin),
1.0e-30, parameters.tol_lin);
timer.leave_subsection();
timer.enter_subsection(
"Linear solver postprocessing");
timer.leave_subsection();
dof_handler.begin_active()),
for (; cell !=
endc; ++cell)
const std::vector<std::shared_ptr<PointHistory<dim> > >
lqph =
:
auto cell = this->dofs->begin_active();
while ((cell != this->dofs->end()) &&
(cell->subdomain_id() != subdomain_id))
{
const double current_time)
const
typename DataOut_DoFData< dim, dim, spacedim, spacedim >::cell_iterator cell_iterator
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
unsigned int subdomain_id
Output -> Ghosted vector
— Additional data —
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(dim,
@ component_is_part_of_vector
Can't use filtered iterators here because the cell count "c" is incorrect for the parallel case
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=
endc; ++cell, ++c)
if (cell->subdomain_id() != this_mpi_process)
continue;
material_id(c) =
static_cast<int>(cell->material_id());
std::vector<std::string>
solution_name(n_components,
"solution_");
std::vector<std::string>
residual_name(n_components,
"residual_");
for (
unsigned int c=0; c<n_components; ++c)
data_out.attach_dof_handler(dof_handler);
data_component_interpretation);
data_out.add_data_vector(residual,
data_component_interpretation);
data_out.add_data_vector (material_id,
"material_id");
data_out.build_patches(degree);
<< (std::to_string(dim) +
"d")
<< "."
<< ".vtu";
<< (std::to_string(dim) +
"d")
<< ".pvtu";
{
<< (std::to_string(dim) +
"d")
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Write out main data file
data_out.write_vtu(output);
Collection of files written in parallel This next set of steps should only be performed by master process
if (this_mpi_process == 0)
List of all files written out at this timestep by all processors
for (
unsigned int p=0; p < n_mpi_processes; ++p)
Time dependent data master file
dof_handler.begin_active()),
for (; cell !=
endc; ++cell)
Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());
if (cell->vertex(v).distance(
pt_ref) < 1e-6*parameters.scale)
for (
unsigned int d=0;
d<dim; ++
d)
for (
unsigned int d=0;
d<dim; ++
d)
{
const unsigned int dim = 2;
catch (std::exception &exc)
std::cerr << std::endl << std::endl
<<
"----------------------------------------------------"
std::cerr <<
"Exception on processing: " << std::endl << exc.what()
<< std::endl <<
"Aborting!" << std::endl
<<
"----------------------------------------------------"
std::cerr << std::endl << std::endl
<<
"----------------------------------------------------"
std::cerr <<
"Unknown exception!" << std::endl <<
"Aborting!"
<<
"----------------------------------------------------"
void write_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > ×_and_names)