600 *
const unsigned int component = 0)
const override
641 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
645 *
else if (p[0] >= -0.5 && p[0] < 0.0)
647 *
else if (p[0] >= 0.0 && p[0] < 0.5)
659 * <a name=
"step_41-ImplementationofthecodeObstacleProblemcodeclass"></a>
668 * <a name=
"step_41-ObstacleProblemObstacleProblem"></a>
687 * <a name=
"step_41-ObstacleProblemmake_grid"></a>
703 *
std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
713 * <a name=
"step_41-ObstacleProblemsetup_system"></a>
727 *
dof_handler.distribute_dofs(fe);
730 *
std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
738 *
constraints.close();
743 *
system_matrix.reinit(
dsp);
757 * diagonal, and in the following then first compute all of this as a
758 * matrix and then extract the diagonal elements for later use:
761 * TrilinosWrappers::SparseMatrix mass_matrix;
762 * mass_matrix.reinit(dsp);
763 * assemble_mass_matrix_diagonal(mass_matrix);
764 * diagonal_of_mass_matrix.reinit(solution_index_set);
765 * for (unsigned int j = 0; j < solution.size(); ++j)
766 * diagonal_of_mass_matrix(j) = mass_matrix.diag_element(j);
773 * <a name="step_41-ObstacleProblemassemble_system"></a>
774 * <h4>ObstacleProblem::assemble_system</h4>
778 * This function at once assembles the system matrix and right-hand-side and
779 * applied the constraints (both due to the active set as well as from
780 * boundary values) to our system. Otherwise, it is functionally equivalent
781 * to the corresponding function in, for example, @ref step_4 "step-4".
785 * void ObstacleProblem<dim>::assemble_system()
787 * std::cout << " Assembling system..." << std::endl;
792 * const QGauss<dim> quadrature_formula(fe.degree + 1);
793 * RightHandSide<dim> right_hand_side;
795 * FEValues<dim> fe_values(fe,
796 * quadrature_formula,
797 * update_values | update_gradients |
798 * update_quadrature_points | update_JxW_values);
800 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
801 * const unsigned int n_q_points = quadrature_formula.size();
803 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
804 * Vector<double> cell_rhs(dofs_per_cell);
806 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
808 * for (const auto &cell : dof_handler.active_cell_iterators())
810 * fe_values.reinit(cell);
814 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
815 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
817 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
818 * cell_matrix(i, j) +=
819 * (fe_values.shape_grad(i, q_point) *
820 * fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
823 * (fe_values.shape_value(i, q_point) *
824 * right_hand_side.value(fe_values.quadrature_point(q_point)) *
825 * fe_values.JxW(q_point));
828 * cell->get_dof_indices(local_dof_indices);
830 * constraints.distribute_local_to_global(cell_matrix,
844 * <a name="step_41-ObstacleProblemassemble_mass_matrix_diagonal"></a>
845 * <h4>ObstacleProblem::assemble_mass_matrix_diagonal</h4>
849 * The next function is used in the computation of the diagonal mass matrix
850 * @f$B@f$ used to scale variables in the active set method. As discussed in the
851 * introduction, we get the mass matrix to be diagonal by choosing the
870 * really difficult, but not the point here, and so we simply assert at the
871 * top of the function that our implicit assumption about the finite element
872 * is in fact satisfied.
876 * void ObstacleProblem<dim>::assemble_mass_matrix_diagonal(
877 * TrilinosWrappers::SparseMatrix &mass_matrix)
879 * Assert(fe.degree == 1, ExcNotImplemented());
881 * const QTrapezoid<dim> quadrature_formula;
882 * FEValues<dim> fe_values(fe,
883 * quadrature_formula,
884 * update_values | update_JxW_values);
886 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
887 * const unsigned int n_q_points = quadrature_formula.size();
889 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
890 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
892 * for (const auto &cell : dof_handler.active_cell_iterators())
894 * fe_values.reinit(cell);
897 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
898 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
899 * cell_matrix(i, i) +=
900 * (fe_values.shape_value(i, q_point) *
901 * fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
903 * cell->get_dof_indices(local_dof_indices);
905 * constraints.distribute_local_to_global(cell_matrix,
915 * <a name="step_41-ObstacleProblemupdate_solution_and_constraints"></a>
916 * <h4>ObstacleProblem::update_solution_and_constraints</h4>
920 * In a sense, this is the central function of this program. It updates the
921 * active set of constrained degrees of freedom as discussed in the
922 * introduction and computes an AffineConstraints object from it that can then
923 * be used to eliminate constrained degrees of freedom from the solution of
924 * the next iteration. At the same time we set the constrained degrees of
925 * freedom of the solution to the correct value, namely the height of the
930 * Fundamentally, the function is rather simple: We have to loop over all
931 * degrees of freedom and check the sign of the function @f$\Lambda^k_i +
932 * c([BU^k]_i - G_i) = \Lambda^k_i + cB_i(U^k_i - [g_h]_i)@f$ because in our
933 * case @f$G_i = B_i[g_h]_i@f$. To this end, we use the formula given in the
934 * introduction by which we can compute the Lagrange multiplier as the
935 * residual of the original linear system (given via the variables
936 * <code>complete_system_matrix</code> and <code>complete_system_rhs</code>.
937 * At the top of this function, we compute this residual using a function
938 * that is part of the matrix classes.
942 * void ObstacleProblem<dim>::update_solution_and_constraints()
944 * std::cout << " Updating active set..." << std::endl;
946 * const double penalty_parameter = 100.0;
948 * TrilinosWrappers::MPI::Vector lambda(
949 * complete_index_set(dof_handler.n_dofs()));
950 * complete_system_matrix.residual(lambda, solution, complete_system_rhs);
954 * compute contact_force[i] = - lambda[i] * diagonal_of_mass_matrix[i]
957 * contact_force = lambda;
958 * contact_force.scale(diagonal_of_mass_matrix);
959 * contact_force *= -1;
963 * The next step is to reset the active set and constraints objects and to
964 * start the loop over all degrees of freedom. This is made slightly more
965 * complicated by the fact that we can't
just loop over all elements
of
979 * and so we add an assertion that makes sure that we only deal with
980 * elements for which all degrees of freedom are located in vertices to
981 * avoid tripping ourselves with non-functional code in case someone wants
982 * to play with increasing the polynomial degree of the solution.
986 * The price to pay for having to loop over cells rather than DoFs is that
987 * we may encounter some degrees of freedom more than once, namely each
988 * time we visit one of the cells adjacent to a given vertex. We will
989 * therefore have to keep track which vertices we have already touched and
994 *
constraints.
clear();
998 *
std::vector<bool>
dof_touched(dof_handler.n_dofs(),
false);
1000 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1003 *
Assert(dof_handler.get_fe().n_dofs_per_cell() == cell->n_vertices(),
1004 *
ExcNotImplemented());
1006 *
const unsigned int dof_index = cell->vertex_dof_index(v, 0);
1054 *
std::cout <<
" Size of active set: " <<
active_set.n_elements()
1057 *
std::cout <<
" Residual of the non-contact part of the system: "
1058 *
<<
lambda.l2_norm() << std::endl;
1071 *
constraints.close();
1077 * <a name=
"step_41-ObstacleProblemsolve"></a>
1078 * <
h4>ObstacleProblem::solve</
h4>
1094 *
std::cout <<
" Solving system..." << std::endl;
1099 *
precondition.initialize(system_matrix);
1101 *
solver.solve(system_matrix, solution,
system_rhs, precondition);
1102 *
constraints.distribute(solution);
1114 * <a name=
"step_41-ObstacleProblemoutput_results"></a>
1123 *
template <
int dim>
1126 *
std::cout <<
" Writing graphical output..." << std::endl;
1135 *
data_out.attach_dof_handler(dof_handler);
1136 *
data_out.add_data_vector(solution,
"displacement");
1140 *
data_out.build_patches();
1152 * <a name=
"step_41-ObstacleProblemrun"></a>
1153 * <
h4>ObstacleProblem::run</
h4>
1177 *
template <
int dim>
1186 *
std::cout <<
"Newton iteration " <<
iteration << std::endl;
1205 *
std::cout << std::endl;
1214 * <a name=
"step_41-Thecodemaincodefunction"></a>
1228 *
using namespace dealii;
1229 *
using namespace Step41;
1241 *
"This program can only be run in serial, use ./step-41"));
1246 *
catch (std::exception &exc)
1248 *
std::cerr << std::endl
1250 *
<<
"----------------------------------------------------"
1252 *
std::cerr <<
"Exception on processing: " << std::endl
1253 *
<< exc.what() << std::endl
1254 *
<<
"Aborting!" << std::endl
1255 *
<<
"----------------------------------------------------"
1262 *
std::cerr << std::endl
1264 *
<<
"----------------------------------------------------"
1266 *
std::cerr <<
"Unknown exception!" << std::endl
1267 *
<<
"Aborting!" << std::endl
1268 *
<<
"----------------------------------------------------"
13255,399 constrained degrees of freedom at that point). The algebraic
1326precondition is apparently working nicely since we only need 4-6 CG
1327iterations to solve the linear system (although this also has a lot to
1328do with the fact that we are not asking for very high accuracy of the
1331More revealing is to look at a sequence of graphical output files
1332(every third step is shown, with the number of the iteration in the
1335<table align="center">
1341 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
1344 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
1347 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
1355 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
1358 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
1361 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
1369 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
1372 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
1375 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
1383 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
1386 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
1389 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
1397 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
1400 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
1403 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
1411 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
1414 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
1417 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
1425 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
1428 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
1431 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
1436The pictures show that in the first step, the solution (which has been
1437computed without any of the constraints active) bends through so much
1438that pretty much every interior point has to be bounced back to the
1439stairstep function, producing a discontinuous solution. Over the
1440course of the active set iterations, this unphysical membrane shape is
1441smoothed out, the contact with the lower-most stair step disappears,
1442and the solution stabilizes.
1444In addition to this, the program also outputs the values of the
1445Lagrange multipliers. Remember that these are the contact forces and
1446so should only be positive on the contact set, and zero outside. If,
1447on the other hand, a Lagrange multiplier is negative in the active
1448set, then this degree of freedom must be removed from the active
1449set. The following pictures show the multipliers in iterations 1, 9
1450and 18, where we use red and browns to indicate positive values, and
1451blue for negative values.
1453<table align="center">
1456 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
1459 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
1462 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
1478It is easy to see that the positive values converge nicely to moderate
1479values in the interior of the contact set and large upward forces at
1480the edges of the steps, as one would expect (to support the large
1481curvature of the membrane there); at the fringes of the active set,
1482multipliers are initially negative, causing the set to shrink until,
1483in iteration 18, there are no more negative multipliers and the
1484algorithm has converged.
1488<a name="step-41-extensions"></a>
1489<a name="step_41-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1492As with any of the programs of this tutorial, there are a number of
1493obvious possibilities for extensions and experiments. The first one is
1494clear: introduce adaptivity. Contact problems are prime candidates for
1495adaptive meshes because the solution has lines along which it is less
1496regular (the places where contact is established between membrane and
1497obstacle) and other areas where the solution is very smooth (or, in
1498the present context, constant wherever it is in contact with the
1499obstacle). Adding this to the current program should not pose too many
1500difficulties, but it is not trivial to find a good error estimator for
1503A more challenging task would be an extension to 3d. The problem here
1504is not so much to simply make everything run in 3d. Rather, it is that
1505when a 3d body is deformed and gets into contact with an obstacle,
1506then the obstacle does not act as a constraining body force within the
1507domain as is the case here. Rather, the contact force only acts on the
1508boundary of the object. The inequality then is not in the differential
1509equation but in fact in the (Neumann-type) boundary conditions, though
1510this leads to a similar kind of variational
1511inequality. Mathematically, this means that the Lagrange multiplier
1512only lives on the surface, though it can of course be extended by zero
1513into the domain if that is convenient. As in the current program, one
1514does not need to form and store this Lagrange multiplier explicitly.
1516A further interesting problem for the 3d case is to consider contact problems
1517with friction. In almost every mechanical process friction has a big influence.
1518For the modelling we have to take into account tangential stresses at the contact
1519surface. Also we have to observe that friction adds another nonlinearity to
1522Another nontrivial modification is to implement a more complex constitutive
1523law like nonlinear elasticity or elasto-plastic material behavior.
1524The difficulty here is to handle the additional nonlinearity arising
1525through the nonlinear constitutive law.
1528<a name="step_41-PlainProg"></a>
1529<h1> The plain program</h1>
1530@include "step-41.cc"
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
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)
std::vector< value_type > preserve(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
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)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
constexpr unsigned int invalid_unsigned_int
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation