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