Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10:02+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,
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="step_41-ImplementationofthecodeObstacleProblemcodeclass"></a>
660 * <h3>Implementation of the <code>ObstacleProblem</code> class</h3>
661 *
662
663 *
664 *
665
666 *
667 *
668 * <a name="step_41-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="step_41-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="step_41-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="step_41-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="step_41-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="step_41-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_constraint(dof_index, {}, obstacle_value);
1048 *  
1049 *   solution(dof_index) = obstacle_value;
1050 *  
1051 *   lambda(dof_index) = 0;
1052 *   }
1053 *   }
1054 *   std::cout << " Size of active set: " << active_set.n_elements()
1055 *   << std::endl;
1056 *  
1057 *   std::cout << " Residual of the non-contact part of the system: "
1058 *   << lambda.l2_norm() << std::endl;
1059 *  
1060 * @endcode
1061 *
1062 * In a final step, we add to the set of constraints on DoFs we have so
1063 * far from the active set those that result from Dirichlet boundary
1064 * values, and close the constraints object:
1065 *
1066 * @code
1068 *   0,
1069 *   BoundaryValues<dim>(),
1070 *   constraints);
1071 *   constraints.close();
1072 *   }
1073 *  
1074 * @endcode
1075 *
1076 *
1077 * <a name="step_41-ObstacleProblemsolve"></a>
1078 * <h4>ObstacleProblem::solve</h4>
1079 *
1080
1081 *
1082 * There is nothing to say really about the solve function. In the context
1083 * of a Newton method, we are not typically interested in very high accuracy
1084 * (why ask for a highly accurate solution of a linear problem that we know
1085 * only gives us an approximation of the solution of the nonlinear problem),
1086 * and so we use the ReductionControl class that stops iterations when
1087 * either an absolute tolerance is reached (for which we choose @f$10^{-12}@f$)
1088 * or when the residual is reduced by a certain factor (here, @f$10^{-3}@f$).
1089 *
1090 * @code
1091 *   template <int dim>
1092 *   void ObstacleProblem<dim>::solve()
1093 *   {
1094 *   std::cout << " Solving system..." << std::endl;
1095 *  
1096 *   ReductionControl reduction_control(100, 1e-12, 1e-3);
1097 *   SolverCG<TrilinosWrappers::MPI::Vector> solver(reduction_control);
1098 *   TrilinosWrappers::PreconditionAMG precondition;
1099 *   precondition.initialize(system_matrix);
1100 *  
1101 *   solver.solve(system_matrix, solution, system_rhs, precondition);
1102 *   constraints.distribute(solution);
1103 *  
1104 *   std::cout << " Error: " << reduction_control.initial_value() << " -> "
1105 *   << reduction_control.last_value() << " in "
1106 *   << reduction_control.last_step() << " CG iterations."
1107 *   << std::endl;
1108 *   }
1109 *  
1110 *  
1111 * @endcode
1112 *
1113 *
1114 * <a name="step_41-ObstacleProblemoutput_results"></a>
1115 * <h4>ObstacleProblem::output_results</h4>
1116 *
1117
1118 *
1119 * We use the vtk-format for the output. The file contains the displacement
1120 * and a numerical representation of the active set.
1121 *
1122 * @code
1123 *   template <int dim>
1124 *   void ObstacleProblem<dim>::output_results(const unsigned int iteration) const
1125 *   {
1126 *   std::cout << " Writing graphical output..." << std::endl;
1127 *  
1128 *   TrilinosWrappers::MPI::Vector active_set_vector(
1129 *   dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
1130 *   for (const auto index : active_set)
1131 *   active_set_vector[index] = 1.;
1132 *  
1133 *   DataOut<dim> data_out;
1134 *  
1135 *   data_out.attach_dof_handler(dof_handler);
1136 *   data_out.add_data_vector(solution, "displacement");
1137 *   data_out.add_data_vector(active_set_vector, "active_set");
1138 *   data_out.add_data_vector(contact_force, "lambda");
1139 *  
1140 *   data_out.build_patches();
1141 *  
1142 *   std::ofstream output_vtk("output_" +
1143 *   Utilities::int_to_string(iteration, 3) + ".vtk");
1144 *   data_out.write_vtk(output_vtk);
1145 *   }
1146 *  
1147 *  
1148 *  
1149 * @endcode
1150 *
1151 *
1152 * <a name="step_41-ObstacleProblemrun"></a>
1153 * <h4>ObstacleProblem::run</h4>
1154 *
1155
1156 *
1157 * This is the function which has the top-level control over everything. It
1158 * is not very long, and in fact rather straightforward: in every iteration
1159 * of the active set method, we assemble the linear system, solve it, update
1160 * the active set and project the solution back to the feasible set, and
1161 * then output the results. The iteration is terminated whenever the active
1162 * set has not changed in the previous iteration.
1163 *
1164
1165 *
1166 * The only trickier part is that we have to save the linear system (i.e.,
1167 * the matrix and right hand side) after assembling it in the first
1168 * iteration. The reason is that this is the only step where we can access
1169 * the linear system as built without any of the contact constraints
1170 * active. We need this to compute the residual of the solution at other
1171 * iterations, but in other iterations that linear system we form has the
1172 * rows and columns that correspond to constrained degrees of freedom
1173 * eliminated, and so we can no longer access the full residual of the
1174 * original equation.
1175 *
1176 * @code
1177 *   template <int dim>
1178 *   void ObstacleProblem<dim>::run()
1179 *   {
1180 *   make_grid();
1181 *   setup_system();
1182 *  
1183 *   IndexSet active_set_old(active_set);
1184 *   for (unsigned int iteration = 0; iteration <= solution.size(); ++iteration)
1185 *   {
1186 *   std::cout << "Newton iteration " << iteration << std::endl;
1187 *  
1188 *   assemble_system();
1189 *  
1190 *   if (iteration == 0)
1191 *   {
1192 *   complete_system_matrix.copy_from(system_matrix);
1193 *   complete_system_rhs = system_rhs;
1194 *   }
1195 *  
1196 *   solve();
1197 *   update_solution_and_constraints();
1198 *   output_results(iteration);
1199 *  
1200 *   if (active_set == active_set_old)
1201 *   break;
1202 *  
1203 *   active_set_old = active_set;
1204 *  
1205 *   std::cout << std::endl;
1206 *   }
1207 *   }
1208 *   } // namespace Step41
1209 *  
1210 *  
1211 * @endcode
1212 *
1213 *
1214 * <a name="step_41-Thecodemaincodefunction"></a>
1215 * <h3>The <code>main</code> function</h3>
1216 *
1217
1218 *
1219 * And this is the main function. It follows the pattern of all other main
1220 * functions. The call to initialize MPI exists because the Trilinos library
1221 * upon which we build our linear solvers in this program requires it.
1222 *
1223 * @code
1224 *   int main(int argc, char *argv[])
1225 *   {
1226 *   try
1227 *   {
1228 *   using namespace dealii;
1229 *   using namespace Step41;
1230 *  
1231 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(
1232 *   argc, argv, numbers::invalid_unsigned_int);
1233 *  
1234 * @endcode
1235 *
1236 * This program can only be run in serial. Otherwise, throw an exception.
1237 *
1238 * @code
1239 *   AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
1240 *   ExcMessage(
1241 *   "This program can only be run in serial, use ./step-41"));
1242 *  
1243 *   ObstacleProblem<2> obstacle_problem;
1244 *   obstacle_problem.run();
1245 *   }
1246 *   catch (std::exception &exc)
1247 *   {
1248 *   std::cerr << std::endl
1249 *   << std::endl
1250 *   << "----------------------------------------------------"
1251 *   << std::endl;
1252 *   std::cerr << "Exception on processing: " << std::endl
1253 *   << exc.what() << std::endl
1254 *   << "Aborting!" << std::endl
1255 *   << "----------------------------------------------------"
1256 *   << std::endl;
1257 *  
1258 *   return 1;
1259 *   }
1260 *   catch (...)
1261 *   {
1262 *   std::cerr << std::endl
1263 *   << std::endl
1264 *   << "----------------------------------------------------"
1265 *   << std::endl;
1266 *   std::cerr << "Unknown exception!" << std::endl
1267 *   << "Aborting!" << std::endl
1268 *   << "----------------------------------------------------"
1269 *   << std::endl;
1270 *   return 1;
1271 *   }
1272 *  
1273 *   return 0;
1274 *   }
1275 * @endcode
1276<a name="step_41-Results"></a><h1>Results</h1>
1277
1278
1279Running the program produces output like this:
1280@code
1281Number of active cells: 16384
1282Total number of cells: 21845
1283Number of degrees of freedom: 16641
1284
1285Newton iteration 0
1286 Assembling system...
1287 Solving system...
1288 Error: 0.310059 -> 5.16619e-05 in 5 CG iterations.
1289 Updating active set...
1290 Size of active set: 13164
1291 Residual of the non-contact part of the system: 1.61863e-05
1292 Writing graphical output...
1293
1294Newton iteration 1
1295 Assembling system...
1296 Solving system...
1297 Error: 1.11987 -> 0.00109377 in 6 CG iterations.
1298 Updating active set...
1299 Size of active set: 12363
1300 Residual of the non-contact part of the system: 3.9373
1301 Writing graphical output...
1302
1303...
1304
1305Newton iteration 17
1306 Assembling system...
1307 Solving system...
1308 Error: 0.00713308 -> 2.29249e-06 in 4 CG iterations.
1309 Updating active set...
1310 Size of active set: 5399
1311 Residual of the non-contact part of the system: 0.000957525
1312 Writing graphical output...
1313
1314Newton iteration 18
1315 Assembling system...
1316 Solving system...
1317 Error: 0.000957525 -> 2.8033e-07 in 4 CG iterations.
1318 Updating active set...
1319 Size of active set: 5399
1320 Residual of the non-contact part of the system: 2.8033e-07
1321 Writing graphical output...
1322@endcode
1323
1324The iterations end once the active set doesn't change any more (it has
13255,399 constrained degrees of freedom at that point). The algebraic
1326precondition is apparently working nicely since we only need 4-6 CG
1327iterations to solve the linear system (although this also has a lot to
1328do with the fact that we are not asking for very high accuracy of the
1329linear solver).
1330
1331More revealing is to look at a sequence of graphical output files
1332(every third step is shown, with the number of the iteration in the
1333leftmost column):
1334
1335<table align="center">
1336 <tr>
1337 <td valign="top">
1338 0 &nbsp;
1339 </td>
1340 <td valign="top">
1341 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
1342 </td>
1343 <td valign="top">
1344 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
1345 </td>
1346 <td valign="top">
1347 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
1348 </td>
1349 </tr>
1350 <tr>
1351 <td valign="top">
1352 3 &nbsp;
1353 </td>
1354 <td valign="top">
1355 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
1356 </td>
1357 <td valign="top">
1358 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
1359 </td>
1360 <td valign="top">
1361 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
1362 </td>
1363 </tr>
1364 <tr>
1365 <td valign="top">
1366 6 &nbsp;
1367 </td>
1368 <td valign="top">
1369 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
1370 </td>
1371 <td valign="top">
1372 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
1373 </td>
1374 <td valign="top">
1375 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
1376 </td>
1377 </tr>
1378 <tr>
1379 <td valign="top">
1380 9 &nbsp;
1381 </td>
1382 <td valign="top">
1383 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
1384 </td>
1385 <td valign="top">
1386 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
1387 </td>
1388 <td valign="top">
1389 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
1390 </td>
1391 </tr>
1392 <tr>
1393 <td valign="top">
1394 12 &nbsp;
1395 </td>
1396 <td valign="top">
1397 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
1398 </td>
1399 <td valign="top">
1400 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
1401 </td>
1402 <td valign="top">
1403 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
1404 </td>
1405 </tr>
1406 <tr>
1407 <td valign="top">
1408 15 &nbsp;
1409 </td>
1410 <td valign="top">
1411 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
1412 </td>
1413 <td valign="top">
1414 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
1415 </td>
1416 <td valign="top">
1417 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
1418 </td>
1419 </tr>
1420 <tr>
1421 <td valign="top">
1422 18 &nbsp;
1423 </td>
1424 <td valign="top">
1425 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
1426 </td>
1427 <td valign="top">
1428 <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
1429 </td>
1430 <td valign="top">
1431 <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
1432 </td>
1433 </tr>
1434</table>
1435
1436The pictures show that in the first step, the solution (which has been
1437computed without any of the constraints active) bends through so much
1438that pretty much every interior point has to be bounced back to the
1439stairstep function, producing a discontinuous solution. Over the
1440course of the active set iterations, this unphysical membrane shape is
1441smoothed out, the contact with the lower-most stair step disappears,
1442and the solution stabilizes.
1443
1444In addition to this, the program also outputs the values of the
1445Lagrange multipliers. Remember that these are the contact forces and
1446so should only be positive on the contact set, and zero outside. If,
1447on the other hand, a Lagrange multiplier is negative in the active
1448set, then this degree of freedom must be removed from the active
1449set. The following pictures show the multipliers in iterations 1, 9
1450and 18, where we use red and browns to indicate positive values, and
1451blue for negative values.
1452
1453<table align="center">
1454 <tr>
1455 <td valign="top">
1456 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
1457 </td>
1458 <td valign="top">
1459 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
1460 </td>
1461 <td valign="top">
1462 <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
1463 </td>
1464 </tr>
1465 <tr>
1466 <td align="center">
1467 Iteration 1
1468 </td>
1469 <td align="center">
1470 Iteration 9
1471 </td>
1472 <td align="center">
1473 Iteration 18
1474 </td>
1475 </tr>
1476</table>
1477
1478It is easy to see that the positive values converge nicely to moderate
1479values in the interior of the contact set and large upward forces at
1480the edges of the steps, as one would expect (to support the large
1481curvature of the membrane there); at the fringes of the active set,
1482multipliers are initially negative, causing the set to shrink until,
1483in iteration 18, there are no more negative multipliers and the
1484algorithm has converged.
1485
1486
1487
1488<a name="step-41-extensions"></a>
1489<a name="step_41-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1490
1491
1492As with any of the programs of this tutorial, there are a number of
1493obvious possibilities for extensions and experiments. The first one is
1494clear: introduce adaptivity. Contact problems are prime candidates for
1495adaptive meshes because the solution has lines along which it is less
1496regular (the places where contact is established between membrane and
1497obstacle) and other areas where the solution is very smooth (or, in
1498the present context, constant wherever it is in contact with the
1499obstacle). Adding this to the current program should not pose too many
1500difficulties, but it is not trivial to find a good error estimator for
1501that purpose.
1502
1503A more challenging task would be an extension to 3d. The problem here
1504is not so much to simply make everything run in 3d. Rather, it is that
1505when a 3d body is deformed and gets into contact with an obstacle,
1506then the obstacle does not act as a constraining body force within the
1507domain as is the case here. Rather, the contact force only acts on the
1508boundary of the object. The inequality then is not in the differential
1509equation but in fact in the (Neumann-type) boundary conditions, though
1510this leads to a similar kind of variational
1511inequality. Mathematically, this means that the Lagrange multiplier
1512only lives on the surface, though it can of course be extended by zero
1513into the domain if that is convenient. As in the current program, one
1514does not need to form and store this Lagrange multiplier explicitly.
1515
1516A further interesting problem for the 3d case is to consider contact problems
1517with friction. In almost every mechanical process friction has a big influence.
1518For the modelling we have to take into account tangential stresses at the contact
1519surface. Also we have to observe that friction adds another nonlinearity to
1520our problem.
1521
1522Another nontrivial modification is to implement a more complex constitutive
1523law like nonlinear elasticity or elasto-plastic material behavior.
1524The difficulty here is to handle the additional nonlinearity arising
1525through the nonlinear constitutive law.
1526 *
1527 *
1528<a name="step_41-PlainProg"></a>
1529<h1> The plain program</h1>
1530@include "step-41.cc"
1531*/
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]
__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(IteratorType begin, std_cxx20::type_identity_t< IteratorType > 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, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
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)
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
const std::vector< bool > & used
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
@ 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)
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 >(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:70
static const unsigned int invalid_unsigned_int
Definition types.h:220
const InputIterator OutputIterator out
Definition parallel.h:167
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation