Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-41.h
Go to the documentation of this file.
1 ,
603  * const unsigned int component = 0) const override
604  * {
605  * (void)component;
606  * AssertIndexRange(component, 1);
607  *
608  * return -10;
609  * }
610  * };
611  *
612  *
613  *
614  * template <int dim>
615  * class BoundaryValues : public Function<dim>
616  * {
617  * public:
618  * virtual double value(const Point<dim> & /*p*/,
619  * const unsigned int component = 0) const override
620  * {
621  * (void)component;
622  * AssertIndexRange(component, 1);
623  *
624  * return 0;
625  * }
626  * };
627  *
628  *
629  *
630  * @endcode
631  *
632  * We describe the obstacle function by a cascaded barrier (think: stair
633  * steps):
634  *
635  * @code
636  * template <int dim>
637  * class Obstacle : public Function<dim>
638  * {
639  * public:
640  * virtual double value(const Point<dim> & p,
641  * const unsigned int component = 0) const override
642  * {
643  * (void)component;
644  * Assert(component == 0, ExcIndexRange(component, 0, 1));
645  *
646  * if (p(0) < -0.5)
647  * return -0.2;
648  * else if (p(0) >= -0.5 && p(0) < 0.0)
649  * return -0.4;
650  * else if (p(0) >= 0.0 && p(0) < 0.5)
651  * return -0.6;
652  * else
653  * return -0.8;
654  * }
655  * };
656  *
657  *
658  *
659  * @endcode
660  *
661  *
662  * <a name="ImplementationofthecodeObstacleProblemcodeclass"></a>
663  * <h3>Implementation of the <code>ObstacleProblem</code> class</h3>
664  *
665 
666  *
667  *
668 
669  *
670  *
671  * <a name="ObstacleProblemObstacleProblem"></a>
672  * <h4>ObstacleProblem::ObstacleProblem</h4>
673  *
674 
675  *
676  * To everyone who has taken a look at the first few tutorial programs, the
677  * constructor is completely obvious:
678  *
679  * @code
680  * template <int dim>
681  * ObstacleProblem<dim>::ObstacleProblem()
682  * : fe(1)
683  * , dof_handler(triangulation)
684  * {}
685  *
686  *
687  * @endcode
688  *
689  *
690  * <a name="ObstacleProblemmake_grid"></a>
691  * <h4>ObstacleProblem::make_grid</h4>
692  *
693 
694  *
695  * We solve our obstacle problem on the square @f$[-1,1]\times [-1,1]@f$ in
696  * 2D. This function therefore just sets up one of the simplest possible
697  * meshes.
698  *
699  * @code
700  * template <int dim>
701  * void ObstacleProblem<dim>::make_grid()
702  * {
704  * triangulation.refine_global(7);
705  *
706  * std::cout << "Number of active cells: " << triangulation.n_active_cells()
707  * << std::endl
708  * << "Total number of cells: " << triangulation.n_cells()
709  * << std::endl;
710  * }
711  *
712  *
713  * @endcode
714  *
715  *
716  * <a name="ObstacleProblemsetup_system"></a>
717  * <h4>ObstacleProblem::setup_system</h4>
718  *
719 
720  *
721  * In this first function of note, we set up the degrees of freedom handler,
722  * resize vectors and matrices, and deal with the constraints. Initially,
723  * the constraints are, of course, only given by boundary values, so we
724  * interpolate them towards the top of the function.
725  *
726  * @code
727  * template <int dim>
728  * void ObstacleProblem<dim>::setup_system()
729  * {
730  * dof_handler.distribute_dofs(fe);
731  * active_set.set_size(dof_handler.n_dofs());
732  *
733  * std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
734  * << std::endl
735  * << std::endl;
736  *
738  * 0,
739  * BoundaryValues<dim>(),
740  * constraints);
741  * constraints.close();
742  *
743  * DynamicSparsityPattern dsp(dof_handler.n_dofs());
744  * DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
745  *
746  * system_matrix.reinit(dsp);
747  * complete_system_matrix.reinit(dsp);
748  *
749  * IndexSet solution_index_set = dof_handler.locally_owned_dofs();
750  * solution.reinit(solution_index_set, MPI_COMM_WORLD);
751  * system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
752  * complete_system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
753  * contact_force.reinit(solution_index_set, MPI_COMM_WORLD);
754  *
755  * @endcode
756  *
757  * The only other thing to do here is to compute the factors in the @f$B@f$
758  * matrix which is used to scale the residual. As discussed in the
759  * introduction, we'll use a little trick to make this mass matrix
760  * diagonal, and in the following then first compute all of this as a
761  * matrix and then extract the diagonal elements for later use:
762  *
763  * @code
764  * TrilinosWrappers::SparseMatrix mass_matrix;
765  * mass_matrix.reinit(dsp);
766  * assemble_mass_matrix_diagonal(mass_matrix);
767  * diagonal_of_mass_matrix.reinit(solution_index_set);
768  * for (unsigned int j = 0; j < solution.size(); j++)
769  * diagonal_of_mass_matrix(j) = mass_matrix.diag_element(j);
770  * }
771  *
772  *
773  * @endcode
774  *
775  *
776  * <a name="ObstacleProblemassemble_system"></a>
777  * <h4>ObstacleProblem::assemble_system</h4>
778  *
779 
780  *
781  * This function at once assembles the system matrix and right-hand-side and
782  * applied the constraints (both due to the active set as well as from
783  * boundary values) to our system. Otherwise, it is functionally equivalent
784  * to the corresponding function in, for example, @ref step_4 "step-4".
785  *
786  * @code
787  * template <int dim>
788  * void ObstacleProblem<dim>::assemble_system()
789  * {
790  * std::cout << " Assembling system..." << std::endl;
791  *
792  * system_matrix = 0;
793  * system_rhs = 0;
794  *
795  * const QGauss<dim> quadrature_formula(fe.degree + 1);
796  * RightHandSide<dim> right_hand_side;
797  *
798  * FEValues<dim> fe_values(fe,
799  * quadrature_formula,
800  * update_values | update_gradients |
801  * update_quadrature_points | update_JxW_values);
802  *
803  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
804  * const unsigned int n_q_points = quadrature_formula.size();
805  *
806  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
807  * Vector<double> cell_rhs(dofs_per_cell);
808  *
809  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
810  *
811  * for (const auto &cell : dof_handler.active_cell_iterators())
812  * {
813  * fe_values.reinit(cell);
814  * cell_matrix = 0;
815  * cell_rhs = 0;
816  *
817  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
818  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
819  * {
820  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
821  * cell_matrix(i, j) +=
822  * (fe_values.shape_grad(i, q_point) *
823  * fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
824  *
825  * cell_rhs(i) +=
826  * (fe_values.shape_value(i, q_point) *
827  * right_hand_side.value(fe_values.quadrature_point(q_point)) *
828  * fe_values.JxW(q_point));
829  * }
830  *
831  * cell->get_dof_indices(local_dof_indices);
832  *
833  * constraints.distribute_local_to_global(cell_matrix,
834  * cell_rhs,
835  * local_dof_indices,
836  * system_matrix,
837  * system_rhs,
838  * true);
839  * }
840  * }
841  *
842  *
843  *
844  * @endcode
845  *
846  *
847  * <a name="ObstacleProblemassemble_mass_matrix_diagonal"></a>
848  * <h4>ObstacleProblem::assemble_mass_matrix_diagonal</h4>
849  *
850 
851  *
852  * The next function is used in the computation of the diagonal mass matrix
853  * @f$B@f$ used to scale variables in the active set method. As discussed in the
854  * introduction, we get the mass matrix to be diagonal by choosing the
855  * trapezoidal rule for quadrature. Doing so we don't really need the triple
856  * loop over quadrature points, indices @f$i@f$ and indices @f$j@f$ any more and
857  * can, instead, just use a double loop. The rest of the function is obvious
858  * given what we have discussed in many of the previous tutorial programs.
859  *
860 
861  *
862  * Note that at the time this function is called, the constraints object
863  * only contains boundary value constraints; we therefore do not have to pay
864  * attention in the last copy-local-to-global step to preserve the values of
865  * matrix entries that may later on be constrained by the active set.
866  *
867 
868  *
869  * Note also that the trick with the trapezoidal rule only works if we have
870  * in fact @f$Q_1@f$ elements. For higher order elements, one would need to use
871  * a quadrature formula that has quadrature points at all the support points
872  * of the finite element. Constructing such a quadrature formula isn't
873  * really difficult, but not the point here, and so we simply assert at the
874  * top of the function that our implicit assumption about the finite element
875  * is in fact satisfied.
876  *
877  * @code
878  * template <int dim>
879  * void ObstacleProblem<dim>::assemble_mass_matrix_diagonal(
880  * TrilinosWrappers::SparseMatrix &mass_matrix)
881  * {
882  * Assert(fe.degree == 1, ExcNotImplemented());
883  *
884  * const QTrapez<dim> quadrature_formula;
885  * FEValues<dim> fe_values(fe,
886  * quadrature_formula,
887  * update_values | update_JxW_values);
888  *
889  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
890  * const unsigned int n_q_points = quadrature_formula.size();
891  *
892  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
893  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
894  *
895  * for (const auto &cell : dof_handler.active_cell_iterators())
896  * {
897  * fe_values.reinit(cell);
898  * cell_matrix = 0;
899  *
900  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
901  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
902  * cell_matrix(i, i) +=
903  * (fe_values.shape_value(i, q_point) *
904  * fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
905  *
906  * cell->get_dof_indices(local_dof_indices);
907  *
908  * constraints.distribute_local_to_global(cell_matrix,
909  * local_dof_indices,
910  * mass_matrix);
911  * }
912  * }
913  *
914  *
915  * @endcode
916  *
917  *
918  * <a name="ObstacleProblemupdate_solution_and_constraints"></a>
919  * <h4>ObstacleProblem::update_solution_and_constraints</h4>
920  *
921 
922  *
923  * In a sense, this is the central function of this program. It updates the
924  * active set of constrained degrees of freedom as discussed in the
925  * introduction and computes an AffineConstraints object from it that can then
926  * be used to eliminate constrained degrees of freedom from the solution of
927  * the next iteration. At the same time we set the constrained degrees of
928  * freedom of the solution to the correct value, namely the height of the
929  * obstacle.
930  *
931 
932  *
933  * Fundamentally, the function is rather simple: We have to loop over all
934  * degrees of freedom and check the sign of the function @f$\Lambda^k_i +
935  * c([BU^k]_i - G_i) = \Lambda^k_i + cB_i(U^k_i - [g_h]_i)@f$ because in our
936  * case @f$G_i = B_i[g_h]_i@f$. To this end, we use the formula given in the
937  * introduction by which we can compute the Lagrange multiplier as the
938  * residual of the original linear system (given via the variables
939  * <code>complete_system_matrix</code> and <code>complete_system_rhs</code>.
940  * At the top of this function, we compute this residual using a function
941  * that is part of the matrix classes.
942  *
943  * @code
944  * template <int dim>
945  * void ObstacleProblem<dim>::update_solution_and_constraints()
946  * {
947  * std::cout << " Updating active set..." << std::endl;
948  *
949  * const double penalty_parameter = 100.0;
950  *
951  * TrilinosWrappers::MPI::Vector lambda(
952  * complete_index_set(dof_handler.n_dofs()));
953  * complete_system_matrix.residual(lambda, solution, complete_system_rhs);
954  *
955  * @endcode
956  *
957  * compute contact_force[i] = - lambda[i] * diagonal_of_mass_matrix[i]
958  *
959  * @code
960  * contact_force = lambda;
961  * contact_force.scale(diagonal_of_mass_matrix);
962  * contact_force *= -1;
963  *
964  * @endcode
965  *
966  * The next step is to reset the active set and constraints objects and to
967  * start the loop over all degrees of freedom. This is made slightly more
968  * complicated by the fact that we can't just loop over all elements of
969  * the solution vector since there is no way for us then to find out what
970  * location a DoF is associated with; however, we need this location to
971  * test whether the displacement of a DoF is larger or smaller than the
972  * height of the obstacle at this location.
973  *
974 
975  *
976  * We work around this by looping over all cells and DoFs defined on each
977  * of these cells. We use here that the displacement is described using a
978  * @f$Q_1@f$ function for which degrees of freedom are always located on the
979  * vertices of the cell; thus, we can get the index of each degree of
980  * freedom and its location by asking the vertex for this information. On
981  * the other hand, this clearly wouldn't work for higher order elements,
982  * and so we add an assertion that makes sure that we only deal with
983  * elements for which all degrees of freedom are located in vertices to
984  * avoid tripping ourselves with non-functional code in case someone wants
985  * to play with increasing the polynomial degree of the solution.
986  *
987 
988  *
989  * The price to pay for having to loop over cells rather than DoFs is that
990  * we may encounter some degrees of freedom more than once, namely each
991  * time we visit one of the cells adjacent to a given vertex. We will
992  * therefore have to keep track which vertices we have already touched and
993  * which we haven't so far. We do so by using an array of flags
994  * <code>dof_touched</code>:
995  *
996  * @code
997  * constraints.clear();
998  * active_set.clear();
999  *
1000  * const Obstacle<dim> obstacle;
1001  * std::vector<bool> dof_touched(dof_handler.n_dofs(), false);
1002  *
1003  * for (const auto &cell : dof_handler.active_cell_iterators())
1004  * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1005  * {
1006  * Assert(dof_handler.get_fe().dofs_per_cell ==
1008  * ExcNotImplemented());
1009  *
1010  * const unsigned int dof_index = cell->vertex_dof_index(v, 0);
1011  *
1012  * if (dof_touched[dof_index] == false)
1013  * dof_touched[dof_index] = true;
1014  * else
1015  * continue;
1016  *
1017  * @endcode
1018  *
1019  * Now that we know that we haven't touched this DoF yet, let's get
1020  * the value of the displacement function there as well as the value
1021  * of the obstacle function and use this to decide whether the
1022  * current DoF belongs to the active set. For that we use the
1023  * function given above and in the introduction.
1024  *
1025 
1026  *
1027  * If we decide that the DoF should be part of the active set, we
1028  * add its index to the active set, introduce an inhomogeneous
1029  * equality constraint in the AffineConstraints object, and reset the
1030  * solution value to the height of the obstacle. Finally, the
1031  * residual of the non-contact part of the system serves as an
1032  * additional control (the residual equals the remaining,
1033  * unaccounted forces, and should be zero outside the contact zone),
1034  * so we zero out the components of the residual vector (i.e., the
1035  * Lagrange multiplier lambda) that correspond to the area where the
1036  * body is in contact; at the end of the loop over all cells, the
1037  * residual will therefore only consist of the residual in the
1038  * non-contact zone. We output the norm of this residual along with
1039  * the size of the active set after the loop.
1040  *
1041  * @code
1042  * const double obstacle_value = obstacle.value(cell->vertex(v));
1043  * const double solution_value = solution(dof_index);
1044  *
1045  * if (lambda(dof_index) + penalty_parameter *
1046  * diagonal_of_mass_matrix(dof_index) *
1047  * (solution_value - obstacle_value) <
1048  * 0)
1049  * {
1050  * active_set.add_index(dof_index);
1051  * constraints.add_line(dof_index);
1052  * constraints.set_inhomogeneity(dof_index, obstacle_value);
1053  *
1054  * solution(dof_index) = obstacle_value;
1055  *
1056  * lambda(dof_index) = 0;
1057  * }
1058  * }
1059  * std::cout << " Size of active set: " << active_set.n_elements()
1060  * << std::endl;
1061  *
1062  * std::cout << " Residual of the non-contact part of the system: "
1063  * << lambda.l2_norm() << std::endl;
1064  *
1065  * @endcode
1066  *
1067  * In a final step, we add to the set of constraints on DoFs we have so
1068  * far from the active set those that result from Dirichlet boundary
1069  * values, and close the constraints object:
1070  *
1071  * @code
1073  * 0,
1074  * BoundaryValues<dim>(),
1075  * constraints);
1076  * constraints.close();
1077  * }
1078  *
1079  * @endcode
1080  *
1081  *
1082  * <a name="ObstacleProblemsolve"></a>
1083  * <h4>ObstacleProblem::solve</h4>
1084  *
1085 
1086  *
1087  * There is nothing to say really about the solve function. In the context
1088  * of a Newton method, we are not typically interested in very high accuracy
1089  * (why ask for a highly accurate solution of a linear problem that we know
1090  * only gives us an approximation of the solution of the nonlinear problem),
1091  * and so we use the ReductionControl class that stops iterations when
1092  * either an absolute tolerance is reached (for which we choose @f$10^{-12}@f$)
1093  * or when the residual is reduced by a certain factor (here, @f$10^{-3}@f$).
1094  *
1095  * @code
1096  * template <int dim>
1097  * void ObstacleProblem<dim>::solve()
1098  * {
1099  * std::cout << " Solving system..." << std::endl;
1100  *
1101  * ReductionControl reduction_control(100, 1e-12, 1e-3);
1102  * SolverCG<TrilinosWrappers::MPI::Vector> solver(reduction_control);
1103  * TrilinosWrappers::PreconditionAMG precondition;
1104  * precondition.initialize(system_matrix);
1105  *
1106  * solver.solve(system_matrix, solution, system_rhs, precondition);
1107  * constraints.distribute(solution);
1108  *
1109  * std::cout << " Error: " << reduction_control.initial_value() << " -> "
1110  * << reduction_control.last_value() << " in "
1111  * << reduction_control.last_step() << " CG iterations."
1112  * << std::endl;
1113  * }
1114  *
1115  *
1116  * @endcode
1117  *
1118  *
1119  * <a name="ObstacleProblemoutput_results"></a>
1120  * <h4>ObstacleProblem::output_results</h4>
1121  *
1122 
1123  *
1124  * We use the vtk-format for the output. The file contains the displacement
1125  * and a numerical representation of the active set.
1126  *
1127  * @code
1128  * template <int dim>
1129  * void ObstacleProblem<dim>::output_results(const unsigned int iteration) const
1130  * {
1131  * std::cout << " Writing graphical output..." << std::endl;
1132  *
1133  * TrilinosWrappers::MPI::Vector active_set_vector(
1134  * dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
1135  * for (const auto index : active_set)
1136  * active_set_vector[index] = 1.;
1137  *
1138  * DataOut<dim> data_out;
1139  *
1140  * data_out.attach_dof_handler(dof_handler);
1141  * data_out.add_data_vector(solution, "displacement");
1142  * data_out.add_data_vector(active_set_vector, "active_set");
1143  * data_out.add_data_vector(contact_force, "lambda");
1144  *
1145  * data_out.build_patches();
1146  *
1147  * std::ofstream output_vtk("output_" +
1148  * Utilities::int_to_string(iteration, 3) + ".vtk");
1149  * data_out.write_vtk(output_vtk);
1150  * }
1151  *
1152  *
1153  *
1154  * @endcode
1155  *
1156  *
1157  * <a name="ObstacleProblemrun"></a>
1158  * <h4>ObstacleProblem::run</h4>
1159  *
1160 
1161  *
1162  * This is the function which has the top-level control over everything. It
1163  * is not very long, and in fact rather straightforward: in every iteration
1164  * of the active set method, we assemble the linear system, solve it, update
1165  * the active set and project the solution back to the feasible set, and
1166  * then output the results. The iteration is terminated whenever the active
1167  * set has not changed in the previous iteration.
1168  *
1169 
1170  *
1171  * The only trickier part is that we have to save the linear system (i.e.,
1172  * the matrix and right hand side) after assembling it in the first
1173  * iteration. The reason is that this is the only step where we can access
1174  * the linear system as built without any of the contact constraints
1175  * active. We need this to compute the residual of the solution at other
1176  * iterations, but in other iterations that linear system we form has the
1177  * rows and columns that correspond to constrained degrees of freedom
1178  * eliminated, and so we can no longer access the full residual of the
1179  * original equation.
1180  *
1181  * @code
1182  * template <int dim>
1183  * void ObstacleProblem<dim>::run()
1184  * {
1185  * make_grid();
1186  * setup_system();
1187  *
1188  * IndexSet active_set_old(active_set);
1189  * for (unsigned int iteration = 0; iteration <= solution.size(); ++iteration)
1190  * {
1191  * std::cout << "Newton iteration " << iteration << std::endl;
1192  *
1193  * assemble_system();
1194  *
1195  * if (iteration == 0)
1196  * {
1197  * complete_system_matrix.copy_from(system_matrix);
1198  * complete_system_rhs = system_rhs;
1199  * }
1200  *
1201  * solve();
1202  * update_solution_and_constraints();
1203  * output_results(iteration);
1204  *
1205  * if (active_set == active_set_old)
1206  * break;
1207  *
1208  * active_set_old = active_set;
1209  *
1210  * std::cout << std::endl;
1211  * }
1212  * }
1213  * } // namespace Step41
1214  *
1215  *
1216  * @endcode
1217  *
1218  *
1219  * <a name="Thecodemaincodefunction"></a>
1220  * <h3>The <code>main</code> function</h3>
1221  *
1222 
1223  *
1224  * And this is the main function. It follows the pattern of all other main
1225  * functions. The call to initialize MPI exists because the Trilinos library
1226  * upon which we build our linear solvers in this program requires it.
1227  *
1228  * @code
1229  * int main(int argc, char *argv[])
1230  * {
1231  * try
1232  * {
1233  * using namespace dealii;
1234  * using namespace Step41;
1235  *
1236  * Utilities::MPI::MPI_InitFinalize mpi_initialization(
1237  * argc, argv, numbers::invalid_unsigned_int);
1238  *
1239  * @endcode
1240  *
1241  * This program can only be run in serial. Otherwise, throw an exception.
1242  *
1243  * @code
1244  * AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
1245  * ExcMessage(
1246  * "This program can only be run in serial, use ./step-41"));
1247  *
1248  * ObstacleProblem<2> obstacle_problem;
1249  * obstacle_problem.run();
1250  * }
1251  * catch (std::exception &exc)
1252  * {
1253  * std::cerr << std::endl
1254  * << std::endl
1255  * << "----------------------------------------------------"
1256  * << std::endl;
1257  * std::cerr << "Exception on processing: " << std::endl
1258  * << exc.what() << std::endl
1259  * << "Aborting!" << std::endl
1260  * << "----------------------------------------------------"
1261  * << std::endl;
1262  *
1263  * return 1;
1264  * }
1265  * catch (...)
1266  * {
1267  * std::cerr << std::endl
1268  * << std::endl
1269  * << "----------------------------------------------------"
1270  * << std::endl;
1271  * std::cerr << "Unknown exception!" << std::endl
1272  * << "Aborting!" << std::endl
1273  * << "----------------------------------------------------"
1274  * << std::endl;
1275  * return 1;
1276  * }
1277  *
1278  * return 0;
1279  * }
1280  * @endcode
1281 <a name="Results"></a><h1>Results</h1>
1282 
1283 
1284 Running the program produces output like this:
1285 @code
1286 Number of active cells: 16384
1287 Total number of cells: 21845
1288 Number of degrees of freedom: 16641
1289 
1290 Newton iteration 0
1291  Assembling system...
1292  Solving system...
1293  Error: 0.310059 -> 5.16619e-05 in 5 CG iterations.
1294  Updating active set...
1295  Size of active set: 13164
1296  Residual of the non-contact part of the system: 1.61863e-05
1297  Writing graphical output...
1298 
1299 Newton iteration 1
1300  Assembling system...
1301  Solving system...
1302  Error: 1.11987 -> 0.00109377 in 6 CG iterations.
1303  Updating active set...
1304  Size of active set: 12363
1305  Residual of the non-contact part of the system: 3.9373
1306  Writing graphical output...
1307 
1308 ...
1309 
1310 Newton iteration 17
1311  Assembling system...
1312  Solving system...
1313  Error: 0.00713308 -> 2.29249e-06 in 4 CG iterations.
1314  Updating active set...
1315  Size of active set: 5399
1316  Residual of the non-contact part of the system: 0.000957525
1317  Writing graphical output...
1318 
1319 Newton iteration 18
1320  Assembling system...
1321  Solving system...
1322  Error: 0.000957525 -> 2.8033e-07 in 4 CG iterations.
1323  Updating active set...
1324  Size of active set: 5399
1325  Residual of the non-contact part of the system: 2.8033e-07
1326  Writing graphical output...
1327 @endcode
1328 
1329 The iterations end once the active set doesn't change any more (it has
1330 5,399 constrained degrees of freedom at that point). The algebraic
1331 precondition is apparently working nicely since we only need 4-6 CG
1332 iterations to solve the linear system (although this also has a lot to
1333 do with the fact that we are not asking for very high accuracy of the
1334 linear solver).
1335 
1336 More revealing is to look at a sequence of graphical output files
1337 (every third step is shown, with the number of the iteration in the
1338 leftmost column):
1339 
1340 <table align="center">
1341  <tr>
1342  <td valign="top">
1343  0 &nbsp;
1344  </td>
1345  <td valign="top">
1346  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
1347  </td>
1348  <td valign="top">
1349  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
1350  </td>
1351  <td valign="top">
1352  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
1353  </td>
1354  </tr>
1355  <tr>
1356  <td valign="top">
1357  3 &nbsp;
1358  </td>
1359  <td valign="top">
1360  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
1361  </td>
1362  <td valign="top">
1363  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
1364  </td>
1365  <td valign="top">
1366  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
1367  </td>
1368  </tr>
1369  <tr>
1370  <td valign="top">
1371  6 &nbsp;
1372  </td>
1373  <td valign="top">
1374  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
1375  </td>
1376  <td valign="top">
1377  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
1378  </td>
1379  <td valign="top">
1380  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
1381  </td>
1382  </tr>
1383  <tr>
1384  <td valign="top">
1385  9 &nbsp;
1386  </td>
1387  <td valign="top">
1388  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
1389  </td>
1390  <td valign="top">
1391  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
1392  </td>
1393  <td valign="top">
1394  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
1395  </td>
1396  </tr>
1397  <tr>
1398  <td valign="top">
1399  12 &nbsp;
1400  </td>
1401  <td valign="top">
1402  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
1403  </td>
1404  <td valign="top">
1405  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
1406  </td>
1407  <td valign="top">
1408  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
1409  </td>
1410  </tr>
1411  <tr>
1412  <td valign="top">
1413  15 &nbsp;
1414  </td>
1415  <td valign="top">
1416  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
1417  </td>
1418  <td valign="top">
1419  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
1420  </td>
1421  <td valign="top">
1422  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
1423  </td>
1424  </tr>
1425  <tr>
1426  <td valign="top">
1427  18 &nbsp;
1428  </td>
1429  <td valign="top">
1430  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
1431  </td>
1432  <td valign="top">
1433  <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
1434  </td>
1435  <td valign="top">
1436  <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
1437  </td>
1438  </tr>
1439 </table>
1440 
1441 The pictures show that in the first step, the solution (which has been
1442 computed without any of the constraints active) bends through so much
1443 that pretty much every interior point has to be bounced back to the
1444 stairstep function, producing a discontinuous solution. Over the
1445 course of the active set iterations, this unphysical membrane shape is
1446 smoothed out, the contact with the lower-most stair step disappears,
1447 and the solution stabilizes.
1448 
1449 In addition to this, the program also outputs the values of the
1450 Lagrange multipliers. Remember that these are the contact forces and
1451 so should only be positive on the contact set, and zero outside. If,
1452 on the other hand, a Lagrange multiplier is negative in the active
1453 set, then this degree of freedom must be removed from the active
1454 set. The following pictures show the multipliers in iterations 1, 9
1455 and 18, where we use red and browns to indicate positive values, and
1456 blue for negative values.
1457 
1458 <table align="center">
1459  <tr>
1460  <td valign="top">
1461  <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
1462  </td>
1463  <td valign="top">
1464  <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
1465  </td>
1466  <td valign="top">
1467  <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
1468  </td>
1469  </tr>
1470  <tr>
1471  <td align="center">
1472  Iteration 1
1473  </td>
1474  <td align="center">
1475  Iteration 9
1476  </td>
1477  <td align="center">
1478  Iteration 18
1479  </td>
1480  </tr>
1481 </table>
1482 
1483 It is easy to see that the positive values converge nicely to moderate
1484 values in the interior of the contact set and large upward forces at
1485 the edges of the steps, as one would expect (to support the large
1486 curvature of the membrane there); at the fringes of the active set,
1487 multipliers are initially negative, causing the set to shrink until,
1488 in iteration 18, there are no more negative multipliers and the
1489 algorithm has converged.
1490 
1491 
1492 
1493 <a name="extensions"></a>
1494 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1495 
1496 
1497 As with any of the programs of this tutorial, there are a number of
1498 obvious possibilities for extensions and experiments. The first one is
1499 clear: introduce adaptivity. Contact problems are prime candidates for
1500 adaptive meshes because the solution has lines along which it is less
1501 regular (the places where contact is established between membrane and
1502 obstacle) and other areas where the solution is very smooth (or, in
1503 the present context, constant wherever it is in contact with the
1504 obstacle). Adding this to the current program should not pose too many
1505 difficulties, but it is not trivial to find a good error estimator for
1506 that purpose.
1507 
1508 A more challenging task would be an extension to 3d. The problem here
1509 is not so much to simply make everything run in 3d. Rather, it is that
1510 when a 3d body is deformed and gets into contact with an obstacle,
1511 then the obstacle does not act as a constraining body force within the
1512 domain as is the case here. Rather, the contact force only acts on the
1513 boundary of the object. The inequality then is not in the differential
1514 equation but in fact in the (Neumann-type) boundary conditions, though
1515 this leads to a similar kind of variational
1516 inequality. Mathematically, this means that the Lagrange multiplier
1517 only lives on the surface, though it can of course be extended by zero
1518 into the domain if that is convenient. As in the current program, one
1519 does not need to form and store this Lagrange multiplier explicitly.
1520 
1521 A further interesting problem for the 3d case is to consider contact problems
1522 with friction. In almost every mechanical process friction has a big influence.
1523 For the modelling we have to take into account tangential stresses at the contact
1524 surface. Also we have to observe that friction adds another nonlinearity to
1525 our problem.
1526 
1527 Another nontrivial modification is to implement a more complex constitutive
1528 law like nonlinear elasticity or elasto-plastic material behavior.
1529 The difficulty here is to handle the additional nonlinearity arising
1530 through the nonlinear constitutive law.
1531  *
1532  *
1533 <a name="PlainProg"></a>
1534 <h1> The plain program</h1>
1535 @include "step-41.cc"
1536 */
AdaptationStrategies::Refinement::preserve
std::vector< value_type > preserve(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
IndexSet
Definition: index_set.h:74
TrilinosWrappers::MPI::Vector
Definition: trilinos_vector.h:400
SolverCG
Definition: solver_cg.h:98
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
dealii
Definition: namespace_dealii.h:25
Differentiation::SD::OptimizerType::lambda
@ lambda
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
DataOutBase::vtk
@ vtk
Definition: data_out_base.h:1605
GeometryInfo
Definition: geometry_info.h:1224
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >(0)), const bool project_to_boundary_first=false)
ReductionControl
Definition: solver_control.h:427
DoFTools::always
@ always
Definition: dof_tools.h:236
DataOutInterface::write_vtk
void write_vtk(std::ostream &out) const
Definition: data_out_base.cc:6853
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
TrilinosWrappers::PreconditionAMG
Definition: trilinos_precondition.h:1361
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Tensor::clear
constexpr void clear()
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
value
static const bool value
Definition: dof_tools_constraints.cc:433
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
AffineConstraints
Definition: affine_constraints.h:180
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
Point< dim >
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
TrilinosWrappers::PreconditionAMG::initialize
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: trilinos_precondition_ml.cc:222
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
int
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:828
DataOut_DoFData::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090