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