Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-15.h
Go to the documentation of this file.
1 ) const
533  * {
534  * return std::sin(2 * numbers::PI * (p[0] + p[1]));
535  * }
536  *
537  * @endcode
538  *
539  *
540  * <a name="ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
541  * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
542  *
543 
544  *
545  *
546  * <a name="MinimalSurfaceProblemMinimalSurfaceProblem"></a>
547  * <h4>MinimalSurfaceProblem::MinimalSurfaceProblem</h4>
548  *
549 
550  *
551  * The constructor and destructor of the class are the same as in the first
552  * few tutorials.
553  *
554 
555  *
556  *
557  * @code
558  * template <int dim>
559  * MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
560  * : dof_handler(triangulation)
561  * , fe(2)
562  * {}
563  *
564  *
565  *
566  * template <int dim>
567  * MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem()
568  * {
569  * dof_handler.clear();
570  * }
571  *
572  * @endcode
573  *
574  *
575  * <a name="MinimalSurfaceProblemsetup_system"></a>
576  * <h4>MinimalSurfaceProblem::setup_system</h4>
577  *
578 
579  *
580  * As always in the setup-system function, we setup the variables of the
581  * finite element method. There are same differences to @ref step_6 "step-6", because
582  * there we start solving the PDE from scratch in every refinement cycle
583  * whereas here we need to take the solution from the previous mesh onto the
584  * current mesh. Consequently, we can't just reset solution vectors. The
585  * argument passed to this function thus indicates whether we can
586  * distributed degrees of freedom (plus compute constraints) and set the
587  * solution vector to zero or whether this has happened elsewhere already
588  * (specifically, in <code>refine_mesh()</code>).
589  *
590 
591  *
592  *
593  * @code
594  * template <int dim>
595  * void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
596  * {
597  * if (initial_step)
598  * {
599  * dof_handler.distribute_dofs(fe);
600  * present_solution.reinit(dof_handler.n_dofs());
601  *
602  * hanging_node_constraints.clear();
603  * DoFTools::make_hanging_node_constraints(dof_handler,
604  * hanging_node_constraints);
605  * hanging_node_constraints.close();
606  * }
607  *
608  *
609  * @endcode
610  *
611  * The remaining parts of the function are the same as in @ref step_6 "step-6".
612  *
613 
614  *
615  *
616  * @code
617  * newton_update.reinit(dof_handler.n_dofs());
618  * system_rhs.reinit(dof_handler.n_dofs());
619  *
620  * DynamicSparsityPattern dsp(dof_handler.n_dofs());
621  * DoFTools::make_sparsity_pattern(dof_handler, dsp);
622  *
623  * hanging_node_constraints.condense(dsp);
624  *
625  * sparsity_pattern.copy_from(dsp);
626  * system_matrix.reinit(sparsity_pattern);
627  * }
628  *
629  * @endcode
630  *
631  *
632  * <a name="MinimalSurfaceProblemassemble_system"></a>
633  * <h4>MinimalSurfaceProblem::assemble_system</h4>
634  *
635 
636  *
637  * This function does the same as in the previous tutorials except that now,
638  * of course, the matrix and right hand side functions depend on the
639  * previous iteration's solution. As discussed in the introduction, we need
640  * to use zero boundary values for the Newton updates; we compute them at
641  * the end of this function.
642  *
643 
644  *
645  * The top of the function contains the usual boilerplate code, setting up
646  * the objects that allow us to evaluate shape functions at quadrature
647  * points and temporary storage locations for the local matrices and
648  * vectors, as well as for the gradients of the previous solution at the
649  * quadrature points. We then start the loop over all cells:
650  *
651  * @code
652  * template <int dim>
653  * void MinimalSurfaceProblem<dim>::assemble_system()
654  * {
655  * const QGauss<dim> quadrature_formula(fe.degree + 1);
656  *
657  * system_matrix = 0;
658  * system_rhs = 0;
659  *
660  * FEValues<dim> fe_values(fe,
661  * quadrature_formula,
664  *
665  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
666  * const unsigned int n_q_points = quadrature_formula.size();
667  *
668  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
669  * Vector<double> cell_rhs(dofs_per_cell);
670  *
671  * std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
672  *
673  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
674  *
675  * for (const auto &cell : dof_handler.active_cell_iterators())
676  * {
677  * cell_matrix = 0;
678  * cell_rhs = 0;
679  *
680  * fe_values.reinit(cell);
681  *
682  * @endcode
683  *
684  * For the assembly of the linear system, we have to obtain the values
685  * of the previous solution's gradients at the quadrature
686  * points. There is a standard way of doing this: the
687  * FEValues::get_function_gradients function takes a vector that
688  * represents a finite element field defined on a DoFHandler, and
689  * evaluates the gradients of this field at the quadrature points of the
690  * cell with which the FEValues object has last been reinitialized.
691  * The values of the gradients at all quadrature points are then written
692  * into the second argument:
693  *
694  * @code
695  * fe_values.get_function_gradients(present_solution,
696  * old_solution_gradients);
697  *
698  * @endcode
699  *
700  * With this, we can then do the integration loop over all quadrature
701  * points and shape functions. Having just computed the gradients of
702  * the old solution in the quadrature points, we are able to compute
703  * the coefficients @f$a_{n}@f$ in these points. The assembly of the
704  * system itself then looks similar to what we always do with the
705  * exception of the nonlinear terms, as does copying the results from
706  * the local objects into the global ones:
707  *
708  * @code
709  * for (unsigned int q = 0; q < n_q_points; ++q)
710  * {
711  * const double coeff =
712  * 1.0 / std::sqrt(1 + old_solution_gradients[q] *
713  * old_solution_gradients[q]);
714  *
715  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
716  * {
717  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
718  * cell_matrix(i, j) +=
719  * (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
720  * * coeff // * a_n
721  * * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
722  * - // -
723  * (fe_values.shape_grad(i, q) // (\nabla \phi_i
724  * * coeff * coeff * coeff // * a_n^3
725  * * (fe_values.shape_grad(j, q) // * (\nabla \phi_j
726  * * old_solution_gradients[q]) // * \nabla u_n)
727  * * old_solution_gradients[q])) // * \nabla u_n)))
728  * * fe_values.JxW(q)); // * dx
729  *
730  * cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
731  * * coeff // * a_n
732  * * old_solution_gradients[q] // * u_n
733  * * fe_values.JxW(q)); // * dx
734  * }
735  * }
736  *
737  * cell->get_dof_indices(local_dof_indices);
738  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
739  * {
740  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
741  * system_matrix.add(local_dof_indices[i],
742  * local_dof_indices[j],
743  * cell_matrix(i, j));
744  *
745  * system_rhs(local_dof_indices[i]) += cell_rhs(i);
746  * }
747  * }
748  *
749  * @endcode
750  *
751  * Finally, we remove hanging nodes from the system and apply zero
752  * boundary values to the linear system that defines the Newton updates
753  * @f$\delta u^n@f$:
754  *
755  * @code
756  * hanging_node_constraints.condense(system_matrix);
757  * hanging_node_constraints.condense(system_rhs);
758  *
759  * std::map<types::global_dof_index, double> boundary_values;
760  * VectorTools::interpolate_boundary_values(dof_handler,
761  * 0,
762  * Functions::ZeroFunction<dim>(),
763  * boundary_values);
764  * MatrixTools::apply_boundary_values(boundary_values,
765  * system_matrix,
766  * newton_update,
767  * system_rhs);
768  * }
769  *
770  *
771  *
772  * @endcode
773  *
774  *
775  * <a name="MinimalSurfaceProblemsolve"></a>
776  * <h4>MinimalSurfaceProblem::solve</h4>
777  *
778 
779  *
780  * The solve function is the same as always. At the end of the solution
781  * process we update the current solution by setting
782  * @f$u^{n+1}=u^n+\alpha^n\;\delta u^n@f$.
783  *
784  * @code
785  * template <int dim>
786  * void MinimalSurfaceProblem<dim>::solve()
787  * {
788  * SolverControl solver_control(system_rhs.size(),
789  * system_rhs.l2_norm() * 1e-6);
790  * SolverCG<Vector<double>> solver(solver_control);
791  *
792  * PreconditionSSOR<SparseMatrix<double>> preconditioner;
793  * preconditioner.initialize(system_matrix, 1.2);
794  *
795  * solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
796  *
797  * hanging_node_constraints.distribute(newton_update);
798  *
799  * const double alpha = determine_step_length();
800  * present_solution.add(alpha, newton_update);
801  * }
802  *
803  *
804  * @endcode
805  *
806  *
807  * <a name="MinimalSurfaceProblemrefine_mesh"></a>
808  * <h4>MinimalSurfaceProblem::refine_mesh</h4>
809  *
810 
811  *
812  * The first part of this function is the same as in @ref step_6 "step-6"... However,
813  * after refining the mesh we have to transfer the old solution to the new
814  * one which we do with the help of the SolutionTransfer class. The process
815  * is slightly convoluted, so let us describe it in detail:
816  *
817  * @code
818  * template <int dim>
819  * void MinimalSurfaceProblem<dim>::refine_mesh()
820  * {
821  * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
822  *
823  * KellyErrorEstimator<dim>::estimate(
824  * dof_handler,
825  * QGauss<dim - 1>(fe.degree + 1),
826  * std::map<types::boundary_id, const Function<dim> *>(),
827  * present_solution,
828  * estimated_error_per_cell);
829  *
830  * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
831  * estimated_error_per_cell,
832  * 0.3,
833  * 0.03);
834  *
835  * @endcode
836  *
837  * Then we need an additional step: if, for example, you flag a cell that
838  * is once more refined than its neighbor, and that neighbor is not
839  * flagged for refinement, we would end up with a jump of two refinement
840  * levels across a cell interface. To avoid these situations, the library
841  * will silently also have to refine the neighbor cell once. It does so by
842  * calling the Triangulation::prepare_coarsening_and_refinement function
843  * before actually doing the refinement and coarsening. This function
844  * flags a set of additional cells for refinement or coarsening, to
845  * enforce rules like the one-hanging-node rule. The cells that are
846  * flagged for refinement and coarsening after calling this function are
847  * exactly the ones that will actually be refined or coarsened. Usually,
848  * you don't have to do this by hand
850  * you). However, we need to initialize the SolutionTransfer class and it
851  * needs to know the final set of cells that will be coarsened or refined
852  * in order to store the data from the old mesh and transfer to the new
853  * one. Thus, we call the function by hand:
854  *
855  * @code
856  * triangulation.prepare_coarsening_and_refinement();
857  *
858  * @endcode
859  *
860  * With this out of the way, we initialize a SolutionTransfer object with
861  * the present DoFHandler and attach the solution vector to it, followed
862  * by doing the actual refinement and distribution of degrees of freedom
863  * on the new mesh
864  *
865  * @code
866  * SolutionTransfer<dim> solution_transfer(dof_handler);
867  * solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
868  *
869  * triangulation.execute_coarsening_and_refinement();
870  *
871  * dof_handler.distribute_dofs(fe);
872  *
873  * @endcode
874  *
875  * Finally, we retrieve the old solution interpolated to the new
876  * mesh. Since the SolutionTransfer function does not actually store the
877  * values of the old solution, but rather indices, we need to preserve the
878  * old solution vector until we have gotten the new interpolated
879  * values. Thus, we have the new values written into a temporary vector,
880  * and only afterwards write them into the solution vector object. Once we
881  * have this solution we have to make sure that the @f$u^n@f$ we now have
882  * actually has the correct boundary values. As explained at the end of
883  * the introduction, this is not automatically the case even if the
884  * solution before refinement had the correct boundary values, and so we
885  * have to explicitly make sure that it now has:
886  *
887  * @code
888  * Vector<double> tmp(dof_handler.n_dofs());
889  * solution_transfer.interpolate(present_solution, tmp);
890  * present_solution = tmp;
891  *
892  * set_boundary_values();
893  *
894  * @endcode
895  *
896  * On the new mesh, there are different hanging nodes, which we have to
897  * compute again. To ensure there are no hanging nodes of the old mesh in
898  * the object, it's first cleared. To be on the safe side, we then also
899  * make sure that the current solution's vector entries satisfy the
900  * hanging node constraints (see the discussion in the documentation of
901  * the SolutionTransfer class for why this is necessary):
902  *
903  * @code
904  * hanging_node_constraints.clear();
905  *
907  * hanging_node_constraints);
908  * hanging_node_constraints.close();
909  *
910  * hanging_node_constraints.distribute(present_solution);
911  *
912  * @endcode
913  *
914  * We end the function by updating all the remaining data structures,
915  * indicating to <code>setup_dofs()</code> that this is not the first
916  * go-around and that it needs to preserve the content of the solution
917  * vector:
918  *
919  * @code
920  * setup_system(false);
921  * }
922  *
923  *
924  *
925  * @endcode
926  *
927  *
928  * <a name="MinimalSurfaceProblemset_boundary_values"></a>
929  * <h4>MinimalSurfaceProblem::set_boundary_values</h4>
930  *
931 
932  *
933  * The next function ensures that the solution vector's entries respect the
934  * boundary values for our problem. Having refined the mesh (or just
935  * started computations), there might be new nodal points on the
936  * boundary. These have values that are simply interpolated from the
937  * previous mesh (or are just zero), instead of the correct boundary
938  * values. This is fixed up by setting all boundary nodes explicit to the
939  * right value:
940  *
941  * @code
942  * template <int dim>
943  * void MinimalSurfaceProblem<dim>::set_boundary_values()
944  * {
945  * std::map<types::global_dof_index, double> boundary_values;
946  * VectorTools::interpolate_boundary_values(dof_handler,
947  * 0,
948  * BoundaryValues<dim>(),
949  * boundary_values);
950  * for (auto &boundary_value : boundary_values)
951  * present_solution(boundary_value.first) = boundary_value.second;
952  * }
953  *
954  *
955  * @endcode
956  *
957  *
958  * <a name="MinimalSurfaceProblemcompute_residual"></a>
959  * <h4>MinimalSurfaceProblem::compute_residual</h4>
960  *
961 
962  *
963  * In order to monitor convergence, we need a way to compute the norm of the
964  * (discrete) residual, i.e., the norm of the vector
965  * @f$\left<F(u^n),\varphi_i\right>@f$ with @f$F(u)=-\nabla \cdot \left(
966  * \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right)@f$ as discussed in the
967  * introduction. It turns out that (although we don't use this feature in
968  * the current version of the program) one needs to compute the residual
969  * @f$\left<F(u^n+\alpha^n\;\delta u^n),\varphi_i\right>@f$ when determining
970  * optimal step lengths, and so this is what we implement here: the function
971  * takes the step length @f$\alpha^n@f$ as an argument. The original
972  * functionality is of course obtained by passing a zero as argument.
973  *
974 
975  *
976  * In the function below, we first set up a vector for the residual, and
977  * then a vector for the evaluation point @f$u^n+\alpha^n\;\delta u^n@f$. This
978  * is followed by the same boilerplate code we use for all integration
979  * operations:
980  *
981  * @code
982  * template <int dim>
983  * double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
984  * {
985  * Vector<double> residual(dof_handler.n_dofs());
986  *
987  * Vector<double> evaluation_point(dof_handler.n_dofs());
988  * evaluation_point = present_solution;
989  * evaluation_point.add(alpha, newton_update);
990  *
991  * const QGauss<dim> quadrature_formula(fe.degree + 1);
992  * FEValues<dim> fe_values(fe,
993  * quadrature_formula,
996  *
997  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
998  * const unsigned int n_q_points = quadrature_formula.size();
999  *
1000  * Vector<double> cell_residual(dofs_per_cell);
1001  * std::vector<Tensor<1, dim>> gradients(n_q_points);
1002  *
1003  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1004  *
1005  * for (const auto &cell : dof_handler.active_cell_iterators())
1006  * {
1007  * cell_residual = 0;
1008  * fe_values.reinit(cell);
1009  *
1010  * @endcode
1011  *
1012  * The actual computation is much as in
1013  * <code>assemble_system()</code>. We first evaluate the gradients of
1014  * @f$u^n+\alpha^n\,\delta u^n@f$ at the quadrature points, then compute
1015  * the coefficient @f$a_n@f$, and then plug it all into the formula for
1016  * the residual:
1017  *
1018  * @code
1019  * fe_values.get_function_gradients(evaluation_point, gradients);
1020  *
1021  *
1022  * for (unsigned int q = 0; q < n_q_points; ++q)
1023  * {
1024  * const double coeff = 1 / std::sqrt(1 + gradients[q] * gradients[q]);
1025  *
1026  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1027  * cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
1028  * * coeff // * a_n
1029  * * gradients[q] // * u_n
1030  * * fe_values.JxW(q)); // * dx
1031  * }
1032  *
1033  * cell->get_dof_indices(local_dof_indices);
1034  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1035  * residual(local_dof_indices[i]) += cell_residual(i);
1036  * }
1037  *
1038  * @endcode
1039  *
1040  * At the end of this function we also have to deal with the hanging node
1041  * constraints and with the issue of boundary values. With regard to the
1042  * latter, we have to set to zero the elements of the residual vector for
1043  * all entries that correspond to degrees of freedom that sit at the
1044  * boundary. The reason is that because the value of the solution there is
1045  * fixed, they are of course no "real" degrees of freedom and so, strictly
1046  * speaking, we shouldn't have assembled entries in the residual vector
1047  * for them. However, as we always do, we want to do exactly the same
1048  * thing on every cell and so we didn't not want to deal with the question
1049  * of whether a particular degree of freedom sits at the boundary in the
1050  * integration above. Rather, we will simply set to zero these entries
1051  * after the fact. To this end, we first need to determine which degrees
1052  * of freedom do in fact belong to the boundary and then loop over all of
1053  * those and set the residual entry to zero. This happens in the following
1054  * lines which we have already seen used in @ref step_11 "step-11":
1055  *
1056  * @code
1057  * hanging_node_constraints.condense(residual);
1058  *
1059  * std::vector<bool> boundary_dofs(dof_handler.n_dofs());
1060  * DoFTools::extract_boundary_dofs(dof_handler,
1061  * ComponentMask(),
1062  * boundary_dofs);
1063  * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1064  * if (boundary_dofs[i] == true)
1065  * residual(i) = 0;
1066  *
1067  * @endcode
1068  *
1069  * At the end of the function, we return the norm of the residual:
1070  *
1071  * @code
1072  * return residual.l2_norm();
1073  * }
1074  *
1075  *
1076  *
1077  * @endcode
1078  *
1079  *
1080  * <a name="MinimalSurfaceProblemdetermine_step_length"></a>
1081  * <h4>MinimalSurfaceProblem::determine_step_length</h4>
1082  *
1083 
1084  *
1085  * As discussed in the introduction, Newton's method frequently does not
1086  * converge if we always take full steps, i.e., compute @f$u^{n+1}=u^n+\delta
1087  * u^n@f$. Rather, one needs a damping parameter (step length) @f$\alpha^n@f$ and
1088  * set @f$u^{n+1}=u^n+\alpha^n\delta u^n@f$. This function is the one called
1089  * to compute @f$\alpha^n@f$.
1090  *
1091 
1092  *
1093  * Here, we simply always return 0.1. This is of course a sub-optimal
1094  * choice: ideally, what one wants is that the step size goes to one as we
1095  * get closer to the solution, so that we get to enjoy the rapid quadratic
1096  * convergence of Newton's method. We will discuss better strategies below
1097  * in the results section.
1098  *
1099  * @code
1100  * template <int dim>
1101  * double MinimalSurfaceProblem<dim>::determine_step_length() const
1102  * {
1103  * return 0.1;
1104  * }
1105  *
1106  *
1107  *
1108  * @endcode
1109  *
1110  *
1111  * <a name="MinimalSurfaceProblemrun"></a>
1112  * <h4>MinimalSurfaceProblem::run</h4>
1113  *
1114 
1115  *
1116  * In the run function, we build the first grid and then have the top-level
1117  * logic for the Newton iteration. The function has two variables, one that
1118  * indicates whether this is the first time we solve for a Newton update and
1119  * one that indicates the refinement level of the mesh:
1120  *
1121  * @code
1122  * template <int dim>
1124  * {
1125  * unsigned int refinement = 0;
1126  * bool first_step = true;
1127  *
1128  * @endcode
1129  *
1130  * As described in the introduction, the domain is the unit disk around
1131  * the origin, created in the same way as shown in @ref step_6 "step-6". The mesh is
1132  * globally refined twice followed later on by several adaptive cycles:
1133  *
1134  * @code
1136  * triangulation.refine_global(2);
1137  *
1138  * @endcode
1139  *
1140  * The Newton iteration starts next. During the first step we do not have
1141  * information about the residual prior to this step and so we continue
1142  * the Newton iteration until we have reached at least one iteration and
1143  * until residual is less than @f$10^{-3}@f$.
1144  *
1145 
1146  *
1147  * At the beginning of the loop, we do a bit of setup work. In the first
1148  * go around, we compute the solution on the twice globally refined mesh
1149  * after setting up the basic data structures and ensuring that the first
1150  * Newton iterate already has the correct boundary values. In all
1151  * following mesh refinement loops, the mesh will be refined adaptively.
1152  *
1153  * @code
1154  * double previous_res = 0;
1155  * while (first_step || (previous_res > 1e-3))
1156  * {
1157  * if (first_step == true)
1158  * {
1159  * std::cout << "******** Initial mesh "
1160  * << " ********" << std::endl;
1161  *
1162  * setup_system(true);
1163  * set_boundary_values();
1164  *
1165  * first_step = false;
1166  * }
1167  * else
1168  * {
1169  * ++refinement;
1170  * std::cout << "******** Refined mesh " << refinement << " ********"
1171  * << std::endl;
1172  *
1173  * refine_mesh();
1174  * }
1175  *
1176  * @endcode
1177  *
1178  * On every mesh we do exactly five Newton steps. We print the initial
1179  * residual here and then start the iterations on this mesh.
1180  *
1181 
1182  *
1183  * In every Newton step the system matrix and the right hand side have
1184  * to be computed first, after which we store the norm of the right
1185  * hand side as the residual to check against when deciding whether to
1186  * stop the iterations. We then solve the linear system (the function
1187  * also updates @f$u^{n+1}=u^n+\alpha^n\;\delta u^n@f$) and output the
1188  * residual at the end of this Newton step:
1189  *
1190  * @code
1191  * std::cout << " Initial residual: " << compute_residual(0) << std::endl;
1192  *
1193  * for (unsigned int inner_iteration = 0; inner_iteration < 5;
1194  * ++inner_iteration)
1195  * {
1196  * assemble_system();
1197  * previous_res = system_rhs.l2_norm();
1198  *
1199  * solve();
1200  *
1201  * std::cout << " Residual: " << compute_residual(0) << std::endl;
1202  * }
1203  *
1204  * @endcode
1205  *
1206  * Every fifth iteration, i.e., just before we refine the mesh again,
1207  * we output the solution as well as the Newton update. This happens
1208  * as in all programs before:
1209  *
1210  * @code
1211  * DataOut<dim> data_out;
1212  *
1213  * data_out.attach_dof_handler(dof_handler);
1214  * data_out.add_data_vector(present_solution, "solution");
1215  * data_out.add_data_vector(newton_update, "update");
1216  * data_out.build_patches();
1217  * const std::string filename =
1218  * "solution-" + Utilities::int_to_string(refinement, 2) + ".vtk";
1219  * std::ofstream output(filename);
1220  * DataOutBase::VtkFlags vtk_flags;
1221  * vtk_flags.compression_level =
1222  * DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
1223  * data_out.set_flags(vtk_flags);
1224  * data_out.write_vtu(output);
1225  * }
1226  * }
1227  * } // namespace Step15
1228  *
1229  * @endcode
1230  *
1231  *
1232  * <a name="Themainfunction"></a>
1233  * <h4>The main function</h4>
1234  *
1235 
1236  *
1237  * Finally the main function. This follows the scheme of all other main
1238  * functions:
1239  *
1240  * @code
1241  * int main()
1242  * {
1243  * try
1244  * {
1245  * using namespace Step15;
1246  *
1247  * MinimalSurfaceProblem<2> laplace_problem_2d;
1248  * laplace_problem_2d.run();
1249  * }
1250  * catch (std::exception &exc)
1251  * {
1252  * std::cerr << std::endl
1253  * << std::endl
1254  * << "----------------------------------------------------"
1255  * << std::endl;
1256  * std::cerr << "Exception on processing: " << std::endl
1257  * << exc.what() << std::endl
1258  * << "Aborting!" << std::endl
1259  * << "----------------------------------------------------"
1260  * << std::endl;
1261  *
1262  * return 1;
1263  * }
1264  * catch (...)
1265  * {
1266  * std::cerr << std::endl
1267  * << std::endl
1268  * << "----------------------------------------------------"
1269  * << std::endl;
1270  * std::cerr << "Unknown exception!" << std::endl
1271  * << "Aborting!" << std::endl
1272  * << "----------------------------------------------------"
1273  * << std::endl;
1274  * return 1;
1275  * }
1276  * return 0;
1277  * }
1278  * @endcode
1279 <a name="Results"></a><h1>Results</h1>
1280 
1281 
1282 
1283 The output of the program looks as follows:
1284 @code
1285 * ******** Initial mesh ********
1286  Initial residual: 1.53143
1287  Residual: 1.08746
1288  Residual: 0.966748
1289  Residual: 0.859602
1290  Residual: 0.766462
1291  Residual: 0.685475
1292 * ******** Refined mesh 1 ********
1293  Initial residual: 0.865774
1294  Residual: 0.759295
1295  Residual: 0.675281
1296  Residual: 0.603523
1297  Residual: 0.540744
1298  Residual: 0.485238
1299 * ******** Refined mesh 2 ********
1300  Initial residual: 0.425581
1301  Residual: 0.382042
1302  Residual: 0.343307
1303  Residual: 0.308718
1304 ....
1305 @endcode
1306 
1307 Obviously, the scheme converges, if not very fast. We will come back to
1308 strategies for accelerating the method below.
1309 
1310 One can visualize the solution after each set of five Newton
1311 iterations, i.e., on each of the meshes on which we approximate the
1312 solution. This yields the following set of images:
1313 
1314 <div class="twocolumn" style="width: 80%">
1315  <div>
1316  <img src="https://www.dealii.org/images/steps/developer/step_15_solution_1.png"
1317  alt="Solution after zero cycles with contour lines." width="230" height="273">
1318  </div>
1319  <div>
1320  <img src="https://www.dealii.org/images/steps/developer/step_15_solution_2.png"
1321  alt="Solution after one cycle with contour lines." width="230" height="273">
1322  </div>
1323  <div>
1324  <img src="https://www.dealii.org/images/steps/developer/step_15_solution_3.png"
1325  alt="Solution after two cycles with contour lines." width="230" height="273">
1326  </div>
1327  <div>
1328  <img src="https://www.dealii.org/images/steps/developer/step_15_solution_4.png"
1329  alt="Solution after three cycles with contour lines." width="230" height="273">
1330  </div>
1331 </div>
1332 
1333 It is clearly visible, that the solution minimizes the surface
1334 after each refinement. The solution converges to a picture one
1335 would imagine a soap bubble to be that is located inside a wire loop
1336 that is bent like
1337 the boundary. Also it is visible, how the boundary
1338 is smoothed out after each refinement. On the coarse mesh,
1339 the boundary doesn't look like a sine, whereas it does the
1340 finer the mesh gets.
1341 
1342 The mesh is mostly refined near the boundary, where the solution
1343 increases or decreases strongly, whereas it is coarsened on
1344 the inside of the domain, where nothing interesting happens,
1345 because there isn't much change in the solution. The ninth
1346 solution and mesh are shown here:
1347 
1348 <div class="onecolumn" style="width: 60%">
1349  <div>
1350  <img src="https://www.dealii.org/images/steps/developer/step_15_solution_9.png"
1351  alt="Grid and solution of the ninth cycle with contour lines." width="507" height="507">
1352  </div>
1353 </div>
1354 
1355 
1356 
1357 <a name="extensions"></a>
1358 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1359 
1360 
1361 The program shows the basic structure of a solver for a nonlinear, stationary
1362 problem. However, it does not converge particularly fast, for good reasons:
1363 
1364 - The program always takes a step size of 0.1. This precludes the rapid,
1365  quadratic convergence for which Newton's method is typically chosen.
1366 - It does not connect the nonlinear iteration with the mesh refinement
1367  iteration.
1368 
1369 Obviously, a better program would have to address these two points.
1370 We will discuss them in the following.
1371 
1372 
1373 <a name="Steplengthcontrol"></a><h4> Step length control </h4>
1374 
1375 
1376 Newton's method has two well known properties:
1377 - It does not converge from arbitrarily chosen starting points. Rather, a
1378  starting point has to be close enough to the solution to guarantee
1379  convergence. However, we can enlarge the area from which Newton's method
1380  converges by damping the iteration using a <i>step length</i> 0<@f$\alpha^n\le
1381  1@f$.
1382 - It exhibits rapid convergence of quadratic order if (i) the step length is
1383  chosen as @f$\alpha^n=1@f$, and (ii) it does in fact converge with this choice
1384  of step length.
1385 
1386 A consequence of these two observations is that a successful strategy is to
1387 choose @f$\alpha^n<1@f$ for the initial iterations until the iterate has come
1388 close enough to allow for convergence with full step length, at which point we
1389 want to switch to @f$\alpha^n=1@f$. The question is how to choose @f$\alpha^n@f$ in an
1390 automatic fashion that satisfies these criteria.
1391 
1392 We do not want to review the literature on this topic here, but only briefly
1393 mention that there are two fundamental approaches to the problem: backtracking
1394 line search and trust region methods. The former is more widely used for
1395 partial differential equations and essentially does the following:
1396 - Compute a search direction
1397 - See if the resulting residual of @f$u^n + \alpha^n\;\delta u^n@f$ with
1398  @f$\alpha^n=1@f$ is "substantially smaller" than that of @f$u^n@f$ alone.
1399 - If so, then take @f$\alpha^n=1@f$.
1400 - If not, try whether the residual is "substantially smaller" with
1401  @f$\alpha^n=2/3@f$.
1402 - If so, then take @f$\alpha^n=2/3@f$.
1403 - If not, try whether the residual is "substantially smaller" with
1404  @f$\alpha^n=(2/3)^2@f$.
1405 - Etc.
1406 One can of course choose other factors @f$r, r^2, \ldots@f$ than the @f$2/3,
1407 (2/3)^2, \ldots@f$ chosen above, for @f$0<r<1@f$. It is obvious where the term
1408 "backtracking" comes from: we try a long step, but if that doesn't work we try
1409 a shorter step, and ever shorter step, etc. The function
1410 <code>determine_step_length()</code> is written the way it is to support
1411 exactly this kind of use case.
1412 
1413 Whether we accept a particular step length @f$\alpha^n@f$ depends on how we define
1414 "substantially smaller". There are a number of ways to do so, but without
1415 going into detail let us just mention that the most common ones are to use the
1416 Wolfe and Armijo-Goldstein conditions. For these, one can show the following:
1417 - There is always a step length @f$\alpha^n@f$ for which the conditions are
1418  satisfied, i.e., the iteration never gets stuck as long as the problem is
1419  convex.
1420 - If we are close enough to the solution, then the conditions allow for
1421  @f$\alpha^n=1@f$, thereby enabling quadratic convergence.
1422 
1423 We will not dwell on this here any further but leave the implementation of
1424 such algorithms as an exercise. We note, however, that when implemented
1425 correctly then it is a common observation that most reasonably nonlinear
1426 problems can be solved in anywhere between 5 and 15 Newton iterations to
1427 engineering accuracy &mdash; substantially fewer than we need with the current
1428 version of the program.
1429 
1430 More details on globalization methods including backtracking can be found,
1431 for example, in Griva, Nash, Sofer: Linear and nonlinear optimization (2009).
1432 
1433 
1434 <a name="Integratingmeshrefinementandnonlinearandlinearsolvers"></a><h4> Integrating mesh refinement and nonlinear and linear solvers </h4>
1435 
1436 
1437 We currently do exactly 5 iterations on each mesh. But is this optimal? One
1438 could ask the following questions:
1439 - Maybe it is worthwhile doing more iterations on the initial meshes since
1440  there, computations are cheap.
1441 - On the other hand, we do not want to do too many iterations on every mesh:
1442  yes, we could drive the residual to zero on every mesh, but that would only
1443  mean that the nonlinear iteration error is far smaller than the
1444  discretization error.
1445 - Should we use solve the linear systems in each Newton step with higher or
1446  lower accuracy?
1447 
1448 Ultimately, what this boils down to is that we somehow need to couple the
1449 discretization error on the current mesh with the nonlinear residual we want
1450 to achieve with the Newton iterations on a given mesh, and to the linear
1451 iteration we want to achieve with the CG method within each Newton
1452 iterations.
1453 
1454 How to do this is, again, not entirely trivial, and we again leave it as a
1455 future exercise.
1456  *
1457  *
1458 <a name="PlainProg"></a>
1459 <h1> The plain program</h1>
1460 @include "step-15.cc"
1461 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
AdaptationStrategies::Refinement::preserve
std::vector< value_type > preserve(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
ComponentMask
Definition: component_mask.h:83
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
DataOutBase::VtkFlags::compression_level
ZlibCompressionLevel compression_level
Definition: data_out_base.h:1159
DoFTools::always
@ always
Definition: dof_tools.h:236
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)
Definition: advection.h:80
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler
Definition: dof_handler.h:205
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
LocalIntegrators::Advection::cell_residual
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim >> &input, const ArrayView< const std::vector< double >> &velocity, double factor=1.)
Definition: advection.h:136
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
Triangulation::execute_coarsening_and_refinement
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:13385
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
SolutionTransfer
Definition: solution_transfer.h:340
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
VectorTools::mean
@ mean
Definition: vector_tools_common.h:79
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
DataOutBase::VtkFlags
Definition: data_out_base.h:1095
QGauss
Definition: quadrature_lib.h:40
Tensor::clear
constexpr void clear()
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
value
static const bool value
Definition: dof_tools_constraints.cc:433
GridGenerator::hyper_ball
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
std::sqrt
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
DoFTools::extract_boundary_dofs
void extract_boundary_dofs(const DoFHandlerType &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
Definition: dof_tools.cc:637
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
DerivativeApproximation::internal::approximate
void approximate(SynchronousIterators< std::tuple< TriaActiveIterator< ::DoFCellAccessor< DoFHandlerType< dim, spacedim >, false >>, Vector< float >::iterator >> const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Definition: derivative_approximation.cc:924
FullMatrix< double >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector< double >
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41
std::sin
inline ::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5297