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-14.h
Go to the documentation of this file.
1 ) const
752  * {
753  * std::ofstream out(output_name_base + "-" +
754  * std::to_string(this->refinement_cycle) + ".svg");
755  * GridOut().write_svg(dof_handler.get_triangulation(), out);
756  * }
757  * } // namespace Evaluation
758  *
759  *
760  * @endcode
761  *
762  *
763  * <a name="TheLaplacesolverclasses"></a>
764  * <h3>The Laplace solver classes</h3>
765  *
766 
767  *
768  * Next are the actual solver classes. Again, we discuss only the
769  * differences to the previous program.
770  *
771  * @code
772  * namespace LaplaceSolver
773  * {
774  * @endcode
775  *
776  *
777  * <a name="TheLaplacesolverbaseclass"></a>
778  * <h4>The Laplace solver base class</h4>
779  *
780 
781  *
782  * This class is almost unchanged, with the exception that it declares two
783  * more functions: <code>output_solution</code> will be used to generate
784  * output files from the actual solutions computed by derived classes, and
785  * the <code>set_refinement_cycle</code> function by which the testing
786  * framework sets the number of the refinement cycle to a local variable
787  * in this class; this number is later used to generate filenames for the
788  * solution output.
789  *
790  * @code
791  * template <int dim>
792  * class Base
793  * {
794  * public:
795  * Base(Triangulation<dim> &coarse_grid);
796  * virtual ~Base() = default;
797  *
798  * virtual void solve_problem() = 0;
799  * virtual void postprocess(
800  * const Evaluation::EvaluationBase<dim> &postprocessor) const = 0;
801  * virtual void refine_grid() = 0;
802  * virtual unsigned int n_dofs() const = 0;
803  *
804  * virtual void set_refinement_cycle(const unsigned int cycle);
805  *
806  * virtual void output_solution() const = 0;
807  *
808  * protected:
810  *
811  * unsigned int refinement_cycle;
812  * };
813  *
814  *
815  * template <int dim>
816  * Base<dim>::Base(Triangulation<dim> &coarse_grid)
817  * : triangulation(&coarse_grid)
818  * , refinement_cycle(numbers::invalid_unsigned_int)
819  * {}
820  *
821  *
822  *
823  * template <int dim>
824  * void Base<dim>::set_refinement_cycle(const unsigned int cycle)
825  * {
826  * refinement_cycle = cycle;
827  * }
828  *
829  *
830  * @endcode
831  *
832  *
833  * <a name="TheLaplaceSolverclass"></a>
834  * <h4>The Laplace Solver class</h4>
835  *
836 
837  *
838  * Likewise, the <code>Solver</code> class is entirely unchanged and will
839  * thus not be discussed.
840  *
841  * @code
842  * template <int dim>
843  * class Solver : public virtual Base<dim>
844  * {
845  * public:
847  * const FiniteElement<dim> & fe,
848  * const Quadrature<dim> & quadrature,
849  * const Quadrature<dim - 1> &face_quadrature,
850  * const Function<dim> & boundary_values);
851  * virtual ~Solver() override;
852  *
853  * virtual void solve_problem() override;
854  *
855  * virtual void postprocess(
856  * const Evaluation::EvaluationBase<dim> &postprocessor) const override;
857  *
858  * virtual unsigned int n_dofs() const override;
859  *
860  * protected:
862  * const SmartPointer<const Quadrature<dim>> quadrature;
863  * const SmartPointer<const Quadrature<dim - 1>> face_quadrature;
864  * DoFHandler<dim> dof_handler;
865  * Vector<double> solution;
866  * const SmartPointer<const Function<dim>> boundary_values;
867  *
868  * virtual void assemble_rhs(Vector<double> &rhs) const = 0;
869  *
870  * private:
871  * struct LinearSystem
872  * {
873  * LinearSystem(const DoFHandler<dim> &dof_handler);
874  *
875  * void solve(Vector<double> &solution) const;
876  *
877  * AffineConstraints<double> hanging_node_constraints;
878  * SparsityPattern sparsity_pattern;
880  * Vector<double> rhs;
881  * };
882  *
883  *
884  * @endcode
885  *
886  * The remainder of the class is essentially a copy of @ref step_13 "step-13"
887  * as well, including the data structures and functions
888  * necessary to compute the linear system in parallel using the
889  * WorkStream framework:
890  *
891  * @code
892  * struct AssemblyScratchData
893  * {
894  * AssemblyScratchData(const FiniteElement<dim> &fe,
895  * const Quadrature<dim> & quadrature);
896  * AssemblyScratchData(const AssemblyScratchData &scratch_data);
897  *
898  * FEValues<dim> fe_values;
899  * };
900  *
901  * struct AssemblyCopyData
902  * {
904  * std::vector<types::global_dof_index> local_dof_indices;
905  * };
906  *
907  *
908  * void assemble_linear_system(LinearSystem &linear_system);
909  *
910  * void local_assemble_matrix(
911  * const typename DoFHandler<dim>::active_cell_iterator &cell,
912  * AssemblyScratchData & scratch_data,
913  * AssemblyCopyData & copy_data) const;
914  *
915  *
916  * void copy_local_to_global(const AssemblyCopyData &copy_data,
917  * LinearSystem & linear_system) const;
918  * };
919  *
920  *
921  *
922  * template <int dim>
924  * const FiniteElement<dim> & fe,
925  * const Quadrature<dim> & quadrature,
926  * const Quadrature<dim - 1> &face_quadrature,
927  * const Function<dim> & boundary_values)
928  * : Base<dim>(triangulation)
929  * , fe(&fe)
930  * , quadrature(&quadrature)
931  * , face_quadrature(&face_quadrature)
932  * , dof_handler(triangulation)
933  * , boundary_values(&boundary_values)
934  * {}
935  *
936  *
937  * template <int dim>
939  * {
940  * dof_handler.clear();
941  * }
942  *
943  *
944  * template <int dim>
946  * {
947  * dof_handler.distribute_dofs(*fe);
948  * solution.reinit(dof_handler.n_dofs());
949  *
950  * LinearSystem linear_system(dof_handler);
951  * assemble_linear_system(linear_system);
952  * linear_system.solve(solution);
953  * }
954  *
955  *
956  * template <int dim>
958  * const Evaluation::EvaluationBase<dim> &postprocessor) const
959  * {
960  * postprocessor(dof_handler, solution);
961  * }
962  *
963  *
964  * template <int dim>
965  * unsigned int Solver<dim>::n_dofs() const
966  * {
967  * return dof_handler.n_dofs();
968  * }
969  *
970  *
971  * @endcode
972  *
973  * The following few functions and constructors are verbatim
974  * copies taken from @ref step_13 "step-13":
975  *
976  * @code
977  * template <int dim>
978  * void Solver<dim>::assemble_linear_system(LinearSystem &linear_system)
979  * {
980  * Threads::Task<void> rhs_task =
981  * Threads::new_task(&Solver<dim>::assemble_rhs, *this, linear_system.rhs);
982  *
983  * auto worker =
984  * [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
985  * AssemblyScratchData &scratch_data,
986  * AssemblyCopyData & copy_data) {
987  * this->local_assemble_matrix(cell, scratch_data, copy_data);
988  * };
989  *
990  * auto copier = [this, &linear_system](const AssemblyCopyData &copy_data) {
991  * this->copy_local_to_global(copy_data, linear_system);
992  * };
993  *
994  * WorkStream::run(dof_handler.begin_active(),
995  * dof_handler.end(),
996  * worker,
997  * copier,
998  * AssemblyScratchData(*fe, *quadrature),
999  * AssemblyCopyData());
1000  * linear_system.hanging_node_constraints.condense(linear_system.matrix);
1001  *
1002  * std::map<types::global_dof_index, double> boundary_value_map;
1004  * 0,
1005  * *boundary_values,
1006  * boundary_value_map);
1007  *
1008  * rhs_task.join();
1009  * linear_system.hanging_node_constraints.condense(linear_system.rhs);
1010  *
1011  * MatrixTools::apply_boundary_values(boundary_value_map,
1012  * linear_system.matrix,
1013  * solution,
1014  * linear_system.rhs);
1015  * }
1016  *
1017  *
1018  * template <int dim>
1020  * const FiniteElement<dim> &fe,
1021  * const Quadrature<dim> & quadrature)
1022  * : fe_values(fe, quadrature, update_gradients | update_JxW_values)
1023  * {}
1024  *
1025  *
1026  * template <int dim>
1028  * const AssemblyScratchData &scratch_data)
1029  * : fe_values(scratch_data.fe_values.get_fe(),
1030  * scratch_data.fe_values.get_quadrature(),
1032  * {}
1033  *
1034  *
1035  * template <int dim>
1037  * const typename DoFHandler<dim>::active_cell_iterator &cell,
1038  * AssemblyScratchData & scratch_data,
1039  * AssemblyCopyData & copy_data) const
1040  * {
1041  * const unsigned int dofs_per_cell = fe->dofs_per_cell;
1042  * const unsigned int n_q_points = quadrature->size();
1043  *
1044  * copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1045  *
1046  * copy_data.local_dof_indices.resize(dofs_per_cell);
1047  *
1048  * scratch_data.fe_values.reinit(cell);
1049  *
1050  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1051  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1052  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1053  * copy_data.cell_matrix(i, j) +=
1054  * (scratch_data.fe_values.shape_grad(i, q_point) *
1055  * scratch_data.fe_values.shape_grad(j, q_point) *
1056  * scratch_data.fe_values.JxW(q_point));
1057  *
1058  * cell->get_dof_indices(copy_data.local_dof_indices);
1059  * }
1060  *
1061  *
1062  *
1063  * template <int dim>
1064  * void Solver<dim>::copy_local_to_global(const AssemblyCopyData &copy_data,
1065  * LinearSystem &linear_system) const
1066  * {
1067  * for (unsigned int i = 0; i < copy_data.local_dof_indices.size(); ++i)
1068  * for (unsigned int j = 0; j < copy_data.local_dof_indices.size(); ++j)
1069  * linear_system.matrix.add(copy_data.local_dof_indices[i],
1070  * copy_data.local_dof_indices[j],
1071  * copy_data.cell_matrix(i, j));
1072  * }
1073  *
1074  *
1075  * @endcode
1076  *
1077  * Now for the functions that implement actions in the linear
1078  * system class. First, the constructor initializes all data
1079  * elements to their correct sizes, and sets up a number of
1080  * additional data structures, such as constraints due to hanging
1081  * nodes. Since setting up the hanging nodes and finding out about
1082  * the nonzero elements of the matrix is independent, we do that
1083  * in parallel (if the library was configured to use concurrency,
1084  * at least; otherwise, the actions are performed
1085  * sequentially). Note that we start only one thread, and do the
1086  * second action in the main thread. Since only one thread is
1087  * generated, we don't use the <code>Threads::ThreadGroup</code>
1088  * class here, but rather use the one created thread object
1089  * directly to wait for this particular thread's exit. The
1090  * approach is generally the same as the one we have used in
1091  * <code>Solver::assemble_linear_system()</code> above.
1092  *
1093 
1094  *
1095  * Note that taking the address of the
1096  * <code>DoFTools::make_hanging_node_constraints</code> function
1097  * is a little tricky, since there are actually three functions of
1098  * this name, one for each supported space dimension. Taking
1099  * addresses of overloaded functions is somewhat complicated in
1100  * C++, since the address-of operator <code>&</code> in that case
1101  * returns a set of values (the addresses of all
1102  * functions with that name), and selecting the right one is then
1103  * the next step. If the context dictates which one to take (for
1104  * example by assigning to a function pointer of known type), then
1105  * the compiler can do that by itself, but if this set of pointers
1106  * shall be given as the argument to a function that takes a
1107  * template, the compiler could choose all without having a
1108  * preference for one. We therefore have to make it clear to the
1109  * compiler which one we would like to have; for this, we could
1110  * use a cast, but for more clarity, we assign it to a temporary
1111  * <code>mhnc_p</code> (short for <code>pointer to
1112  * make_hanging_node_constraints</code>) with the right type, and
1113  * using this pointer instead.
1114  *
1115  * @code
1116  * template <int dim>
1118  * {
1119  * hanging_node_constraints.clear();
1120  *
1121  * void (*mhnc_p)(const DoFHandler<dim> &, AffineConstraints<double> &) =
1123  *
1124  * @endcode
1125  *
1126  * Start a side task then continue on the main thread
1127  *
1128  * @code
1129  * Threads::Task<void> side_task =
1130  * Threads::new_task(mhnc_p, dof_handler, hanging_node_constraints);
1131  *
1132  * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1133  * DoFTools::make_sparsity_pattern(dof_handler, dsp);
1134  *
1135  *
1136  *
1137  * @endcode
1138  *
1139  * Wait for the side task to be done before going further
1140  *
1141  * @code
1142  * side_task.join();
1143  *
1144  * hanging_node_constraints.close();
1145  * hanging_node_constraints.condense(dsp);
1146  * sparsity_pattern.copy_from(dsp);
1147  *
1148  * matrix.reinit(sparsity_pattern);
1149  * rhs.reinit(dof_handler.n_dofs());
1150  * }
1151  *
1152  *
1153  *
1154  * template <int dim>
1155  * void Solver<dim>::LinearSystem::solve(Vector<double> &solution) const
1156  * {
1157  * SolverControl solver_control(5000, 1e-12);
1158  * SolverCG<Vector<double>> cg(solver_control);
1159  *
1160  * PreconditionSSOR<SparseMatrix<double>> preconditioner;
1161  * preconditioner.initialize(matrix, 1.2);
1162  *
1163  * cg.solve(matrix, solution, rhs, preconditioner);
1164  *
1165  * hanging_node_constraints.distribute(solution);
1166  * }
1167  *
1168  *
1169  *
1170  * @endcode
1171  *
1172  *
1173  * <a name="ThePrimalSolverclass"></a>
1174  * <h4>The PrimalSolver class</h4>
1175  *
1176 
1177  *
1178  * The <code>PrimalSolver</code> class is also mostly unchanged except for
1179  * implementing the <code>output_solution</code> function. We keep the
1180  * <code>GlobalRefinement</code> and <code>RefinementKelly</code> classes
1181  * in this program, and they can then rely on the default implementation
1182  * of this function which simply outputs the primal solution. The class
1183  * implementing dual weighted error estimators will overload this function
1184  * itself, to also output the dual solution.
1185  *
1186  * @code
1187  * template <int dim>
1188  * class PrimalSolver : public Solver<dim>
1189  * {
1190  * public:
1191  * PrimalSolver(Triangulation<dim> & triangulation,
1192  * const FiniteElement<dim> & fe,
1193  * const Quadrature<dim> & quadrature,
1194  * const Quadrature<dim - 1> &face_quadrature,
1195  * const Function<dim> & rhs_function,
1196  * const Function<dim> & boundary_values);
1197  *
1198  * virtual void output_solution() const override;
1199  *
1200  * protected:
1201  * const SmartPointer<const Function<dim>> rhs_function;
1202  * virtual void assemble_rhs(Vector<double> &rhs) const override;
1203  * };
1204  *
1205  *
1206  * template <int dim>
1207  * PrimalSolver<dim>::PrimalSolver(Triangulation<dim> & triangulation,
1208  * const FiniteElement<dim> & fe,
1209  * const Quadrature<dim> & quadrature,
1210  * const Quadrature<dim - 1> &face_quadrature,
1211  * const Function<dim> & rhs_function,
1212  * const Function<dim> & boundary_values)
1213  * : Base<dim>(triangulation)
1215  * fe,
1216  * quadrature,
1217  * face_quadrature,
1218  * boundary_values)
1219  * , rhs_function(&rhs_function)
1220  * {}
1221  *
1222  *
1223  *
1224  * template <int dim>
1225  * void PrimalSolver<dim>::output_solution() const
1226  * {
1227  * DataOut<dim> data_out;
1228  * data_out.attach_dof_handler(this->dof_handler);
1229  * data_out.add_data_vector(this->solution, "solution");
1230  * data_out.build_patches();
1231  *
1232  * std::ofstream out("solution-" + std::to_string(this->refinement_cycle) +
1233  * ".vtu");
1234  * data_out.write(out, DataOutBase::vtu);
1235  * }
1236  *
1237  *
1238  *
1239  * template <int dim>
1240  * void PrimalSolver<dim>::assemble_rhs(Vector<double> &rhs) const
1241  * {
1242  * FEValues<dim> fe_values(*this->fe,
1243  * *this->quadrature,
1245  * update_JxW_values);
1246  *
1247  * const unsigned int dofs_per_cell = this->fe->dofs_per_cell;
1248  * const unsigned int n_q_points = this->quadrature->size();
1249  *
1250  * Vector<double> cell_rhs(dofs_per_cell);
1251  * std::vector<double> rhs_values(n_q_points);
1252  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1253  *
1254  * for (const auto &cell : this->dof_handler.active_cell_iterators())
1255  * {
1256  * cell_rhs = 0;
1257  *
1258  * fe_values.reinit(cell);
1259  *
1260  * rhs_function->value_list(fe_values.get_quadrature_points(),
1261  * rhs_values);
1262  *
1263  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1264  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1265  * cell_rhs(i) += (fe_values.shape_value(i, q_point) * // phi_i(x_q)
1266  * rhs_values[q_point] * // f((x_q)
1267  * fe_values.JxW(q_point)); // dx
1268  *
1269  * cell->get_dof_indices(local_dof_indices);
1270  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1271  * rhs(local_dof_indices[i]) += cell_rhs(i);
1272  * }
1273  * }
1274  *
1275  *
1276  * @endcode
1277  *
1278  *
1279  * <a name="TheRefinementGlobalandRefinementKellyclasses"></a>
1280  * <h4>The RefinementGlobal and RefinementKelly classes</h4>
1281  *
1282 
1283  *
1284  * For the following two classes, the same applies as for most of the
1285  * above: the class is taken from the previous example as-is:
1286  *
1287  * @code
1288  * template <int dim>
1289  * class RefinementGlobal : public PrimalSolver<dim>
1290  * {
1291  * public:
1292  * RefinementGlobal(Triangulation<dim> & coarse_grid,
1293  * const FiniteElement<dim> & fe,
1294  * const Quadrature<dim> & quadrature,
1295  * const Quadrature<dim - 1> &face_quadrature,
1296  * const Function<dim> & rhs_function,
1297  * const Function<dim> & boundary_values);
1298  *
1299  * virtual void refine_grid() override;
1300  * };
1301  *
1302  *
1303  *
1304  * template <int dim>
1305  * RefinementGlobal<dim>::RefinementGlobal(
1306  * Triangulation<dim> & coarse_grid,
1307  * const FiniteElement<dim> & fe,
1308  * const Quadrature<dim> & quadrature,
1309  * const Quadrature<dim - 1> &face_quadrature,
1310  * const Function<dim> & rhs_function,
1311  * const Function<dim> & boundary_values)
1312  * : Base<dim>(coarse_grid)
1313  * , PrimalSolver<dim>(coarse_grid,
1314  * fe,
1315  * quadrature,
1316  * face_quadrature,
1317  * rhs_function,
1318  * boundary_values)
1319  * {}
1320  *
1321  *
1322  *
1323  * template <int dim>
1324  * void RefinementGlobal<dim>::refine_grid()
1325  * {
1326  * this->triangulation->refine_global(1);
1327  * }
1328  *
1329  *
1330  *
1331  * template <int dim>
1332  * class RefinementKelly : public PrimalSolver<dim>
1333  * {
1334  * public:
1335  * RefinementKelly(Triangulation<dim> & coarse_grid,
1336  * const FiniteElement<dim> & fe,
1337  * const Quadrature<dim> & quadrature,
1338  * const Quadrature<dim - 1> &face_quadrature,
1339  * const Function<dim> & rhs_function,
1340  * const Function<dim> & boundary_values);
1341  *
1342  * virtual void refine_grid() override;
1343  * };
1344  *
1345  *
1346  *
1347  * template <int dim>
1348  * RefinementKelly<dim>::RefinementKelly(
1349  * Triangulation<dim> & coarse_grid,
1350  * const FiniteElement<dim> & fe,
1351  * const Quadrature<dim> & quadrature,
1352  * const Quadrature<dim - 1> &face_quadrature,
1353  * const Function<dim> & rhs_function,
1354  * const Function<dim> & boundary_values)
1355  * : Base<dim>(coarse_grid)
1356  * , PrimalSolver<dim>(coarse_grid,
1357  * fe,
1358  * quadrature,
1359  * face_quadrature,
1360  * rhs_function,
1361  * boundary_values)
1362  * {}
1363  *
1364  *
1365  *
1366  * template <int dim>
1367  * void RefinementKelly<dim>::refine_grid()
1368  * {
1369  * Vector<float> estimated_error_per_cell(
1370  * this->triangulation->n_active_cells());
1372  * this->dof_handler,
1373  * QGauss<dim - 1>(this->fe->degree + 1),
1374  * std::map<types::boundary_id, const Function<dim> *>(),
1375  * this->solution,
1376  * estimated_error_per_cell);
1378  * estimated_error_per_cell,
1379  * 0.3,
1380  * 0.03);
1381  * this->triangulation->execute_coarsening_and_refinement();
1382  * }
1383  *
1384  *
1385  *
1386  * @endcode
1387  *
1388  *
1389  * <a name="TheRefinementWeightedKellyclass"></a>
1390  * <h4>The RefinementWeightedKelly class</h4>
1391  *
1392 
1393  *
1394  * This class is a variant of the previous one, in that it allows to
1395  * weight the refinement indicators we get from the library's Kelly
1396  * indicator by some function. We include this class since the goal of
1397  * this example program is to demonstrate automatic refinement criteria
1398  * even for complex output quantities such as point values or stresses. If
1399  * we did not solve a dual problem and compute the weights thereof, we
1400  * would probably be tempted to give a hand-crafted weighting to the
1401  * indicators to account for the fact that we are going to evaluate these
1402  * quantities. This class accepts such a weighting function as argument to
1403  * its constructor:
1404  *
1405  * @code
1406  * template <int dim>
1407  * class RefinementWeightedKelly : public PrimalSolver<dim>
1408  * {
1409  * public:
1410  * RefinementWeightedKelly(Triangulation<dim> & coarse_grid,
1411  * const FiniteElement<dim> & fe,
1412  * const Quadrature<dim> & quadrature,
1413  * const Quadrature<dim - 1> &face_quadrature,
1414  * const Function<dim> & rhs_function,
1415  * const Function<dim> & boundary_values,
1416  * const Function<dim> & weighting_function);
1417  *
1418  * virtual void refine_grid() override;
1419  *
1420  * private:
1421  * const SmartPointer<const Function<dim>> weighting_function;
1422  * };
1423  *
1424  *
1425  *
1426  * template <int dim>
1427  * RefinementWeightedKelly<dim>::RefinementWeightedKelly(
1428  * Triangulation<dim> & coarse_grid,
1429  * const FiniteElement<dim> & fe,
1430  * const Quadrature<dim> & quadrature,
1431  * const Quadrature<dim - 1> &face_quadrature,
1432  * const Function<dim> & rhs_function,
1433  * const Function<dim> & boundary_values,
1434  * const Function<dim> & weighting_function)
1435  * : Base<dim>(coarse_grid)
1436  * , PrimalSolver<dim>(coarse_grid,
1437  * fe,
1438  * quadrature,
1439  * face_quadrature,
1440  * rhs_function,
1441  * boundary_values)
1442  * , weighting_function(&weighting_function)
1443  * {}
1444  *
1445  *
1446  *
1447  * @endcode
1448  *
1449  * Now, here comes the main function, including the weighting:
1450  *
1451  * @code
1452  * template <int dim>
1453  * void RefinementWeightedKelly<dim>::refine_grid()
1454  * {
1455  * @endcode
1456  *
1457  * First compute some residual based error indicators for all cells by a
1458  * method already implemented in the library. What exactly we compute
1459  * here is described in more detail in the documentation of that class.
1460  *
1461  * @code
1462  * Vector<float> estimated_error_per_cell(
1463  * this->triangulation->n_active_cells());
1464  * std::map<types::boundary_id, const Function<dim> *> dummy_function_map;
1465  * KellyErrorEstimator<dim>::estimate(this->dof_handler,
1466  * *this->face_quadrature,
1467  * dummy_function_map,
1468  * this->solution,
1469  * estimated_error_per_cell);
1470  *
1471  * @endcode
1472  *
1473  * Next weigh each entry in the vector of indicators by the value of the
1474  * function given to the constructor, evaluated at the cell center. We
1475  * need to write the result into the vector entry that corresponds to the
1476  * current cell, which we can obtain by asking the cell what its index
1477  * among all active cells is using CellAccessor::active_cell_index(). (In
1478  * reality, this index is zero for the first cell we handle in the loop,
1479  * one for the second cell, etc., and we could as well just keep track of
1480  * this index using an integer counter; but using
1481  * CellAccessor::active_cell_index() makes this more explicit.)
1482  *
1483  * @code
1484  * for (const auto &cell : this->dof_handler.active_cell_iterators())
1485  * estimated_error_per_cell(cell->active_cell_index()) *=
1486  * weighting_function->value(cell->center());
1487  *
1488  * GridRefinement::refine_and_coarsen_fixed_number(*this->triangulation,
1489  * estimated_error_per_cell,
1490  * 0.3,
1491  * 0.03);
1492  * this->triangulation->execute_coarsening_and_refinement();
1493  * }
1494  *
1495  * } // namespace LaplaceSolver
1496  *
1497  *
1498  * @endcode
1499  *
1500  *
1501  * <a name="Equationdata"></a>
1502  * <h3>Equation data</h3>
1503  *
1504 
1505  *
1506  * In this example program, we work with the same data sets as in the
1507  * previous one, but as it may so happen that someone wants to run the
1508  * program with different boundary values and right hand side functions, or
1509  * on a different grid, we show a simple technique to do exactly that. For
1510  * more clarity, we furthermore pack everything that has to do with equation
1511  * data into a namespace of its own.
1512  *
1513 
1514  *
1515  * The underlying assumption is that this is a research program, and that
1516  * there we often have a number of test cases that consist of a domain, a
1517  * right hand side, boundary values, possibly a specified coefficient, and a
1518  * number of other parameters. They often vary all at the same time when
1519  * shifting from one example to another. To make handling such sets of
1520  * problem description parameters simple is the goal of the following.
1521  *
1522 
1523  *
1524  * Basically, the idea is this: let us have a structure for each set of
1525  * data, in which we pack everything that describes a test case: here, these
1526  * are two subclasses, one called <code>BoundaryValues</code> for the
1527  * boundary values of the exact solution, and one called
1528  * <code>RightHandSide</code>, and then a way to generate the coarse
1529  * grid. Since the solution of the previous example program looked like
1530  * curved ridges, we use this name here for the enclosing class. Note that
1531  * the names of the two inner classes have to be the same for all enclosing
1532  * test case classes, and also that we have attached the dimension template
1533  * argument to the enclosing class rather than to the inner ones, to make
1534  * further processing simpler. (From a language viewpoint, a namespace
1535  * would be better to encapsulate these inner classes, rather than a
1536  * structure. However, namespaces cannot be given as template arguments, so
1537  * we use a structure to allow a second object to select from within its
1538  * given argument. The enclosing structure, of course, has no member
1539  * variables apart from the classes it declares, and a static function to
1540  * generate the coarse mesh; it will in general never be instantiated.)
1541  *
1542 
1543  *
1544  * The idea is then the following (this is the right time to also take a
1545  * brief look at the code below): we can generate objects for boundary
1546  * values and right hand side by simply giving the name of the outer class
1547  * as a template argument to a class which we call here
1548  * <code>Data::SetUp</code>, and it then creates objects for the inner
1549  * classes. In this case, to get all that characterizes the curved ridge
1550  * solution, we would simply generate an instance of
1551  * <code>Data::SetUp@<Data::CurvedRidge@></code>, and everything we need to
1552  * know about the solution would be static member variables and functions of
1553  * that object.
1554  *
1555 
1556  *
1557  * This approach might seem like overkill in this case, but will become very
1558  * handy once a certain set up is not only characterized by Dirichlet
1559  * boundary values and a right hand side function, but in addition by
1560  * material properties, Neumann values, different boundary descriptors,
1561  * etc. In that case, the <code>SetUp</code> class might consist of a dozen
1562  * or more objects, and each descriptor class (like the
1563  * <code>CurvedRidges</code> class below) would have to provide them. Then,
1564  * you will be happy to be able to change from one set of data to another by
1565  * only changing the template argument to the <code>SetUp</code> class at
1566  * one place, rather than at many.
1567  *
1568 
1569  *
1570  * With this framework for different test cases, we are almost finished, but
1571  * one thing remains: by now we can select statically, by changing one
1572  * template argument, which data set to choose. In order to be able to do
1573  * that dynamically, i.e. at run time, we need a base class. This we provide
1574  * in the obvious way, see below, with virtual abstract functions. It forces
1575  * us to introduce a second template parameter <code>dim</code> which we
1576  * need for the base class (which could be avoided using some template
1577  * magic, but we omit that), but that's all.
1578  *
1579 
1580  *
1581  * Adding new testcases is now simple, you don't have to touch the framework
1582  * classes, only a structure like the <code>CurvedRidges</code> one is
1583  * needed.
1584  *
1585  * @code
1586  * namespace Data
1587  * {
1588  * @endcode
1589  *
1590  *
1591  * <a name="TheSetUpBaseandSetUpclasses"></a>
1592  * <h4>The SetUpBase and SetUp classes</h4>
1593  *
1594 
1595  *
1596  * Based on the above description, the <code>SetUpBase</code> class then
1597  * looks as follows. To allow using the <code>SmartPointer</code> class
1598  * with this class, we derived from the <code>Subscriptor</code> class.
1599  *
1600  * @code
1601  * template <int dim>
1602  * struct SetUpBase : public Subscriptor
1603  * {
1604  * virtual const Function<dim> &get_boundary_values() const = 0;
1605  *
1606  * virtual const Function<dim> &get_right_hand_side() const = 0;
1607  *
1608  * virtual void
1609  * create_coarse_grid(Triangulation<dim> &coarse_grid) const = 0;
1610  * };
1611  *
1612  *
1613  * @endcode
1614  *
1615  * And now for the derived class that takes the template argument as
1616  * explained above.
1617  *
1618 
1619  *
1620  * Here we pack the data elements into private variables, and allow access
1621  * to them through the methods of the base class.
1622  *
1623  * @code
1624  * template <class Traits, int dim>
1625  * struct SetUp : public SetUpBase<dim>
1626  * {
1627  * virtual const Function<dim> &get_boundary_values() const override;
1628  *
1629  * virtual const Function<dim> &get_right_hand_side() const override;
1630  *
1631  *
1632  * virtual void
1633  * create_coarse_grid(Triangulation<dim> &coarse_grid) const override;
1634  *
1635  * private:
1636  * static const typename Traits::BoundaryValues boundary_values;
1637  * static const typename Traits::RightHandSide right_hand_side;
1638  * };
1639  *
1640  * @endcode
1641  *
1642  * We have to provide definitions for the static member variables of the
1643  * above class:
1644  *
1645  * @code
1646  * template <class Traits, int dim>
1647  * const typename Traits::BoundaryValues SetUp<Traits, dim>::boundary_values;
1648  * template <class Traits, int dim>
1649  * const typename Traits::RightHandSide SetUp<Traits, dim>::right_hand_side;
1650  *
1651  * @endcode
1652  *
1653  * And definitions of the member functions:
1654  *
1655  * @code
1656  * template <class Traits, int dim>
1657  * const Function<dim> &SetUp<Traits, dim>::get_boundary_values() const
1658  * {
1659  * return boundary_values;
1660  * }
1661  *
1662  *
1663  * template <class Traits, int dim>
1664  * const Function<dim> &SetUp<Traits, dim>::get_right_hand_side() const
1665  * {
1666  * return right_hand_side;
1667  * }
1668  *
1669  *
1670  * template <class Traits, int dim>
1671  * void SetUp<Traits, dim>::create_coarse_grid(
1672  * Triangulation<dim> &coarse_grid) const
1673  * {
1674  * Traits::create_coarse_grid(coarse_grid);
1675  * }
1676  *
1677  *
1678  * @endcode
1679  *
1680  *
1681  * <a name="TheCurvedRidgesclass"></a>
1682  * <h4>The CurvedRidges class</h4>
1683  *
1684 
1685  *
1686  * The class that is used to describe the boundary values and right hand
1687  * side of the <code>curved ridge</code> problem already used in the
1688  * @ref step_13 "step-13" example program is then like so:
1689  *
1690  * @code
1691  * template <int dim>
1692  * struct CurvedRidges
1693  * {
1694  * class BoundaryValues : public Function<dim>
1695  * {
1696  * public:
1697  * virtual double value(const Point<dim> & p,
1698  * const unsigned int component) const;
1699  * };
1700  *
1701  *
1702  * class RightHandSide : public Function<dim>
1703  * {
1704  * public:
1705  * virtual double value(const Point<dim> & p,
1706  * const unsigned int component) const;
1707  * };
1708  *
1709  * static void create_coarse_grid(Triangulation<dim> &coarse_grid);
1710  * };
1711  *
1712  *
1713  * template <int dim>
1714  * double CurvedRidges<dim>::BoundaryValues::value(
1715  * const Point<dim> &p,
1716  * const unsigned int /*component*/) const
1717  * {
1718  * double q = p(0);
1719  * for (unsigned int i = 1; i < dim; ++i)
1720  * q += std::sin(10 * p(i) + 5 * p(0) * p(0));
1721  * const double exponential = std::exp(q);
1722  * return exponential;
1723  * }
1724  *
1725  *
1726  *
1727  * template <int dim>
1728  * double CurvedRidges<dim>::RightHandSide::value(
1729  * const Point<dim> &p,
1730  * const unsigned int /*component*/) const
1731  * {
1732  * double q = p(0);
1733  * for (unsigned int i = 1; i < dim; ++i)
1734  * q += std::sin(10 * p(i) + 5 * p(0) * p(0));
1735  * const double u = std::exp(q);
1736  * double t1 = 1, t2 = 0, t3 = 0;
1737  * for (unsigned int i = 1; i < dim; ++i)
1738  * {
1739  * t1 += std::cos(10 * p(i) + 5 * p(0) * p(0)) * 10 * p(0);
1740  * t2 += 10 * std::cos(10 * p(i) + 5 * p(0) * p(0)) -
1741  * 100 * std::sin(10 * p(i) + 5 * p(0) * p(0)) * p(0) * p(0);
1742  * t3 += 100 * std::cos(10 * p(i) + 5 * p(0) * p(0)) *
1743  * std::cos(10 * p(i) + 5 * p(0) * p(0)) -
1744  * 100 * std::sin(10 * p(i) + 5 * p(0) * p(0));
1745  * }
1746  * t1 = t1 * t1;
1747  *
1748  * return -u * (t1 + t2 + t3);
1749  * }
1750  *
1751  *
1752  * template <int dim>
1753  * void CurvedRidges<dim>::create_coarse_grid(Triangulation<dim> &coarse_grid)
1754  * {
1755  * GridGenerator::hyper_cube(coarse_grid, -1, 1);
1756  * coarse_grid.refine_global(2);
1757  * }
1758  *
1759  *
1760  * @endcode
1761  *
1762  *
1763  * <a name="TheExercise_2_3class"></a>
1764  * <h4>The Exercise_2_3 class</h4>
1765  *
1766 
1767  *
1768  * This example program was written while giving practical courses for a
1769  * lecture on adaptive finite element methods and duality based error
1770  * estimates. For these courses, we had one exercise, which required to
1771  * solve the Laplace equation with constant right hand side on a square
1772  * domain with a square hole in the center, and zero boundary
1773  * values. Since the implementation of the properties of this problem is
1774  * so particularly simple here, lets do it. As the number of the exercise
1775  * was 2.3, we take the liberty to retain this name for the class as well.
1776  *
1777  * @code
1778  * template <int dim>
1779  * struct Exercise_2_3
1780  * {
1781  * @endcode
1782  *
1783  * We need a class to denote the boundary values of the problem. In this
1784  * case, this is simple: it's the zero function, so don't even declare a
1785  * class, just an alias:
1786  *
1787  * @code
1788  * using BoundaryValues = Functions::ZeroFunction<dim>;
1789  *
1790  * @endcode
1791  *
1792  * Second, a class that denotes the right hand side. Since they are
1793  * constant, just subclass the corresponding class of the library and be
1794  * done:
1795  *
1796  * @code
1797  * class RightHandSide : public Functions::ConstantFunction<dim>
1798  * {
1799  * public:
1800  * RightHandSide()
1801  * : Functions::ConstantFunction<dim>(1.)
1802  * {}
1803  * };
1804  *
1805  * @endcode
1806  *
1807  * Finally a function to generate the coarse grid. This is somewhat more
1808  * complicated here, see immediately below.
1809  *
1810  * @code
1811  * static void create_coarse_grid(Triangulation<dim> &coarse_grid);
1812  * };
1813  *
1814  *
1815  * @endcode
1816  *
1817  * As stated above, the grid for this example is the square [-1,1]^2 with
1818  * the square [-1/2,1/2]^2 as hole in it. We create the coarse grid as 4
1819  * times 4 cells with the middle four ones missing. To understand how
1820  * exactly the mesh is going to look, it may be simplest to just look
1821  * at the "Results" section of this tutorial program first. In general,
1822  * if you'd like to understand more about creating meshes either from
1823  * scratch by hand, as we do here, or using other techniques, you
1824  * should take a look at @ref step_49 "step-49".
1825  *
1826 
1827  *
1828  * Of course, the example has an extension to 3d, but since this function
1829  * cannot be written in a dimension independent way we choose not to
1830  * implement this here, but rather only specialize the template for
1831  * dim=2. If you compile the program for 3d, you'll get a message from the
1832  * linker that this function is not implemented for 3d, and needs to be
1833  * provided.
1834  *
1835 
1836  *
1837  * For the creation of this geometry, the library has no predefined
1838  * method. In this case, the geometry is still simple enough to do the
1839  * creation by hand, rather than using a mesh generator.
1840  *
1841  * @code
1842  * template <>
1843  * void Exercise_2_3<2>::create_coarse_grid(Triangulation<2> &coarse_grid)
1844  * {
1845  * @endcode
1846  *
1847  * We first define the space dimension, to allow those parts of the
1848  * function that are actually dimension independent to use this
1849  * variable. That makes it simpler if you later take this as a starting
1850  * point to implement a 3d version of this mesh. The next step is then
1851  * to have a list of vertices. Here, they are 24 (5 times 5, with the
1852  * middle one omitted). It is probably best to draw a sketch here.
1853  *
1854  * @code
1855  * const unsigned int dim = 2;
1856  *
1857  * const std::vector<Point<2>> vertices = {
1858  * {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0},
1859  * {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5},
1860  * {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0},
1861  * {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5},
1862  * {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}};
1863  *
1864  * @endcode
1865  *
1866  * Next, we have to define the cells and the vertices they contain.
1867  *
1868  * @code
1869  * const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
1870  * cell_vertices = {{{0, 1, 5, 6}},
1871  * {{1, 2, 6, 7}},
1872  * {{2, 3, 7, 8}},
1873  * {{3, 4, 8, 9}},
1874  * {{5, 6, 10, 11}},
1875  * {{8, 9, 12, 13}},
1876  * {{10, 11, 14, 15}},
1877  * {{12, 13, 17, 18}},
1878  * {{14, 15, 19, 20}},
1879  * {{15, 16, 20, 21}},
1880  * {{16, 17, 21, 22}},
1881  * {{17, 18, 22, 23}}};
1882  *
1883  * const unsigned int n_cells = cell_vertices.size();
1884  *
1885  * @endcode
1886  *
1887  * Again, we generate a C++ vector type from this, but this time by
1888  * looping over the cells (yes, this is boring). Additionally, we set
1889  * the material indicator to zero for all the cells:
1890  *
1891  * @code
1892  * std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
1893  * for (unsigned int i = 0; i < n_cells; ++i)
1894  * {
1895  * for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell;
1896  * ++j)
1897  * cells[i].vertices[j] = cell_vertices[i][j];
1898  * cells[i].material_id = 0;
1899  * }
1900  *
1901  * @endcode
1902  *
1903  * Finally pass all this information to the library to generate a
1904  * triangulation. The last parameter may be used to pass information
1905  * about non-zero boundary indicators at certain faces of the
1906  * triangulation to the library, but we don't want that here, so we give
1907  * an empty object:
1908  *
1909  * @code
1910  * coarse_grid.create_triangulation(vertices, cells, SubCellData());
1911  *
1912  * @endcode
1913  *
1914  * And since we want that the evaluation point (3/4,3/4) in this example
1915  * is a grid point, we refine once globally:
1916  *
1917  * @code
1918  * coarse_grid.refine_global(1);
1919  * }
1920  * } // namespace Data
1921  *
1922  * @endcode
1923  *
1924  *
1925  * <a name="Discussion"></a>
1926  * <h4>Discussion</h4>
1927  *
1928 
1929  *
1930  * As you have now read through this framework, you may be wondering why we
1931  * have not chosen to implement the classes implementing a certain setup
1932  * (like the <code>CurvedRidges</code> class) directly as classes derived
1933  * from <code>Data::SetUpBase</code>. Indeed, we could have done very well
1934  * so. The only reason is that then we would have to have member variables
1935  * for the solution and right hand side classes in the
1936  * <code>CurvedRidges</code> class, as well as member functions overloading
1937  * the abstract functions of the base class giving access to these member
1938  * variables. The <code>SetUp</code> class has the sole reason to relieve us
1939  * from the need to reiterate these member variables and functions that
1940  * would be necessary in all such classes. In some way, the template
1941  * mechanism here only provides a way to have default implementations for a
1942  * number of functions that depend on external quantities and can thus not
1943  * be provided using normal virtual functions, at least not without the help
1944  * of templates.
1945  *
1946 
1947  *
1948  * However, there might be good reasons to actually implement classes
1949  * derived from <code>Data::SetUpBase</code>, for example if the solution or
1950  * right hand side classes require constructors that take arguments, which
1951  * the <code>Data::SetUpBase</code> class cannot provide. In that case,
1952  * subclassing is a worthwhile strategy. Other possibilities for special
1953  * cases are to derive from <code>Data::SetUp@<SomeSetUp@></code> where
1954  * <code>SomeSetUp</code> denotes a class, or even to explicitly specialize
1955  * <code>Data::SetUp@<SomeSetUp@></code>. The latter allows to transparently
1956  * use the way the <code>SetUp</code> class is used for other set-ups, but
1957  * with special actions taken for special arguments.
1958  *
1959 
1960  *
1961  * A final observation favoring the approach taken here is the following: we
1962  * have found numerous times that when starting a project, the number of
1963  * parameters (usually boundary values, right hand side, coarse grid, just
1964  * as here) was small, and the number of test cases was small as well. One
1965  * then starts out by handcoding them into a number of <code>switch</code>
1966  * statements. Over time, projects grow, and so does the number of test
1967  * cases. The number of <code>switch</code> statements grows with that, and
1968  * their length as well, and one starts to find ways to consider impossible
1969  * examples where domains, boundary values, and right hand sides do not fit
1970  * together any more, and starts losing the overview over the whole
1971  * structure. Encapsulating everything belonging to a certain test case into
1972  * a structure of its own has proven worthwhile for this, as it keeps
1973  * everything that belongs to one test case in one place. Furthermore, it
1974  * allows to put these things all in one or more files that are only devoted
1975  * to test cases and their data, without having to bring their actual
1976  * implementation into contact with the rest of the program.
1977  *
1978 
1979  *
1980  *
1981 
1982  *
1983  *
1984  * <a name="Dualfunctionals"></a>
1985  * <h3>Dual functionals</h3>
1986  *
1987 
1988  *
1989  * As with the other components of the program, we put everything we need to
1990  * describe dual functionals into a namespace of its own, and define an
1991  * abstract base class that provides the interface the class solving the
1992  * dual problem needs for its work.
1993  *
1994 
1995  *
1996  * We will then implement two such classes, for the evaluation of a point
1997  * value and of the derivative of the solution at that point. For these
1998  * functionals we already have the corresponding evaluation objects, so they
1999  * are complementary.
2000  *
2001  * @code
2002  * namespace DualFunctional
2003  * {
2004  * @endcode
2005  *
2006  *
2007  * <a name="TheDualFunctionalBaseclass"></a>
2008  * <h4>The DualFunctionalBase class</h4>
2009  *
2010 
2011  *
2012  * First start with the base class for dual functionals. Since for linear
2013  * problems the characteristics of the dual problem play a role only in
2014  * the right hand side, we only need to provide for a function that
2015  * assembles the right hand side for a given discretization:
2016  *
2017  * @code
2018  * template <int dim>
2019  * class DualFunctionalBase : public Subscriptor
2020  * {
2021  * public:
2022  * virtual void assemble_rhs(const DoFHandler<dim> &dof_handler,
2023  * Vector<double> & rhs) const = 0;
2024  * };
2025  *
2026  *
2027  * @endcode
2028  *
2029  *
2030  * <a name="ThedualfunctionalPointValueEvaluationclass"></a>
2031  * <h4>The dual functional PointValueEvaluation class</h4>
2032  *
2033 
2034  *
2035  * As a first application, we consider the functional corresponding to the
2036  * evaluation of the solution's value at a given point which again we
2037  * assume to be a vertex. Apart from the constructor that takes and stores
2038  * the evaluation point, this class consists only of the function that
2039  * implements assembling the right hand side.
2040  *
2041  * @code
2042  * template <int dim>
2043  * class PointValueEvaluation : public DualFunctionalBase<dim>
2044  * {
2045  * public:
2046  * PointValueEvaluation(const Point<dim> &evaluation_point);
2047  *
2048  * virtual void assemble_rhs(const DoFHandler<dim> &dof_handler,
2049  * Vector<double> & rhs) const override;
2050  *
2051  * DeclException1(
2052  * ExcEvaluationPointNotFound,
2053  * Point<dim>,
2054  * << "The evaluation point " << arg1
2055  * << " was not found among the vertices of the present grid.");
2056  *
2057  * protected:
2058  * const Point<dim> evaluation_point;
2059  * };
2060  *
2061  *
2062  * template <int dim>
2063  * PointValueEvaluation<dim>::PointValueEvaluation(
2064  * const Point<dim> &evaluation_point)
2065  * : evaluation_point(evaluation_point)
2066  * {}
2067  *
2068  *
2069  * @endcode
2070  *
2071  * As for doing the main purpose of the class, assembling the right hand
2072  * side, let us first consider what is necessary: The right hand side of
2073  * the dual problem is a vector of values J(phi_i), where J is the error
2074  * functional, and phi_i is the i-th shape function. Here, J is the
2075  * evaluation at the point x0, i.e. J(phi_i)=phi_i(x0).
2076  *
2077 
2078  *
2079  * Now, we have assumed that the evaluation point is a vertex. Thus, for
2080  * the usual finite elements we might be using in this program, we can
2081  * take for granted that at such a point exactly one shape function is
2082  * nonzero, and in particular has the value one. Thus, we set the right
2083  * hand side vector to all-zeros, then seek for the shape function
2084  * associated with that point and set the corresponding value of the right
2085  * hand side vector to one:
2086  *
2087  * @code
2088  * template <int dim>
2089  * void
2090  * PointValueEvaluation<dim>::assemble_rhs(const DoFHandler<dim> &dof_handler,
2091  * Vector<double> & rhs) const
2092  * {
2093  * @endcode
2094  *
2095  * So, first set everything to zeros...
2096  *
2097  * @code
2098  * rhs.reinit(dof_handler.n_dofs());
2099  *
2100  * @endcode
2101  *
2102  * ...then loop over cells and find the evaluation point among the
2103  * vertices (or very close to a vertex, which may happen due to floating
2104  * point round-off):
2105  *
2106  * @code
2107  * for (const auto &cell : dof_handler.active_cell_iterators())
2108  * for (unsigned int vertex = 0;
2109  * vertex < GeometryInfo<dim>::vertices_per_cell;
2110  * ++vertex)
2111  * if (cell->vertex(vertex).distance(evaluation_point) <
2112  * cell->diameter() * 1e-8)
2113  * {
2114  * @endcode
2115  *
2116  * Ok, found, so set corresponding entry, and leave function
2117  * since we are finished:
2118  *
2119  * @code
2120  * rhs(cell->vertex_dof_index(vertex, 0)) = 1;
2121  * return;
2122  * }
2123  *
2124  * @endcode
2125  *
2126  * Finally, a sanity check: if we somehow got here, then we must have
2127  * missed the evaluation point, so raise an exception unconditionally:
2128  *
2129  * @code
2130  * AssertThrow(false, ExcEvaluationPointNotFound(evaluation_point));
2131  * }
2132  *
2133  *
2134  * @endcode
2135  *
2136  *
2137  * <a name="ThedualfunctionalPointXDerivativeEvaluationclass"></a>
2138  * <h4>The dual functional PointXDerivativeEvaluation class</h4>
2139  *
2140 
2141  *
2142  * As second application, we again consider the evaluation of the
2143  * x-derivative of the solution at one point. Again, the declaration of
2144  * the class, and the implementation of its constructor is not too
2145  * interesting:
2146  *
2147  * @code
2148  * template <int dim>
2149  * class PointXDerivativeEvaluation : public DualFunctionalBase<dim>
2150  * {
2151  * public:
2152  * PointXDerivativeEvaluation(const Point<dim> &evaluation_point);
2153  *
2154  * virtual void assemble_rhs(const DoFHandler<dim> &dof_handler,
2155  * Vector<double> & rhs) const;
2156  *
2157  * DeclException1(
2158  * ExcEvaluationPointNotFound,
2159  * Point<dim>,
2160  * << "The evaluation point " << arg1
2161  * << " was not found among the vertices of the present grid.");
2162  *
2163  * protected:
2164  * const Point<dim> evaluation_point;
2165  * };
2166  *
2167  *
2168  * template <int dim>
2169  * PointXDerivativeEvaluation<dim>::PointXDerivativeEvaluation(
2170  * const Point<dim> &evaluation_point)
2171  * : evaluation_point(evaluation_point)
2172  * {}
2173  *
2174  *
2175  * @endcode
2176  *
2177  * What is interesting is the implementation of this functional: here,
2178  * J(phi_i)=d/dx phi_i(x0).
2179  *
2180 
2181  *
2182  * We could, as in the implementation of the respective evaluation object
2183  * take the average of the gradients of each shape function phi_i at this
2184  * evaluation point. However, we take a slightly different approach: we
2185  * simply take the average over all cells that surround this point. The
2186  * question which cells <code>surrounds</code> the evaluation point is
2187  * made dependent on the mesh width by including those cells for which the
2188  * distance of the cell's midpoint to the evaluation point is less than
2189  * the cell's diameter.
2190  *
2191 
2192  *
2193  * Taking the average of the gradient over the area/volume of these cells
2194  * leads to a dual solution which is very close to the one which would
2195  * result from the point evaluation of the gradient. It is simple to
2196  * justify theoretically that this does not change the method
2197  * significantly.
2198  *
2199  * @code
2200  * template <int dim>
2201  * void PointXDerivativeEvaluation<dim>::assemble_rhs(
2202  * const DoFHandler<dim> &dof_handler,
2203  * Vector<double> & rhs) const
2204  * {
2205  * @endcode
2206  *
2207  * Again, first set all entries to zero:
2208  *
2209  * @code
2210  * rhs.reinit(dof_handler.n_dofs());
2211  *
2212  * @endcode
2213  *
2214  * Initialize a <code>FEValues</code> object with a quadrature formula,
2215  * have abbreviations for the number of quadrature points and shape
2216  * functions...
2217  *
2218  * @code
2219  * QGauss<dim> quadrature(dof_handler.get_fe().degree + 1);
2220  * FEValues<dim> fe_values(dof_handler.get_fe(),
2221  * quadrature,
2222  * update_gradients | update_quadrature_points |
2223  * update_JxW_values);
2224  * const unsigned int n_q_points = fe_values.n_quadrature_points;
2225  * const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2226  *
2227  * @endcode
2228  *
2229  * ...and have two objects that are used to store the global indices of
2230  * the degrees of freedom on a cell, and the values of the gradients of
2231  * the shape functions at the quadrature points:
2232  *
2233  * @code
2234  * Vector<double> cell_rhs(dofs_per_cell);
2235  * std::vector<unsigned int> local_dof_indices(dofs_per_cell);
2236  *
2237  * @endcode
2238  *
2239  * Finally have a variable in which we will sum up the area/volume of
2240  * the cells over which we integrate, by integrating the unit functions
2241  * on these cells:
2242  *
2243  * @code
2244  * double total_volume = 0;
2245  *
2246  * @endcode
2247  *
2248  * Then start the loop over all cells, and select those cells which are
2249  * close enough to the evaluation point:
2250  *
2251  * @code
2252  * for (const auto &cell : dof_handler.active_cell_iterators())
2253  * if (cell->center().distance(evaluation_point) <= cell->diameter())
2254  * {
2255  * @endcode
2256  *
2257  * If we have found such a cell, then initialize the
2258  * <code>FEValues</code> object and integrate the x-component of
2259  * the gradient of each shape function, as well as the unit
2260  * function for the total area/volume.
2261  *
2262  * @code
2263  * fe_values.reinit(cell);
2264  * cell_rhs = 0;
2265  *
2266  * for (unsigned int q = 0; q < n_q_points; ++q)
2267  * {
2268  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2269  * cell_rhs(i) +=
2270  * fe_values.shape_grad(i, q)[0] // (d/dx phi_i(x_q))
2271  * * fe_values.JxW(q); // * dx
2272  * total_volume += fe_values.JxW(q);
2273  * }
2274  *
2275  * @endcode
2276  *
2277  * If we have the local contributions, distribute them to the
2278  * global vector:
2279  *
2280  * @code
2281  * cell->get_dof_indices(local_dof_indices);
2282  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2283  * rhs(local_dof_indices[i]) += cell_rhs(i);
2284  * }
2285  *
2286  * @endcode
2287  *
2288  * After we have looped over all cells, check whether we have found any
2289  * at all, by making sure that their volume is non-zero. If not, then
2290  * the results will be botched, as the right hand side should then still
2291  * be zero, so throw an exception:
2292  *
2293  * @code
2294  * AssertThrow(total_volume > 0,
2295  * ExcEvaluationPointNotFound(evaluation_point));
2296  *
2297  * @endcode
2298  *
2299  * Finally, we have by now only integrated the gradients of the shape
2300  * functions, not taking their mean value. We fix this by dividing by
2301  * the measure of the volume over which we have integrated:
2302  *
2303  * @code
2304  * rhs /= total_volume;
2305  * }
2306  *
2307  *
2308  * } // namespace DualFunctional
2309  *
2310  *
2311  * @endcode
2312  *
2313  *
2314  * <a name="ExtendingtheLaplaceSolvernamespace"></a>
2315  * <h3>Extending the LaplaceSolver namespace</h3>
2316  *
2317  * @code
2318  * namespace LaplaceSolver
2319  * {
2320  * @endcode
2321  *
2322  *
2323  * <a name="TheDualSolverclass"></a>
2324  * <h4>The DualSolver class</h4>
2325  *
2326 
2327  *
2328  * In the same way as the <code>PrimalSolver</code> class above, we now
2329  * implement a <code>DualSolver</code>. It has all the same features, the
2330  * only difference is that it does not take a function object denoting a
2331  * right hand side object, but now takes a <code>DualFunctionalBase</code>
2332  * object that will assemble the right hand side vector of the dual
2333  * problem. The rest of the class is rather trivial.
2334  *
2335 
2336  *
2337  * Since both primal and dual solver will use the same triangulation, but
2338  * different discretizations, it now becomes clear why we have made the
2339  * <code>Base</code> class a virtual one: since the final class will be
2340  * derived from both <code>PrimalSolver</code> as well as
2341  * <code>DualSolver</code>, it would have two <code>Base</code> instances,
2342  * would we not have marked the inheritance as virtual. Since in many
2343  * applications the base class would store much more information than just
2344  * the triangulation which needs to be shared between primal and dual
2345  * solvers, we do not usually want to use two such base classes.
2346  *
2347  * @code
2348  * template <int dim>
2349  * class DualSolver : public Solver<dim>
2350  * {
2351  * public:
2352  * DualSolver(
2353  * Triangulation<dim> & triangulation,
2354  * const FiniteElement<dim> & fe,
2355  * const Quadrature<dim> & quadrature,
2356  * const Quadrature<dim - 1> & face_quadrature,
2357  * const DualFunctional::DualFunctionalBase<dim> &dual_functional);
2358  *
2359  * protected:
2360  * const SmartPointer<const DualFunctional::DualFunctionalBase<dim>>
2361  * dual_functional;
2362  * virtual void assemble_rhs(Vector<double> &rhs) const override;
2363  *
2364  * static const Functions::ZeroFunction<dim> boundary_values;
2365  * };
2366  *
2367  * template <int dim>
2368  * const Functions::ZeroFunction<dim> DualSolver<dim>::boundary_values;
2369  *
2370  * template <int dim>
2371  * DualSolver<dim>::DualSolver(
2372  * Triangulation<dim> & triangulation,
2373  * const FiniteElement<dim> & fe,
2374  * const Quadrature<dim> & quadrature,
2375  * const Quadrature<dim - 1> & face_quadrature,
2376  * const DualFunctional::DualFunctionalBase<dim> &dual_functional)
2377  * : Base<dim>(triangulation)
2378  * , Solver<dim>(triangulation,
2379  * fe,
2380  * quadrature,
2381  * face_quadrature,
2382  * boundary_values)
2383  * , dual_functional(&dual_functional)
2384  * {}
2385  *
2386  *
2387  *
2388  * template <int dim>
2389  * void DualSolver<dim>::assemble_rhs(Vector<double> &rhs) const
2390  * {
2391  * dual_functional->assemble_rhs(this->dof_handler, rhs);
2392  * }
2393  *
2394  *
2395  * @endcode
2396  *
2397  *
2398  * <a name="TheWeightedResidualclass"></a>
2399  * <h4>The WeightedResidual class</h4>
2400  *
2401 
2402  *
2403  * Here finally comes the main class of this program, the one that
2404  * implements the dual weighted residual error estimator. It joins the
2405  * primal and dual solver classes to use them for the computation of
2406  * primal and dual solutions, and implements the error representation
2407  * formula for use as error estimate and mesh refinement.
2408  *
2409 
2410  *
2411  * The first few of the functions of this class are mostly overriders of
2412  * the respective functions of the base class:
2413  *
2414  * @code
2415  * template <int dim>
2416  * class WeightedResidual : public PrimalSolver<dim>, public DualSolver<dim>
2417  * {
2418  * public:
2419  * WeightedResidual(
2420  * Triangulation<dim> & coarse_grid,
2421  * const FiniteElement<dim> & primal_fe,
2422  * const FiniteElement<dim> & dual_fe,
2423  * const Quadrature<dim> & quadrature,
2424  * const Quadrature<dim - 1> & face_quadrature,
2425  * const Function<dim> & rhs_function,
2426  * const Function<dim> & boundary_values,
2427  * const DualFunctional::DualFunctionalBase<dim> &dual_functional);
2428  *
2429  * virtual void solve_problem() override;
2430  *
2431  * virtual void postprocess(
2432  * const Evaluation::EvaluationBase<dim> &postprocessor) const override;
2433  *
2434  * virtual unsigned int n_dofs() const override;
2435  *
2436  * virtual void refine_grid() override;
2437  *
2438  * virtual void output_solution() const override;
2439  *
2440  * private:
2441  * @endcode
2442  *
2443  * In the private section, we have two functions that are used to call
2444  * the <code>solve_problem</code> functions of the primal and dual base
2445  * classes. These two functions will be called in parallel by the
2446  * <code>solve_problem</code> function of this class.
2447  *
2448  * @code
2449  * void solve_primal_problem();
2450  * void solve_dual_problem();
2451  * @endcode
2452  *
2453  * Then declare abbreviations for active cell iterators, to avoid that
2454  * we have to write this lengthy name over and over again:
2455  *
2456 
2457  *
2458  *
2459  * @code
2460  * using active_cell_iterator =
2461  * typename DoFHandler<dim>::active_cell_iterator;
2462  *
2463  * @endcode
2464  *
2465  * Next, declare a data type that we will us to store the contribution
2466  * of faces to the error estimator. The idea is that we can compute the
2467  * face terms from each of the two cells to this face, as they are the
2468  * same when viewed from both sides. What we will do is to compute them
2469  * only once, based on some rules explained below which of the two
2470  * adjacent cells will be in charge to do so. We then store the
2471  * contribution of each face in a map mapping faces to their values, and
2472  * only collect the contributions for each cell by looping over the
2473  * cells a second time and grabbing the values from the map.
2474  *
2475 
2476  *
2477  * The data type of this map is declared here:
2478  *
2479  * @code
2480  * using FaceIntegrals =
2481  * typename std::map<typename DoFHandler<dim>::face_iterator, double>;
2482  *
2483  * @endcode
2484  *
2485  * In the computation of the error estimates on cells and faces, we need
2486  * a number of helper objects, such as <code>FEValues</code> and
2487  * <code>FEFaceValues</code> functions, but also temporary objects
2488  * storing the values and gradients of primal and dual solutions, for
2489  * example. These fields are needed in the three functions that do the
2490  * integration on cells, and regular and irregular faces, respectively.
2491  *
2492 
2493  *
2494  * There are three reasonable ways to provide these fields: first, as
2495  * local variables in the function that needs them; second, as member
2496  * variables of this class; third, as arguments passed to that function.
2497  *
2498 
2499  *
2500  * These three alternatives all have drawbacks: the third that their
2501  * number is not negligible and would make calling these functions a
2502  * lengthy enterprise. The second has the drawback that it disallows
2503  * parallelization, since the threads that will compute the error
2504  * estimate have to have their own copies of these variables each, so
2505  * member variables of the enclosing class will not work. The first
2506  * approach, although straightforward, has a subtle but important
2507  * drawback: we will call these functions over and over again, many
2508  * thousands of times maybe; it now turns out that allocating
2509  * vectors and other objects that need memory from the heap is an
2510  * expensive business in terms of run-time, since memory allocation is
2511  * expensive when several threads are involved. It is thus
2512  * significantly better to allocate the memory only once, and recycle
2513  * the objects as often as possible.
2514  *
2515 
2516  *
2517  * What to do? Our answer is to use a variant of the third strategy.
2518  * In fact, this is exactly what the WorkStream concept is supposed to
2519  * do (we have already introduced it above, but see also @ref threads).
2520  * To avoid that we have to give these functions a dozen or so
2521  * arguments, we pack all these variables into two structures, one which
2522  * is used for the computations on cells, the other doing them on the
2523  * faces. Both are then joined into the WeightedResidualScratchData class
2524  * that will serve as the "scratch data" class of the WorkStream concept:
2525  *
2526  * @code
2527  * struct CellData
2528  * {
2529  * FEValues<dim> fe_values;
2530  * const SmartPointer<const Function<dim>> right_hand_side;
2531  *
2532  * std::vector<double> cell_residual;
2533  * std::vector<double> rhs_values;
2534  * std::vector<double> dual_weights;
2535  * std::vector<double> cell_laplacians;
2536  * CellData(const FiniteElement<dim> &fe,
2537  * const Quadrature<dim> & quadrature,
2538  * const Function<dim> & right_hand_side);
2539  * CellData(const CellData &cell_data);
2540  * };
2541  *
2542  * struct FaceData
2543  * {
2544  * FEFaceValues<dim> fe_face_values_cell;
2545  * FEFaceValues<dim> fe_face_values_neighbor;
2546  * FESubfaceValues<dim> fe_subface_values_cell;
2547  *
2548  * std::vector<double> jump_residual;
2549  * std::vector<double> dual_weights;
2550  * typename std::vector<Tensor<1, dim>> cell_grads;
2551  * typename std::vector<Tensor<1, dim>> neighbor_grads;
2552  * FaceData(const FiniteElement<dim> & fe,
2553  * const Quadrature<dim - 1> &face_quadrature);
2554  * FaceData(const FaceData &face_data);
2555  * };
2556  *
2557  * struct WeightedResidualScratchData
2558  * {
2559  * WeightedResidualScratchData(
2560  * const FiniteElement<dim> & primal_fe,
2561  * const Quadrature<dim> & primal_quadrature,
2562  * const Quadrature<dim - 1> &primal_face_quadrature,
2563  * const Function<dim> & rhs_function,
2564  * const Vector<double> & primal_solution,
2565  * const Vector<double> & dual_weights);
2566  *
2567  * WeightedResidualScratchData(
2568  * const WeightedResidualScratchData &scratch_data);
2569  *
2570  * CellData cell_data;
2571  * FaceData face_data;
2572  * Vector<double> primal_solution;
2573  * Vector<double> dual_weights;
2574  * };
2575  *
2576  *
2577  * @endcode
2578  *
2579  * WorkStream::run generally wants both a scratch object and a copy
2580  * object. Here, for reasons similar to what we had in @ref step_9 "step-9" when
2581  * discussing the computation of an approximation of the gradient, we
2582  * don't actually need a "copy data" structure. Since WorkStream insists
2583  * on having one of these, we just declare an empty structure that does
2584  * nothing other than being there.
2585  *
2586  * @code
2587  * struct WeightedResidualCopyData
2588  * {};
2589  *
2590  *
2591  *
2592  * @endcode
2593  *
2594  * Regarding the evaluation of the error estimator, we have one driver
2595  * function that uses WorkStream::run() to call the second function on
2596  * every cell:
2597  *
2598  * @code
2599  * void estimate_error(Vector<float> &error_indicators) const;
2600  *
2601  * void estimate_on_one_cell(const active_cell_iterator & cell,
2602  * WeightedResidualScratchData &scratch_data,
2603  * WeightedResidualCopyData & copy_data,
2604  * Vector<float> & error_indicators,
2605  * FaceIntegrals &face_integrals) const;
2606  *
2607  * @endcode
2608  *
2609  * Then we have functions that do the actual integration of the error
2610  * representation formula. They will treat the terms on the cell
2611  * interiors, on those faces that have no hanging nodes, and on those
2612  * faces with hanging nodes, respectively:
2613  *
2614  * @code
2615  * void integrate_over_cell(const active_cell_iterator &cell,
2616  * const Vector<double> & primal_solution,
2617  * const Vector<double> & dual_weights,
2618  * CellData & cell_data,
2619  * Vector<float> &error_indicators) const;
2620  *
2621  * void integrate_over_regular_face(const active_cell_iterator &cell,
2622  * const unsigned int face_no,
2623  * const Vector<double> &primal_solution,
2624  * const Vector<double> &dual_weights,
2625  * FaceData & face_data,
2626  * FaceIntegrals &face_integrals) const;
2627  * void integrate_over_irregular_face(const active_cell_iterator &cell,
2628  * const unsigned int face_no,
2629  * const Vector<double> &primal_solution,
2630  * const Vector<double> &dual_weights,
2631  * FaceData & face_data,
2632  * FaceIntegrals &face_integrals) const;
2633  * };
2634  *
2635  *
2636  *
2637  * @endcode
2638  *
2639  * In the implementation of this class, we first have the constructors of
2640  * the <code>CellData</code> and <code>FaceData</code> member classes, and
2641  * the <code>WeightedResidual</code> constructor. They only initialize
2642  * fields to their correct lengths, so we do not have to discuss them in
2643  * too much detail:
2644  *
2645  * @code
2646  * template <int dim>
2647  * WeightedResidual<dim>::CellData::CellData(
2648  * const FiniteElement<dim> &fe,
2649  * const Quadrature<dim> & quadrature,
2650  * const Function<dim> & right_hand_side)
2651  * : fe_values(fe,
2652  * quadrature,
2655  * , right_hand_side(&right_hand_side)
2656  * , cell_residual(quadrature.size())
2657  * , rhs_values(quadrature.size())
2658  * , dual_weights(quadrature.size())
2659  * , cell_laplacians(quadrature.size())
2660  * {}
2661  *
2662  *
2663  *
2664  * template <int dim>
2665  * WeightedResidual<dim>::CellData::CellData(const CellData &cell_data)
2666  * : fe_values(cell_data.fe_values.get_fe(),
2667  * cell_data.fe_values.get_quadrature(),
2670  * , right_hand_side(cell_data.right_hand_side)
2671  * , cell_residual(cell_data.cell_residual)
2672  * , rhs_values(cell_data.rhs_values)
2673  * , dual_weights(cell_data.dual_weights)
2674  * , cell_laplacians(cell_data.cell_laplacians)
2675  * {}
2676  *
2677  *
2678  *
2679  * template <int dim>
2680  * WeightedResidual<dim>::FaceData::FaceData(
2681  * const FiniteElement<dim> & fe,
2682  * const Quadrature<dim - 1> &face_quadrature)
2683  * : fe_face_values_cell(fe,
2684  * face_quadrature,
2687  * , fe_face_values_neighbor(fe,
2688  * face_quadrature,
2691  * , fe_subface_values_cell(fe, face_quadrature, update_gradients)
2692  * {
2693  * const unsigned int n_face_q_points = face_quadrature.size();
2694  *
2695  * jump_residual.resize(n_face_q_points);
2696  * dual_weights.resize(n_face_q_points);
2697  * cell_grads.resize(n_face_q_points);
2698  * neighbor_grads.resize(n_face_q_points);
2699  * }
2700  *
2701  *
2702  *
2703  * template <int dim>
2704  * WeightedResidual<dim>::FaceData::FaceData(const FaceData &face_data)
2705  * : fe_face_values_cell(face_data.fe_face_values_cell.get_fe(),
2706  * face_data.fe_face_values_cell.get_quadrature(),
2709  * , fe_face_values_neighbor(
2710  * face_data.fe_face_values_neighbor.get_fe(),
2711  * face_data.fe_face_values_neighbor.get_quadrature(),
2714  * , fe_subface_values_cell(
2715  * face_data.fe_subface_values_cell.get_fe(),
2716  * face_data.fe_subface_values_cell.get_quadrature(),
2717  * update_gradients)
2718  * , jump_residual(face_data.jump_residual)
2719  * , dual_weights(face_data.dual_weights)
2720  * , cell_grads(face_data.cell_grads)
2721  * , neighbor_grads(face_data.neighbor_grads)
2722  * {}
2723  *
2724  *
2725  *
2726  * template <int dim>
2727  * WeightedResidual<dim>::WeightedResidualScratchData::
2728  * WeightedResidualScratchData(
2729  * const FiniteElement<dim> & primal_fe,
2730  * const Quadrature<dim> & primal_quadrature,
2731  * const Quadrature<dim - 1> &primal_face_quadrature,
2732  * const Function<dim> & rhs_function,
2733  * const Vector<double> & primal_solution,
2734  * const Vector<double> & dual_weights)
2735  * : cell_data(primal_fe, primal_quadrature, rhs_function)
2736  * , face_data(primal_fe, primal_face_quadrature)
2737  * , primal_solution(primal_solution)
2738  * , dual_weights(dual_weights)
2739  * {}
2740  *
2741  * template <int dim>
2742  * WeightedResidual<dim>::WeightedResidualScratchData::
2743  * WeightedResidualScratchData(
2744  * const WeightedResidualScratchData &scratch_data)
2745  * : cell_data(scratch_data.cell_data)
2746  * , face_data(scratch_data.face_data)
2747  * , primal_solution(scratch_data.primal_solution)
2748  * , dual_weights(scratch_data.dual_weights)
2749  * {}
2750  *
2751  *
2752  *
2753  * template <int dim>
2754  * WeightedResidual<dim>::WeightedResidual(
2755  * Triangulation<dim> & coarse_grid,
2756  * const FiniteElement<dim> & primal_fe,
2757  * const FiniteElement<dim> & dual_fe,
2758  * const Quadrature<dim> & quadrature,
2759  * const Quadrature<dim - 1> & face_quadrature,
2760  * const Function<dim> & rhs_function,
2761  * const Function<dim> & bv,
2762  * const DualFunctional::DualFunctionalBase<dim> &dual_functional)
2763  * : Base<dim>(coarse_grid)
2764  * , PrimalSolver<dim>(coarse_grid,
2765  * primal_fe,
2766  * quadrature,
2767  * face_quadrature,
2768  * rhs_function,
2769  * bv)
2770  * , DualSolver<dim>(coarse_grid,
2771  * dual_fe,
2772  * quadrature,
2773  * face_quadrature,
2774  * dual_functional)
2775  * {}
2776  *
2777  *
2778  * @endcode
2779  *
2780  * The next five functions are boring, as they simply relay their work to
2781  * the base classes. The first calls the primal and dual solvers in
2782  * parallel, while postprocessing the solution and retrieving the number
2783  * of degrees of freedom is done by the primal class.
2784  *
2785  * @code
2786  * template <int dim>
2787  * void WeightedResidual<dim>::solve_problem()
2788  * {
2789  * Threads::TaskGroup<void> tasks;
2790  * tasks +=
2791  * Threads::new_task(&WeightedResidual<dim>::solve_primal_problem, *this);
2792  * tasks +=
2793  * Threads::new_task(&WeightedResidual<dim>::solve_dual_problem, *this);
2794  * tasks.join_all();
2795  * }
2796  *
2797  *
2798  * template <int dim>
2799  * void WeightedResidual<dim>::solve_primal_problem()
2800  * {
2801  * PrimalSolver<dim>::solve_problem();
2802  * }
2803  *
2804  * template <int dim>
2805  * void WeightedResidual<dim>::solve_dual_problem()
2806  * {
2807  * DualSolver<dim>::solve_problem();
2808  * }
2809  *
2810  *
2811  * template <int dim>
2812  * void WeightedResidual<dim>::postprocess(
2813  * const Evaluation::EvaluationBase<dim> &postprocessor) const
2814  * {
2815  * PrimalSolver<dim>::postprocess(postprocessor);
2816  * }
2817  *
2818  *
2819  * template <int dim>
2820  * unsigned int WeightedResidual<dim>::n_dofs() const
2821  * {
2822  * return PrimalSolver<dim>::n_dofs();
2823  * }
2824  *
2825  *
2826  *
2827  * @endcode
2828  *
2829  * Now, it is becoming more interesting: the <code>refine_grid()</code>
2830  * function asks the error estimator to compute the cell-wise error
2831  * indicators, then uses their absolute values for mesh refinement.
2832  *
2833  * @code
2834  * template <int dim>
2835  * void WeightedResidual<dim>::refine_grid()
2836  * {
2837  * @endcode
2838  *
2839  * First call the function that computes the cell-wise and global error:
2840  *
2841  * @code
2842  * Vector<float> error_indicators(this->triangulation->n_active_cells());
2843  * estimate_error(error_indicators);
2844  *
2845  * @endcode
2846  *
2847  * Then note that marking cells for refinement or coarsening only works
2848  * if all indicators are positive, to allow their comparison. Thus, drop
2849  * the signs on all these indicators:
2850  *
2851  * @code
2852  * for (float &error_indicator : error_indicators)
2853  * error_indicator = std::fabs(error_indicator);
2854  *
2855  * @endcode
2856  *
2857  * Finally, we can select between different strategies for
2858  * refinement. The default here is to refine those cells with the
2859  * largest error indicators that make up for a total of 80 per cent of
2860  * the error, while we coarsen those with the smallest indicators that
2861  * make up for the bottom 2 per cent of the error.
2862  *
2863  * @code
2865  * error_indicators,
2866  * 0.8,
2867  * 0.02);
2868  * this->triangulation->execute_coarsening_and_refinement();
2869  * }
2870  *
2871  *
2872  * @endcode
2873  *
2874  * Since we want to output both the primal and the dual solution, we
2875  * overload the <code>output_solution</code> function. The only
2876  * interesting feature of this function is that the primal and dual
2877  * solutions are defined on different finite element spaces, which is not
2878  * the format the <code>DataOut</code> class expects. Thus, we have to
2879  * transfer them to a common finite element space. Since we want the
2880  * solutions only to see them qualitatively, we contend ourselves with
2881  * interpolating the dual solution to the (smaller) primal space. For the
2882  * interpolation, there is a library function, that takes a
2883  * AffineConstraints object including the hanging node
2884  * constraints. The rest is standard.
2885  *
2886  * @code
2887  * template <int dim>
2888  * void WeightedResidual<dim>::output_solution() const
2889  * {
2890  * AffineConstraints<double> primal_hanging_node_constraints;
2891  * DoFTools::make_hanging_node_constraints(PrimalSolver<dim>::dof_handler,
2892  * primal_hanging_node_constraints);
2893  * primal_hanging_node_constraints.close();
2894  * Vector<double> dual_solution(PrimalSolver<dim>::dof_handler.n_dofs());
2895  * FETools::interpolate(DualSolver<dim>::dof_handler,
2896  * DualSolver<dim>::solution,
2897  * PrimalSolver<dim>::dof_handler,
2898  * primal_hanging_node_constraints,
2899  * dual_solution);
2900  *
2901  * DataOut<dim> data_out;
2902  * data_out.attach_dof_handler(PrimalSolver<dim>::dof_handler);
2903  *
2904  * @endcode
2905  *
2906  * Add the data vectors for which we want output. Add them both, the
2907  * <code>DataOut</code> functions can handle as many data vectors as you
2908  * wish to write to output:
2909  *
2910  * @code
2911  * data_out.add_data_vector(PrimalSolver<dim>::solution, "primal_solution");
2912  * data_out.add_data_vector(dual_solution, "dual_solution");
2913  *
2914  * data_out.build_patches();
2915  *
2916  * std::ofstream out("solution-" + std::to_string(this->refinement_cycle) +
2917  * ".vtu");
2918  * data_out.write(out, DataOutBase::vtu);
2919  * }
2920  *
2921  *
2922  * @endcode
2923  *
2924  *
2925  * <a name="Estimatingerrors"></a>
2926  * <h3>Estimating errors</h3>
2927  *
2928 
2929  *
2930  *
2931  * <a name="Errorestimationdriverfunctions"></a>
2932  * <h4>Error estimation driver functions</h4>
2933  *
2934 
2935  *
2936  * As for the actual computation of error estimates, let's start with the
2937  * function that drives all this, i.e. calls those functions that actually
2938  * do the work, and finally collects the results.
2939  *
2940  * @code
2941  * template <int dim>
2942  * void
2943  * WeightedResidual<dim>::estimate_error(Vector<float> &error_indicators) const
2944  * {
2945  * @endcode
2946  *
2947  * The first task in computing the error is to set up vectors that
2948  * denote the primal solution, and the weights (z-z_h)=(z-I_hz), both in
2949  * the finite element space for which we have computed the dual
2950  * solution. For this, we have to interpolate the primal solution to the
2951  * dual finite element space, and to subtract the interpolation of the
2952  * computed dual solution to the primal finite element
2953  * space. Fortunately, the library provides functions for the
2954  * interpolation into larger or smaller finite element spaces, so this
2955  * is mostly obvious.
2956  *
2957 
2958  *
2959  * First, let's do that for the primal solution: it is cell-wise
2960  * interpolated into the finite element space in which we have solved
2961  * the dual problem: But, again as in the
2962  * <code>WeightedResidual::output_solution</code> function we first need
2963  * to create an AffineConstraints object including the hanging node
2964  * constraints, but this time of the dual finite element space.
2965  *
2966  * @code
2967  * AffineConstraints<double> dual_hanging_node_constraints;
2968  * DoFTools::make_hanging_node_constraints(DualSolver<dim>::dof_handler,
2969  * dual_hanging_node_constraints);
2970  * dual_hanging_node_constraints.close();
2971  * Vector<double> primal_solution(DualSolver<dim>::dof_handler.n_dofs());
2972  * FETools::interpolate(PrimalSolver<dim>::dof_handler,
2973  * PrimalSolver<dim>::solution,
2974  * DualSolver<dim>::dof_handler,
2975  * dual_hanging_node_constraints,
2976  * primal_solution);
2977  *
2978  * @endcode
2979  *
2980  * Then for computing the interpolation of the numerically approximated
2981  * dual solution z into the finite element space of the primal solution
2982  * and subtracting it from z: use the
2983  * <code>interpolate_difference</code> function, that gives (z-I_hz) in
2984  * the element space of the dual solution.
2985  *
2986  * @code
2987  * AffineConstraints<double> primal_hanging_node_constraints;
2988  * DoFTools::make_hanging_node_constraints(PrimalSolver<dim>::dof_handler,
2989  * primal_hanging_node_constraints);
2990  * primal_hanging_node_constraints.close();
2991  * Vector<double> dual_weights(DualSolver<dim>::dof_handler.n_dofs());
2992  * FETools::interpolation_difference(DualSolver<dim>::dof_handler,
2993  * dual_hanging_node_constraints,
2994  * DualSolver<dim>::solution,
2995  * PrimalSolver<dim>::dof_handler,
2996  * primal_hanging_node_constraints,
2997  * dual_weights);
2998  *
2999  * @endcode
3000  *
3001  * Note that this could probably have been more efficient since those
3002  * constraints have been used previously when assembling matrix and
3003  * right hand side for the primal problem and writing out the dual
3004  * solution. We leave the optimization of the program in this respect as
3005  * an exercise.
3006  *
3007 
3008  *
3009  * Having computed the dual weights we now proceed with computing the
3010  * cell and face residuals of the primal solution. First we set up a map
3011  * between face iterators and their jump term contributions of faces to
3012  * the error estimator. The reason is that we compute the jump terms
3013  * only once, from one side of the face, and want to collect them only
3014  * afterwards when looping over all cells a second time.
3015  *
3016 
3017  *
3018  * We initialize this map already with a value of -1e20 for all faces,
3019  * since this value will stand out in the results if something should go
3020  * wrong and we fail to compute the value for a face for some
3021  * reason. Secondly, this initialization already makes the std::map
3022  * object allocate all objects it may possibly need. This is important
3023  * since we will write into this structure from parallel threads,
3024  * and doing so would not be thread-safe if the map needed to allocate
3025  * memory and thereby reshape its data structures. In other words, the
3026  * initial initialization relieves us from the necessity to synchronize
3027  * the threads through a mutex each time they write to (and modify the
3028  * structure of) this map.
3029  *
3030  * @code
3031  * FaceIntegrals face_integrals;
3032  * for (const auto &cell :
3033  * DualSolver<dim>::dof_handler.active_cell_iterators())
3034  * for (const auto &face : cell->face_iterators())
3035  * face_integrals[face] = -1e20;
3036  *
3037  * auto worker = [this,
3038  * &error_indicators,
3039  * &face_integrals](const active_cell_iterator & cell,
3040  * WeightedResidualScratchData &scratch_data,
3041  * WeightedResidualCopyData & copy_data) {
3042  * this->estimate_on_one_cell(
3043  * cell, scratch_data, copy_data, error_indicators, face_integrals);
3044  * };
3045  *
3046  * auto do_nothing_copier =
3047  * std::function<void(const WeightedResidualCopyData &)>();
3048  *
3049  * @endcode
3050  *
3051  * Then hand it all off to WorkStream::run() to compute the
3052  * estimators for all cells in parallel:
3053  *
3054  * @code
3055  * WorkStream::run(
3056  * DualSolver<dim>::dof_handler.begin_active(),
3057  * DualSolver<dim>::dof_handler.end(),
3058  * worker,
3059  * do_nothing_copier,
3060  * WeightedResidualScratchData(*DualSolver<dim>::fe,
3061  * *DualSolver<dim>::quadrature,
3062  * *DualSolver<dim>::face_quadrature,
3063  * *this->rhs_function,
3064  * primal_solution,
3065  * dual_weights),
3066  * WeightedResidualCopyData());
3067  *
3068  * @endcode
3069  *
3070  * Once the error contributions are computed, sum them up. For this,
3071  * note that the cell terms are already set, and that only the edge
3072  * terms need to be collected. Thus, loop over all cells and their
3073  * faces, make sure that the contributions of each of the faces are
3074  * there, and add them up. Only take minus one half of the jump term,
3075  * since the other half will be taken by the neighboring cell.
3076  *
3077  * @code
3078  * unsigned int present_cell = 0;
3079  * for (const auto &cell :
3080  * DualSolver<dim>::dof_handler.active_cell_iterators())
3081  * {
3082  * for (const auto &face : cell->face_iterators())
3083  * {
3084  * Assert(face_integrals.find(face) != face_integrals.end(),
3085  * ExcInternalError());
3086  * error_indicators(present_cell) -= 0.5 * face_integrals[face];
3087  * }
3088  * ++present_cell;
3089  * }
3090  * std::cout << " Estimated error="
3091  * << std::accumulate(error_indicators.begin(),
3092  * error_indicators.end(),
3093  * 0.)
3094  * << std::endl;
3095  * }
3096  *
3097  *
3098  * @endcode
3099  *
3100  *
3101  * <a name="Estimatingonasinglecell"></a>
3102  * <h4>Estimating on a single cell</h4>
3103  *
3104 
3105  *
3106  * Next we have the function that is called to estimate the error on a
3107  * single cell. The function may be called multiple times if the library was
3108  * configured to use multithreading. Here it goes:
3109  *
3110  * @code
3111  * template <int dim>
3112  * void WeightedResidual<dim>::estimate_on_one_cell(
3113  * const active_cell_iterator & cell,
3114  * WeightedResidualScratchData &scratch_data,
3115  * WeightedResidualCopyData & copy_data,
3116  * Vector<float> & error_indicators,
3117  * FaceIntegrals & face_integrals) const
3118  * {
3119  * @endcode
3120  *
3121  * Because of WorkStream, estimate_on_one_cell requires a CopyData object
3122  * even if it is no used. The next line silences a warning about this
3123  * unused variable.
3124  *
3125  * @code
3126  * (void)copy_data;
3127  *
3128  * @endcode
3129  *
3130  * First task on each cell is to compute the cell residual
3131  * contributions of this cell, and put them into the
3132  * <code>error_indicators</code> variable:
3133  *
3134  * @code
3135  * integrate_over_cell(cell,
3136  * scratch_data.primal_solution,
3137  * scratch_data.dual_weights,
3138  * scratch_data.cell_data,
3139  * error_indicators);
3140  *
3141  * @endcode
3142  *
3143  * After computing the cell terms, turn to the face terms. For this,
3144  * loop over all faces of the present cell, and see whether
3145  * something needs to be computed on it:
3146  *
3147  * @code
3148  * for (unsigned int face_no : GeometryInfo<dim>::face_indices())
3149  * {
3150  * @endcode
3151  *
3152  * First, if this face is part of the boundary, then there is
3153  * nothing to do. However, to make things easier when summing up
3154  * the contributions of the faces of cells, we enter this face
3155  * into the list of faces with a zero contribution to the error.
3156  *
3157  * @code
3158  * if (cell->face(face_no)->at_boundary())
3159  * {
3160  * face_integrals[cell->face(face_no)] = 0;
3161  * continue;
3162  * }
3163  *
3164  * @endcode
3165  *
3166  * Next, note that since we want to compute the jump terms on
3167  * each face only once although we access it twice (if it is not
3168  * at the boundary), we have to define some rules who is
3169  * responsible for computing on a face:
3170  *
3171 
3172  *
3173  * First, if the neighboring cell is on the same level as this
3174  * one, i.e. neither further refined not coarser, then the one
3175  * with the lower index within this level does the work. In
3176  * other words: if the other one has a lower index, then skip
3177  * work on this face:
3178  *
3179  * @code
3180  * if ((cell->neighbor(face_no)->has_children() == false) &&
3181  * (cell->neighbor(face_no)->level() == cell->level()) &&
3182  * (cell->neighbor(face_no)->index() < cell->index()))
3183  * continue;
3184  *
3185  * @endcode
3186  *
3187  * Likewise, we always work from the coarser cell if this and
3188  * its neighbor differ in refinement. Thus, if the neighboring
3189  * cell is less refined than the present one, then do nothing
3190  * since we integrate over the subfaces when we visit the coarse
3191  * cell.
3192  *
3193  * @code
3194  * if (cell->at_boundary(face_no) == false)
3195  * if (cell->neighbor(face_no)->level() < cell->level())
3196  * continue;
3197  *
3198  *
3199  * @endcode
3200  *
3201  * Now we know that we are in charge here, so actually compute
3202  * the face jump terms. If the face is a regular one, i.e. the
3203  * other side's cell is neither coarser not finer than this
3204  * cell, then call one function, and if the cell on the other
3205  * side is further refined, then use another function. Note that
3206  * the case that the cell on the other side is coarser cannot
3207  * happen since we have decided above that we handle this case
3208  * when we pass over that other cell.
3209  *
3210  * @code
3211  * if (cell->face(face_no)->has_children() == false)
3212  * integrate_over_regular_face(cell,
3213  * face_no,
3214  * scratch_data.primal_solution,
3215  * scratch_data.dual_weights,
3216  * scratch_data.face_data,
3217  * face_integrals);
3218  * else
3219  * integrate_over_irregular_face(cell,
3220  * face_no,
3221  * scratch_data.primal_solution,
3222  * scratch_data.dual_weights,
3223  * scratch_data.face_data,
3224  * face_integrals);
3225  * }
3226  * }
3227  *
3228  *
3229  * @endcode
3230  *
3231  *
3232  * <a name="Computingcelltermerrorcontributions"></a>
3233  * <h4>Computing cell term error contributions</h4>
3234  *
3235 
3236  *
3237  * As for the actual computation of the error contributions, first turn to
3238  * the cell terms:
3239  *
3240  * @code
3241  * template <int dim>
3242  * void WeightedResidual<dim>::integrate_over_cell(
3243  * const active_cell_iterator &cell,
3244  * const Vector<double> & primal_solution,
3245  * const Vector<double> & dual_weights,
3246  * CellData & cell_data,
3247  * Vector<float> & error_indicators) const
3248  * {
3249  * @endcode
3250  *
3251  * The tasks to be done are what appears natural from looking at the
3252  * error estimation formula: first get the right hand side and Laplacian
3253  * of the numerical solution at the quadrature points for the cell
3254  * residual,
3255  *
3256  * @code
3257  * cell_data.fe_values.reinit(cell);
3258  * cell_data.right_hand_side->value_list(
3259  * cell_data.fe_values.get_quadrature_points(), cell_data.rhs_values);
3260  * cell_data.fe_values.get_function_laplacians(primal_solution,
3261  * cell_data.cell_laplacians);
3262  *
3263  * @endcode
3264  *
3265  * ...then get the dual weights...
3266  *
3267  * @code
3268  * cell_data.fe_values.get_function_values(dual_weights,
3269  * cell_data.dual_weights);
3270  *
3271  * @endcode
3272  *
3273  * ...and finally build the sum over all quadrature points and store it
3274  * with the present cell:
3275  *
3276  * @code
3277  * double sum = 0;
3278  * for (unsigned int p = 0; p < cell_data.fe_values.n_quadrature_points; ++p)
3279  * sum += ((cell_data.rhs_values[p] + cell_data.cell_laplacians[p]) *
3280  * cell_data.dual_weights[p] * cell_data.fe_values.JxW(p));
3281  * error_indicators(cell->active_cell_index()) += sum;
3282  * }
3283  *
3284  *
3285  * @endcode
3286  *
3287  *
3288  * <a name="Computingedgetermerrorcontributions1"></a>
3289  * <h4>Computing edge term error contributions -- 1</h4>
3290  *
3291 
3292  *
3293  * On the other hand, computation of the edge terms for the error estimate
3294  * is not so simple. First, we have to distinguish between faces with and
3295  * without hanging nodes. Because it is the simple case, we first consider
3296  * the case without hanging nodes on a face (let's call this the `regular'
3297  * case):
3298  *
3299  * @code
3300  * template <int dim>
3301  * void WeightedResidual<dim>::integrate_over_regular_face(
3302  * const active_cell_iterator &cell,
3303  * const unsigned int face_no,
3304  * const Vector<double> & primal_solution,
3305  * const Vector<double> & dual_weights,
3306  * FaceData & face_data,
3307  * FaceIntegrals & face_integrals) const
3308  * {
3309  * const unsigned int n_q_points =
3310  * face_data.fe_face_values_cell.n_quadrature_points;
3311  *
3312  * @endcode
3313  *
3314  * The first step is to get the values of the gradients at the
3315  * quadrature points of the finite element field on the present
3316  * cell. For this, initialize the <code>FEFaceValues</code> object
3317  * corresponding to this side of the face, and extract the gradients
3318  * using that object.
3319  *
3320  * @code
3321  * face_data.fe_face_values_cell.reinit(cell, face_no);
3322  * face_data.fe_face_values_cell.get_function_gradients(
3323  * primal_solution, face_data.cell_grads);
3324  *
3325  * @endcode
3326  *
3327  * The second step is then to extract the gradients of the finite
3328  * element solution at the quadrature points on the other side of the
3329  * face, i.e. from the neighboring cell.
3330  *
3331 
3332  *
3333  * For this, do a sanity check before: make sure that the neighbor
3334  * actually exists (yes, we should not have come here if the neighbor
3335  * did not exist, but in complicated software there are bugs, so better
3336  * check this), and if this is not the case throw an error.
3337  *
3338  * @code
3339  * Assert(cell->neighbor(face_no).state() == IteratorState::valid,
3340  * ExcInternalError());
3341  * @endcode
3342  *
3343  * If we have that, then we need to find out with which face of the
3344  * neighboring cell we have to work, i.e. the <code>how-many'th</code> the
3345  * neighbor the present cell is of the cell behind the present face. For
3346  * this, there is a function, and we put the result into a variable with
3347  * the name <code>neighbor_neighbor</code>:
3348  *
3349  * @code
3350  * const unsigned int neighbor_neighbor =
3351  * cell->neighbor_of_neighbor(face_no);
3352  * @endcode
3353  *
3354  * Then define an abbreviation for the neighbor cell, initialize the
3355  * <code>FEFaceValues</code> object on that cell, and extract the
3356  * gradients on that cell:
3357  *
3358  * @code
3359  * const active_cell_iterator neighbor = cell->neighbor(face_no);
3360  * face_data.fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
3361  * face_data.fe_face_values_neighbor.get_function_gradients(
3362  * primal_solution, face_data.neighbor_grads);
3363  *
3364  * @endcode
3365  *
3366  * Now that we have the gradients on this and the neighboring cell,
3367  * compute the jump residual by multiplying the jump in the gradient
3368  * with the normal vector:
3369  *
3370  * @code
3371  * for (unsigned int p = 0; p < n_q_points; ++p)
3372  * face_data.jump_residual[p] =
3373  * ((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
3374  * face_data.fe_face_values_cell.normal_vector(p));
3375  *
3376  * @endcode
3377  *
3378  * Next get the dual weights for this face:
3379  *
3380  * @code
3381  * face_data.fe_face_values_cell.get_function_values(dual_weights,
3382  * face_data.dual_weights);
3383  *
3384  * @endcode
3385  *
3386  * Finally, we have to compute the sum over jump residuals, dual
3387  * weights, and quadrature weights, to get the result for this face:
3388  *
3389  * @code
3390  * double face_integral = 0;
3391  * for (unsigned int p = 0; p < n_q_points; ++p)
3392  * face_integral +=
3393  * (face_data.jump_residual[p] * face_data.dual_weights[p] *
3394  * face_data.fe_face_values_cell.JxW(p));
3395  *
3396  * @endcode
3397  *
3398  * Double check that the element already exists and that it was not
3399  * already written to...
3400  *
3401  * @code
3402  * Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
3403  * ExcInternalError());
3404  * Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
3405  *
3406  * @endcode
3407  *
3408  * ...then store computed value at assigned location. Note that the
3409  * stored value does not contain the factor 1/2 that appears in the
3410  * error representation. The reason is that the term actually does not
3411  * have this factor if we loop over all faces in the triangulation, but
3412  * only appears if we write it as a sum over all cells and all faces of
3413  * each cell; we thus visit the same face twice. We take account of this
3414  * by using this factor -1/2 later, when we sum up the contributions for
3415  * each cell individually.
3416  *
3417  * @code
3418  * face_integrals[cell->face(face_no)] = face_integral;
3419  * }
3420  *
3421  *
3422  * @endcode
3423  *
3424  *
3425  * <a name="Computingedgetermerrorcontributions2"></a>
3426  * <h4>Computing edge term error contributions -- 2</h4>
3427  *
3428 
3429  *
3430  * We are still missing the case of faces with hanging nodes. This is what
3431  * is covered in this function:
3432  *
3433  * @code
3434  * template <int dim>
3435  * void WeightedResidual<dim>::integrate_over_irregular_face(
3436  * const active_cell_iterator &cell,
3437  * const unsigned int face_no,
3438  * const Vector<double> & primal_solution,
3439  * const Vector<double> & dual_weights,
3440  * FaceData & face_data,
3441  * FaceIntegrals & face_integrals) const
3442  * {
3443  * @endcode
3444  *
3445  * First again two abbreviations, and some consistency checks whether
3446  * the function is called only on faces for which it is supposed to be
3447  * called:
3448  *
3449  * @code
3450  * const unsigned int n_q_points =
3451  * face_data.fe_face_values_cell.n_quadrature_points;
3452  *
3453  * const typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
3454  * const typename DoFHandler<dim>::cell_iterator neighbor =
3455  * cell->neighbor(face_no);
3456  * Assert(neighbor.state() == IteratorState::valid, ExcInternalError());
3457  * Assert(neighbor->has_children(), ExcInternalError());
3458  * (void)neighbor;
3459  *
3460  * @endcode
3461  *
3462  * Then find out which neighbor the present cell is of the adjacent
3463  * cell. Note that we will operate on the children of this adjacent
3464  * cell, but that their orientation is the same as that of their mother,
3465  * i.e. the neighbor direction is the same.
3466  *
3467  * @code
3468  * const unsigned int neighbor_neighbor =
3469  * cell->neighbor_of_neighbor(face_no);
3470  *
3471  * @endcode
3472  *
3473  * Then simply do everything we did in the previous function for one
3474  * face for all the sub-faces now:
3475  *
3476  * @code
3477  * for (unsigned int subface_no = 0; subface_no < face->n_children();
3478  * ++subface_no)
3479  * {
3480  * @endcode
3481  *
3482  * Start with some checks again: get an iterator pointing to the
3483  * cell behind the present subface and check whether its face is a
3484  * subface of the one we are considering. If that were not the case,
3485  * then there would be either a bug in the
3486  * <code>neighbor_neighbor</code> function called above, or -- worse
3487  * -- some function in the library did not keep to some underlying
3488  * assumptions about cells, their children, and their faces. In any
3489  * case, even though this assertion should not be triggered, it does
3490  * not harm to be cautious, and in optimized mode computations the
3491  * assertion will be removed anyway.
3492  *
3493  * @code
3494  * const active_cell_iterator neighbor_child =
3495  * cell->neighbor_child_on_subface(face_no, subface_no);
3496  * Assert(neighbor_child->face(neighbor_neighbor) ==
3497  * cell->face(face_no)->child(subface_no),
3498  * ExcInternalError());
3499  *
3500  * @endcode
3501  *
3502  * Now start the work by again getting the gradient of the solution
3503  * first at this side of the interface,
3504  *
3505  * @code
3506  * face_data.fe_subface_values_cell.reinit(cell, face_no, subface_no);
3507  * face_data.fe_subface_values_cell.get_function_gradients(
3508  * primal_solution, face_data.cell_grads);
3509  * @endcode
3510  *
3511  * then at the other side,
3512  *
3513  * @code
3514  * face_data.fe_face_values_neighbor.reinit(neighbor_child,
3515  * neighbor_neighbor);
3516  * face_data.fe_face_values_neighbor.get_function_gradients(
3517  * primal_solution, face_data.neighbor_grads);
3518  *
3519  * @endcode
3520  *
3521  * and finally building the jump residuals. Since we take the normal
3522  * vector from the other cell this time, revert the sign of the
3523  * first term compared to the other function:
3524  *
3525  * @code
3526  * for (unsigned int p = 0; p < n_q_points; ++p)
3527  * face_data.jump_residual[p] =
3528  * ((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
3529  * face_data.fe_face_values_neighbor.normal_vector(p));
3530  *
3531  * @endcode
3532  *
3533  * Then get dual weights:
3534  *
3535  * @code
3536  * face_data.fe_face_values_neighbor.get_function_values(
3537  * dual_weights, face_data.dual_weights);
3538  *
3539  * @endcode
3540  *
3541  * At last, sum up the contribution of this sub-face, and set it in
3542  * the global map:
3543  *
3544  * @code
3545  * double face_integral = 0;
3546  * for (unsigned int p = 0; p < n_q_points; ++p)
3547  * face_integral +=
3548  * (face_data.jump_residual[p] * face_data.dual_weights[p] *
3549  * face_data.fe_face_values_neighbor.JxW(p));
3550  * face_integrals[neighbor_child->face(neighbor_neighbor)] =
3551  * face_integral;
3552  * }
3553  *
3554  * @endcode
3555  *
3556  * Once the contributions of all sub-faces are computed, loop over all
3557  * sub-faces to collect and store them with the mother face for simple
3558  * use when later collecting the error terms of cells. Again make safety
3559  * checks that the entries for the sub-faces have been computed and do
3560  * not carry an invalid value.
3561  *
3562  * @code
3563  * double sum = 0;
3564  * for (unsigned int subface_no = 0; subface_no < face->n_children();
3565  * ++subface_no)
3566  * {
3567  * Assert(face_integrals.find(face->child(subface_no)) !=
3568  * face_integrals.end(),
3569  * ExcInternalError());
3570  * Assert(face_integrals[face->child(subface_no)] != -1e20,
3571  * ExcInternalError());
3572  *
3573  * sum += face_integrals[face->child(subface_no)];
3574  * }
3575  * @endcode
3576  *
3577  * Finally store the value with the parent face.
3578  *
3579  * @code
3580  * face_integrals[face] = sum;
3581  * }
3582  *
3583  * } // namespace LaplaceSolver
3584  *
3585  *
3586  * @endcode
3587  *
3588  *
3589  * <a name="Asimulationframework"></a>
3590  * <h3>A simulation framework</h3>
3591  *
3592 
3593  *
3594  * In the previous example program, we have had two functions that were used
3595  * to drive the process of solving on subsequently finer grids. We extend
3596  * this here to allow for a number of parameters to be passed to these
3597  * functions, and put all of that into framework class.
3598  *
3599 
3600  *
3601  * You will have noted that this program is built up of a number of small
3602  * parts (evaluation functions, solver classes implementing various
3603  * refinement methods, different dual functionals, different problem and
3604  * data descriptions), which makes the program relatively simple to extend,
3605  * but also allows to solve a large number of different problems by
3606  * replacing one part by another. We reflect this flexibility by declaring a
3607  * structure in the following framework class that holds a number of
3608  * parameters that may be set to test various combinations of the parts of
3609  * this program, and which can be used to test it at various problems and
3610  * discretizations in a simple way.
3611  *
3612  * @code
3613  * template <int dim>
3614  * struct Framework
3615  * {
3616  * public:
3617  * @endcode
3618  *
3619  * First, we declare two abbreviations for simple use of the respective
3620  * data types:
3621  *
3622  * @code
3623  * using Evaluator = Evaluation::EvaluationBase<dim>;
3624  * using EvaluatorList = std::list<Evaluator *>;
3625  *
3626  *
3627  * @endcode
3628  *
3629  * Then we have the structure which declares all the parameters that may
3630  * be set. In the default constructor of the structure, these values are
3631  * all set to default values, for simple use.
3632  *
3633  * @code
3634  * struct ProblemDescription
3635  * {
3636  * @endcode
3637  *
3638  * First allow for the degrees of the piecewise polynomials by which the
3639  * primal and dual problems will be discretized. They default to (bi-,
3640  * tri-)linear ansatz functions for the primal, and (bi-, tri-)quadratic
3641  * ones for the dual problem. If a refinement criterion is chosen that
3642  * does not need the solution of a dual problem, the value of the dual
3643  * finite element degree is of course ignored.
3644  *
3645  * @code
3646  * unsigned int primal_fe_degree;
3647  * unsigned int dual_fe_degree;
3648  *
3649  * @endcode
3650  *
3651  * Then have an object that describes the problem type, i.e. right hand
3652  * side, domain, boundary values, etc. The pointer needed here defaults
3653  * to the Null pointer, i.e. you will have to set it in actual instances
3654  * of this object to make it useful.
3655  *
3656  * @code
3657  * std::unique_ptr<const Data::SetUpBase<dim>> data;
3658  *
3659  * @endcode
3660  *
3661  * Since we allow to use different refinement criteria (global
3662  * refinement, refinement by the Kelly error indicator, possibly with a
3663  * weight, and using the dual estimator), define a number of enumeration
3664  * values, and subsequently a variable of that type. It will default to
3665  * <code>dual_weighted_error_estimator</code>.
3666  *
3667  * @code
3668  * enum RefinementCriterion
3669  * {
3670  * dual_weighted_error_estimator,
3671  * global_refinement,
3672  * kelly_indicator,
3673  * weighted_kelly_indicator
3674  * };
3675  *
3676  * RefinementCriterion refinement_criterion;
3677  *
3678  * @endcode
3679  *
3680  * Next, an object that describes the dual functional. It is only needed
3681  * if the dual weighted residual refinement is chosen, and also defaults
3682  * to a Null pointer.
3683  *
3684  * @code
3685  * std::unique_ptr<const DualFunctional::DualFunctionalBase<dim>>
3686  * dual_functional;
3687  *
3688  * @endcode
3689  *
3690  * Then a list of evaluation objects. Its default value is empty,
3691  * i.e. no evaluation objects.
3692  *
3693  * @code
3694  * EvaluatorList evaluator_list;
3695  *
3696  * @endcode
3697  *
3698  * Next to last, a function that is used as a weight to the
3699  * <code>RefinementWeightedKelly</code> class. The default value of this
3700  * pointer is zero, but you have to set it to some other value if you
3701  * want to use the <code>weighted_kelly_indicator</code> refinement
3702  * criterion.
3703  *
3704  * @code
3705  * std::unique_ptr<const Function<dim>> kelly_weight;
3706  *
3707  * @endcode
3708  *
3709  * Finally, we have a variable that denotes the maximum number of
3710  * degrees of freedom we allow for the (primal) discretization. If it is
3711  * exceeded, we stop the process of solving and intermittent mesh
3712  * refinement. Its default value is 20,000.
3713  *
3714  * @code
3715  * unsigned int max_degrees_of_freedom;
3716  *
3717  * @endcode
3718  *
3719  * Finally the default constructor of this class:
3720  *
3721  * @code
3722  * ProblemDescription();
3723  * };
3724  *
3725  * @endcode
3726  *
3727  * The driver framework class only has one method which calls solver and
3728  * mesh refinement intermittently, and does some other small tasks in
3729  * between. Since it does not need data besides the parameters given to
3730  * it, we make it static:
3731  *
3732  * @code
3733  * static void run(const ProblemDescription &descriptor);
3734  * };
3735  *
3736  *
3737  * @endcode
3738  *
3739  * As for the implementation, first the constructor of the parameter object,
3740  * setting all values to their defaults:
3741  *
3742  * @code
3743  * template <int dim>
3744  * Framework<dim>::ProblemDescription::ProblemDescription()
3745  * : primal_fe_degree(1)
3746  * , dual_fe_degree(2)
3747  * , refinement_criterion(dual_weighted_error_estimator)
3748  * , max_degrees_of_freedom(20000)
3749  * {}
3750  *
3751  *
3752  *
3753  * @endcode
3754  *
3755  * Then the function which drives the whole process:
3756  *
3757  * @code
3758  * template <int dim>
3759  * void Framework<dim>::run(const ProblemDescription &descriptor)
3760  * {
3761  * @endcode
3762  *
3763  * First create a triangulation from the given data object,
3764  *
3765  * @code
3768  * descriptor.data->create_coarse_grid(triangulation);
3769  *
3770  * @endcode
3771  *
3772  * then a set of finite elements and appropriate quadrature formula:
3773  *
3774  * @code
3775  * const FE_Q<dim> primal_fe(descriptor.primal_fe_degree);
3776  * const FE_Q<dim> dual_fe(descriptor.dual_fe_degree);
3777  * const QGauss<dim> quadrature(descriptor.dual_fe_degree + 1);
3778  * const QGauss<dim - 1> face_quadrature(descriptor.dual_fe_degree + 1);
3779  *
3780  * @endcode
3781  *
3782  * Next, select one of the classes implementing different refinement
3783  * criteria.
3784  *
3785  * @code
3786  * std::unique_ptr<LaplaceSolver::Base<dim>> solver;
3787  * switch (descriptor.refinement_criterion)
3788  * {
3789  * case ProblemDescription::dual_weighted_error_estimator:
3790  * {
3791  * solver =
3792  * std_cxx14::make_unique<LaplaceSolver::WeightedResidual<dim>>(
3793  * triangulation,
3794  * primal_fe,
3795  * dual_fe,
3796  * quadrature,
3797  * face_quadrature,
3798  * descriptor.data->get_right_hand_side(),
3799  * descriptor.data->get_boundary_values(),
3800  * *descriptor.dual_functional);
3801  * break;
3802  * }
3803  *
3804  * case ProblemDescription::global_refinement:
3805  * {
3806  * solver =
3807  * std_cxx14::make_unique<LaplaceSolver::RefinementGlobal<dim>>(
3808  * triangulation,
3809  * primal_fe,
3810  * quadrature,
3811  * face_quadrature,
3812  * descriptor.data->get_right_hand_side(),
3813  * descriptor.data->get_boundary_values());
3814  * break;
3815  * }
3816  *
3817  * case ProblemDescription::kelly_indicator:
3818  * {
3819  * solver =
3820  * std_cxx14::make_unique<LaplaceSolver::RefinementKelly<dim>>(
3821  * triangulation,
3822  * primal_fe,
3823  * quadrature,
3824  * face_quadrature,
3825  * descriptor.data->get_right_hand_side(),
3826  * descriptor.data->get_boundary_values());
3827  * break;
3828  * }
3829  *
3830  * case ProblemDescription::weighted_kelly_indicator:
3831  * {
3832  * solver = std_cxx14::make_unique<
3833  * LaplaceSolver::RefinementWeightedKelly<dim>>(
3834  * triangulation,
3835  * primal_fe,
3836  * quadrature,
3837  * face_quadrature,
3838  * descriptor.data->get_right_hand_side(),
3839  * descriptor.data->get_boundary_values(),
3840  * *descriptor.kelly_weight);
3841  * break;
3842  * }
3843  *
3844  * default:
3845  * AssertThrow(false, ExcInternalError());
3846  * }
3847  *
3848  * @endcode
3849  *
3850  * Now that all objects are in place, run the main loop. The stopping
3851  * criterion is implemented at the bottom of the loop.
3852  *
3853 
3854  *
3855  * In the loop, first set the new cycle number, then solve the problem,
3856  * output its solution(s), apply the evaluation objects to it, then decide
3857  * whether we want to refine the mesh further and solve again on this
3858  * mesh, or jump out of the loop.
3859  *
3860  * @code
3861  * for (unsigned int step = 0; true; ++step)
3862  * {
3863  * std::cout << "Refinement cycle: " << step << std::endl;
3864  *
3865  * solver->set_refinement_cycle(step);
3866  * solver->solve_problem();
3867  * solver->output_solution();
3868  *
3869  * std::cout << " Number of degrees of freedom=" << solver->n_dofs()
3870  * << std::endl;
3871  *
3872  * for (const auto &evaluator : descriptor.evaluator_list)
3873  * {
3874  * evaluator->set_refinement_cycle(step);
3875  * solver->postprocess(*evaluator);
3876  * }
3877  *
3878  *
3879  * if (solver->n_dofs() < descriptor.max_degrees_of_freedom)
3880  * solver->refine_grid();
3881  * else
3882  * break;
3883  * }
3884  *
3885  * @endcode
3886  *
3887  * Clean up the screen after the loop has run:
3888  *
3889  * @code
3890  * std::cout << std::endl;
3891  * }
3892  *
3893  * } // namespace Step14
3894  *
3895  *
3896  *
3897  * @endcode
3898  *
3899  *
3900  * <a name="Themainfunction"></a>
3901  * <h3>The main function</h3>
3902  *
3903 
3904  *
3905  * Here finally comes the main function. It drives the whole process by
3906  * specifying a set of parameters to be used for the simulation (polynomial
3907  * degrees, evaluation and dual functionals, etc), and passes them packed into
3908  * a structure to the frame work class above.
3909  *
3910  * @code
3911  * int main()
3912  * {
3913  * try
3914  * {
3915  * using namespace Step14;
3916  *
3917  * @endcode
3918  *
3919  * Describe the problem we want to solve here by passing a descriptor
3920  * object to the function doing the rest of the work:
3921  *
3922  * @code
3923  * const unsigned int dim = 2;
3924  * Framework<dim>::ProblemDescription descriptor;
3925  *
3926  * @endcode
3927  *
3928  * First set the refinement criterion we wish to use:
3929  *
3930  * @code
3931  * descriptor.refinement_criterion =
3932  * Framework<dim>::ProblemDescription::dual_weighted_error_estimator;
3933  * @endcode
3934  *
3935  * Here, we could as well have used <code>global_refinement</code> or
3936  * <code>weighted_kelly_indicator</code>. Note that the information
3937  * given about dual finite elements, dual functional, etc is only
3938  * important for the given choice of refinement criterion, and is
3939  * ignored otherwise.
3940  *
3941 
3942  *
3943  * Then set the polynomial degrees of primal and dual problem. We choose
3944  * here bi-linear and bi-quadratic ones:
3945  *
3946  * @code
3947  * descriptor.primal_fe_degree = 1;
3948  * descriptor.dual_fe_degree = 2;
3949  *
3950  * @endcode
3951  *
3952  * Then set the description of the test case, i.e. domain, boundary
3953  * values, and right hand side. These are prepackaged in classes. We
3954  * take here the description of <code>Exercise_2_3</code>, but you can
3955  * also use <code>CurvedRidges@<dim@></code>:
3956  *
3957  * @code
3958  * descriptor.data =
3959  * std_cxx14::make_unique<Data::SetUp<Data::Exercise_2_3<dim>, dim>>();
3960  *
3961  * @endcode
3962  *
3963  * Next set first a dual functional, then a list of evaluation
3964  * objects. We choose as default the evaluation of the value at an
3965  * evaluation point, represented by the classes
3966  * <code>PointValueEvaluation</code> in the namespaces of evaluation and
3967  * dual functional classes. You can also set the
3968  * <code>PointXDerivativeEvaluation</code> classes for the x-derivative
3969  * instead of the value at the evaluation point.
3970  *
3971 
3972  *
3973  * Note that dual functional and evaluation objects should
3974  * match. However, you can give as many evaluation functionals as you
3975  * want, so you can have both point value and derivative evaluated after
3976  * each step. One such additional evaluation is to output the grid in
3977  * each step.
3978  *
3979  * @code
3980  * const Point<dim> evaluation_point(0.75, 0.75);
3981  * descriptor.dual_functional =
3982  * std_cxx14::make_unique<DualFunctional::PointValueEvaluation<dim>>(
3983  * evaluation_point);
3984  *
3985  * Evaluation::PointValueEvaluation<dim> postprocessor1(evaluation_point);
3986  * Evaluation::GridOutput<dim> postprocessor2("grid");
3987  *
3988  * descriptor.evaluator_list.push_back(&postprocessor1);
3989  * descriptor.evaluator_list.push_back(&postprocessor2);
3990  *
3991  * @endcode
3992  *
3993  * Set the maximal number of degrees of freedom after which we want the
3994  * program to stop refining the mesh further:
3995  *
3996  * @code
3997  * descriptor.max_degrees_of_freedom = 20000;
3998  *
3999  * @endcode
4000  *
4001  * Finally pass the descriptor object to a function that runs the entire
4002  * solution with it:
4003  *
4004  * @code
4005  * Framework<dim>::run(descriptor);
4006  * }
4007  *
4008  * @endcode
4009  *
4010  * Catch exceptions to give information about things that failed:
4011  *
4012  * @code
4013  * catch (std::exception &exc)
4014  * {
4015  * std::cerr << std::endl
4016  * << std::endl
4017  * << "----------------------------------------------------"
4018  * << std::endl;
4019  * std::cerr << "Exception on processing: " << std::endl
4020  * << exc.what() << std::endl
4021  * << "Aborting!" << std::endl
4022  * << "----------------------------------------------------"
4023  * << std::endl;
4024  * return 1;
4025  * }
4026  * catch (...)
4027  * {
4028  * std::cerr << std::endl
4029  * << std::endl
4030  * << "----------------------------------------------------"
4031  * << std::endl;
4032  * std::cerr << "Unknown exception!" << std::endl
4033  * << "Aborting!" << std::endl
4034  * << "----------------------------------------------------"
4035  * << std::endl;
4036  * return 1;
4037  * }
4038  *
4039  * return 0;
4040  * }
4041  * @endcode
4042 <a name="Results"></a><h1>Results</h1>
4043 
4044 
4045 <a name="Pointvalues"></a><h3>Point values</h3>
4046 
4047 
4048 
4049 This program offers a lot of possibilities to play around. We can thus
4050 only show a small part of all possible results that can be obtained
4051 with the help of this program. However, you are encouraged to just try
4052 it out, by changing the settings in the main program. Here, we start
4053 by simply letting it run, unmodified:
4054 @code
4055 Refinement cycle: 0
4056  Number of degrees of freedom=72
4057  Point value=0.03243
4058  Estimated error=0.000702385
4059 Refinement cycle: 1
4060  Number of degrees of freedom=67
4061  Point value=0.0324827
4062  Estimated error=0.000888953
4063 Refinement cycle: 2
4064  Number of degrees of freedom=130
4065  Point value=0.0329619
4066  Estimated error=0.000454606
4067 Refinement cycle: 3
4068  Number of degrees of freedom=307
4069  Point value=0.0331934
4070  Estimated error=0.000241254
4071 Refinement cycle: 4
4072  Number of degrees of freedom=718
4073  Point value=0.0333675
4074  Estimated error=7.4912e-05
4075 Refinement cycle: 5
4076  Number of degrees of freedom=1665
4077  Point value=0.0334083
4078  Estimated error=3.69111e-05
4079 Refinement cycle: 6
4080  Number of degrees of freedom=3975
4081  Point value=0.033431
4082  Estimated error=1.54218e-05
4083 Refinement cycle: 7
4084  Number of degrees of freedom=8934
4085  Point value=0.0334406
4086  Estimated error=6.28359e-06
4087 Refinement cycle: 8
4088  Number of degrees of freedom=21799
4089  Point value=0.0334444
4090 @endcode
4091 
4092 
4093 First let's look what the program actually computed. On the seventh
4094 grid, primal and dual numerical solutions look like this (using a
4095 color scheme intended to evoke the snow-capped mountains of
4096 Colorado that the original author of this program now calls
4097 home):
4098 <table align="center">
4099  <tr>
4100  <td width="50%">
4101  <img src="https://www.dealii.org/images/steps/developer/step-14.point-value.solution-7.9.2.png" alt="">
4102  </td>
4103  <td width="50%">
4104  <img src="https://www.dealii.org/images/steps/developer/step-14.point-value.solution-7-dual.9.2.png" alt="">
4105  </td>
4106  </tr>
4107 </table>
4108 Apparently, the region at the bottom left is so unimportant for the
4109 point value evaluation at the top right that the grid is left entirely
4110 unrefined there, even though the solution has singularities at the inner
4111 corner of that cell! Due
4112 to the symmetry in right hand side and domain, the solution should
4113 actually look like at the top right in all four corners, but the mesh
4114 refinement criterion involving the dual solution chose to refine them
4115 differently -- because we said that we really only care about a single
4116 function value somewhere at the top right.
4117 
4118 
4119 
4120 Here are some of the meshes that are produced in refinement cycles 0,
4121 2, 4 (top row), and 5, 7, and 8 (bottom row):
4122 
4123 <table width="80%" align="center">
4124  <tr>
4125  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-0.9.2.png" alt="" width="100%"></td>
4126  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-2.9.2.png" alt="" width="100%"></td>
4127  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-4.9.2.png" alt="" width="100%"></td>
4128  </tr>
4129  <tr>
4130  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-5.9.2.png" alt="" width="100%"></td>
4131  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-7.9.2.png" alt="" width="100%"></td>
4132  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-8.9.2.png" alt="" width="100%"></td>
4133  </tr>
4134 </table>
4135 
4136 Note the subtle interplay between resolving the corner singularities,
4137 and resolving around the point of evaluation. It will be rather
4138 difficult to generate such a mesh by hand, as this would involve to
4139 judge quantitatively how much which of the four corner singularities
4140 should be resolved, and to set the weight compared to the vicinity of
4141 the evaluation point.
4142 
4143 
4144 
4145 The program prints the point value and the estimated error in this
4146 quantity. From extrapolating it, we can guess that the exact value is
4147 somewhere close to 0.0334473, plus or minus 0.0000001 (note that we get
4148 almost 6 valid digits from only 22,000 (primal) degrees of
4149 freedom. This number cannot be obtained from the value of the
4150 functional alone, but I have used the assumption that the error
4151 estimator is mostly exact, and extrapolated the computed value plus
4152 the estimated error, to get an approximation of the true
4153 value. Computing with more degrees of freedom shows that this
4154 assumption is indeed valid.
4155 
4156 
4157 
4158 From the computed results, we can generate two graphs: one that shows
4159 the convergence of the error @f$J(u)-J(u_h)@f$ (taking the
4160 extrapolated value as correct) in the point value, and the value that
4161 we get by adding up computed value @f$J(u_h)@f$ and estimated
4162 error eta (if the error estimator @f$eta@f$ were exact, then the value
4163 @f$J(u_h)+\eta@f$ would equal the exact point value, and the error
4164 in this quantity would always be zero; however, since the error
4165 estimator is only a - good - approximation to the true error, we can
4166 by this only reduce the size of the error). In this graph, we also
4167 indicate the complexity @f${\cal O}(1/N)@f$ to show that mesh refinement
4168 acts optimal in this case. The second chart compares
4169 true and estimated error, and shows that the two are actually very
4170 close to each other, even for such a complicated quantity as the point
4171 value:
4172 
4173 
4174 <table width="80%" align="center">
4175  <tr>
4176  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.error.png" alt="" width="100%"></td>
4177  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.error-estimation.png" alt="" width="100%"></td>
4178  </tr>
4179 </table>
4180 
4181 
4182 <a name="Comparingrefinementcriteria"></a><h3>Comparing refinement criteria</h3>
4183 
4184 
4185 
4186 Since we have accepted quite some effort when using the mesh
4187 refinement driven by the dual weighted error estimator (for solving
4188 the dual problem, and for evaluating the error representation), it is
4189 worth while asking whether that effort was successful. To this end, we
4190 first compare the achieved error levels for different mesh refinement
4191 criteria. To generate this data, simply change the value of the mesh
4192 refinement criterion variable in the main program. The results are
4193 thus (for the weight in the Kelly indicator, we have chosen the
4194 function @f$1/(r^2+0.1^2)@f$, where @f$r@f$
4195 is the distance to the evaluation point; it can be shown that this is
4196 the optimal weight if we neglect the effects of boundaries):
4197 
4198 <img src="https://www.dealii.org/images/steps/developer/step-14.point-value.error-comparison.png" alt="">
4199 
4200 
4201 
4202 Checking these numbers, we see that for global refinement, the error
4203 is proportional to @f$O(1/(sqrt(N) log(N)))@f$, and for the dual
4204 estimator @f$O(1/N)@f$. Generally speaking, we see that the dual
4205 weighted error estimator is better than the other refinement
4206 indicators, at least when compared with those that have a similarly
4207 regular behavior. The Kelly indicator produces smaller errors, but
4208 jumps about the picture rather irregularly, with the error also
4209 changing signs sometimes. Therefore, its behavior does not allow to
4210 extrapolate the results to larger values of N. Furthermore, if we
4211 trust the error estimates of the dual weighted error estimator, the
4212 results can be improved by adding the estimated error to the computed
4213 values. In terms of reliability, the weighted estimator is thus better
4214 than the Kelly indicator, although the latter sometimes produces
4215 smaller errors.
4216 
4217 
4218 
4219 <a name="Evaluationofpointstresses"></a><h3>Evaluation of point stresses</h3>
4220 
4221 
4222 
4223 Besides evaluating the values of the solution at a certain point, the
4224 program also offers the possibility to evaluate the x-derivatives at a
4225 certain point, and also to tailor mesh refinement for this. To let the
4226 program compute these quantities, simply replace the two occurrences of
4227 <code>PointValueEvaluation</code> in the main function by
4228 <code>PointXDerivativeEvaluation</code>, and let the program run:
4229 @code
4230 Refinement cycle: 0
4231  Number of degrees of freedom=72
4232  Point x-derivative=-0.0719397
4233  Estimated error=-0.0126173
4234 Refinement cycle: 1
4235  Number of degrees of freedom=61
4236  Point x-derivative=-0.0707956
4237  Estimated error=-0.00774316
4238 Refinement cycle: 2
4239  Number of degrees of freedom=131
4240  Point x-derivative=-0.0568671
4241  Estimated error=-0.00313426
4242 Refinement cycle: 3
4243  Number of degrees of freedom=247
4244  Point x-derivative=-0.053033
4245  Estimated error=-0.00136114
4246 Refinement cycle: 4
4247  Number of degrees of freedom=532
4248  Point x-derivative=-0.0526429
4249  Estimated error=-0.000558868
4250 Refinement cycle: 5
4251  Number of degrees of freedom=1267
4252  Point x-derivative=-0.0526955
4253  Estimated error=-0.000220116
4254 Refinement cycle: 6
4255  Number of degrees of freedom=2864
4256  Point x-derivative=-0.0527495
4257  Estimated error=-9.46731e-05
4258 Refinement cycle: 7
4259  Number of degrees of freedom=6409
4260  Point x-derivative=-0.052785
4261  Estimated error=-4.21543e-05
4262 Refinement cycle: 8
4263  Number of degrees of freedom=14183
4264  Point x-derivative=-0.0528028
4265  Estimated error=-2.04241e-05
4266 Refinement cycle: 9
4267  Number of degrees of freedom=29902
4268  Point x-derivative=-0.052814
4269 @endcode
4270 
4271 
4272 
4273 The solution looks roughly the same as before (the exact solution of
4274 course <em>is</em> the same, only the grid changed a little), but the
4275 dual solution is now different. A close-up around the point of
4276 evaluation shows this:
4277 <table align="center">
4278  <tr>
4279  <td width="50%">
4280  <img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.solution-7-dual.png" alt="">
4281  </td>
4282  <td width="50%">
4283  <img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.solution-7-dual-close-up.png" alt="">
4284  </td>
4285 </table>
4286 This time, the grids in refinement cycles 0, 5, 6, 7, 8, and 9 look
4287 like this:
4288 
4289 <table align="center" width="80%">
4290  <tr>
4291  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-0.9.2.png" alt="" width="100%"></td>
4292  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-5.9.2.png" alt="" width="100%"></td>
4293  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-6.9.2.png" alt="" width="100%"></td>
4294  </tr>
4295  <tr>
4296  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-7.9.2.png" alt="" width="100%"></td>
4297  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-8.9.2.png" alt="" width="100%"></td>
4298  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-9.9.2.png" alt="" width="100%"></td>
4299  </tr>
4300 </table>
4301 
4302 Note the asymmetry of the grids compared with those we obtained for
4303 the point evaluation. This is due to the fact that the domain and the primal
4304 solution may be symmetric about the diagonal, but the @f$x@f$-derivative is
4305 not, and the latter enters the refinement criterion.
4306 
4307 
4308 
4309 Then, it is interesting to compare actually computed values of the
4310 quantity of interest (i.e. the x-derivative of the solution at one
4311 point) with a reference value of -0.0528223... plus or minus
4312 0.0000005. We get this reference value by computing on finer grid after
4313 some more mesh refinements, with approximately 130,000 cells.
4314 Recall that if the error is @f$O(1/N)@f$ in the optimal case, then
4315 taking a mesh with ten times more cells gives us one additional digit
4316 in the result.
4317 
4318 
4319 
4320 In the left part of the following chart, you again see the convergence
4321 of the error towards this extrapolated value, while on the right you
4322 see a comparison of true and estimated error:
4323 
4324 <table width="80%" align="center">
4325  <tr>
4326  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.error.png" alt="" width="100%"></td>
4327  <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.error-estimation.png" alt="" width="100%"></td>
4328  </tr>
4329 </table>
4330 
4331 After an initial phase where the true error changes its sign, the
4332 estimated error matches it quite well, again. Also note the dramatic
4333 improvement in the error when using the estimated error to correct the
4334 computed value of @f$J(u_h)@f$.
4335 
4336 
4337 
4338 <a name="step13revisited"></a><h3>step-13 revisited</h3>
4339 
4340 
4341 
4342 If instead of the <code>Exercise_2_3</code> data set, we choose
4343 <code>CurvedRidges</code> in the main function, and choose @f$(0.5,0.5)@f$
4344 as the evaluation point, then we can redo the
4345 computations of the previous example program, to compare whether the
4346 results obtained with the help of the dual weighted error estimator
4347 are better than those we had previously.
4348 
4349 
4350 
4351 First, the meshes after 9 adaptive refinement cycles obtained with
4352 the point evaluation and derivative evaluation refinement
4353 criteria, respectively, look like this:
4354 
4355 <table width="80%" align="center">
4356  <tr>
4357  <td><img src="https://www.dealii.org/images/steps/developer/step-14.step-13.point-value.png" alt="" width="100%"></td>
4358  <td><img src="https://www.dealii.org/images/steps/developer/step-14.step-13.point-derivative.png" alt="" width="100%"></td>
4359  </tr>
4360 </table>
4361 
4362 The features of the solution can still be seen in the mesh, but since the
4363 solution is smooth, the singularities of the dual solution entirely
4364 dominate the mesh refinement criterion, and lead to strongly
4365 concentrated meshes. The solution after the seventh refinement step looks
4366 like the following:
4367 
4368 <table width="40%" align="center">
4369  <tr>
4370  <td><img src="https://www.dealii.org/images/steps/developer/step-14.step-13.solution-7.9.2.png" alt="" width="100%"></td>
4371  </tr>
4372 </table>
4373 
4374 Obviously, the solution is worse at some places, but the mesh
4375 refinement process should have taken care that these places are not
4376 important for computing the point value.
4377 
4378 
4379 
4380 
4381 The next point is to compare the new (duality based) mesh refinement
4382 criterion with the old ones. These are the results:
4383 
4384 <img src="https://www.dealii.org/images/steps/developer/step-14.step-13.error-comparison.png" alt="">
4385 
4386 
4387 
4388 The results are, well, somewhat mixed. First, the Kelly indicator
4389 disqualifies itself by its unsteady behavior, changing the sign of the
4390 error several times, and with increasing errors under mesh
4391 refinement. The dual weighted error estimator has a monotone decrease
4392 in the error, and is better than the weighted Kelly and global
4393 refinement, but the margin is not as large as expected. This is, here,
4394 due to the fact the global refinement can exploit the regular
4395 structure of the meshes around the point of evaluation, which leads to
4396 a better order of convergence for the point error. However, if we had
4397 a mesh that is not locally rectangular, for example because we had to
4398 approximate curved boundaries, or if the coefficients were not
4399 constant, then this advantage of globally refinement meshes would
4400 vanish, while the good performance of the duality based estimator
4401 would remain.
4402 
4403 
4404 
4405 
4406 <a name="Conclusionsandoutlook"></a><h3>Conclusions and outlook</h3>
4407 
4408 
4409 
4410 The results here are not too clearly indicating the superiority of the
4411 dual weighted error estimation approach for mesh refinement over other
4412 mesh refinement criteria, such as the Kelly indicator. This is due to
4413 the relative simplicity of the shown applications. If you are not
4414 convinced yet that this approach is indeed superior, you are invited
4415 to browse through the literature indicated in the introduction, where
4416 plenty of examples are provided where the dual weighted approach can
4417 reduce the necessary numerical work by orders of magnitude, making
4418 this the only way to compute certain quantities to reasonable
4419 accuracies at all.
4420 
4421 
4422 
4423 Besides the objections you may raise against its use as a mesh
4424 refinement criterion, consider that accurate knowledge of the error in
4425 the quantity one might want to compute is of great use, since we can
4426 stop computations when we are satisfied with the accuracy. Using more
4427 traditional approaches, it is very difficult to get accurate estimates
4428 for arbitrary quantities, except for, maybe, the error in the energy
4429 norm, and we will then have no guarantee that the result we computed
4430 satisfies any requirements on its accuracy. Also, as was shown for the
4431 evaluation of point values and derivatives, the error estimate can be
4432 used to extrapolate the results, yielding much higher accuracy in the
4433 quantity we want to know.
4434 
4435 
4436 
4437 Leaving these mathematical considerations, we tried to write the
4438 program in a modular way, such that implementing another test case, or
4439 another evaluation and dual functional is simple. You are encouraged
4440 to take the program as a basis for your own experiments, and to play a
4441 little.
4442 
4443 
4444 
4445 
4446 
4447  *
4448  *
4449 <a name="PlainProg"></a>
4450 <h1> The plain program</h1>
4451 @include "step-14.cc"
4452 */
FETools::interpolation_difference
void interpolation_difference(const DoFHandler< dim, spacedim > &dof1, const InVector &z1, const FiniteElement< dim, spacedim > &fe2, OutVector &z1_difference)
DoFTools::nonzero
@ nonzero
Definition: dof_tools.h:244
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
SolverCG
Definition: solver_cg.h:98
FE_Q
Definition: fe_q.h:554
CellData
Definition: tria_description.h:67
WorkStream
Definition: work_stream.h:157
Triangulation< dim >
GridRefinement::refine_and_coarsen_fixed_number
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
Definition: grid_refinement.cc:189
Threads::Task< void >
SparseMatrix< double >
GeometryInfo
Definition: geometry_info.h:1224
DataOutBase::vtu
@ vtu
Definition: data_out_base.h:1610
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
MatrixTools::apply_boundary_values
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Definition: matrix_tools.cc:81
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:38
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >(0)), const bool project_to_boundary_first=false)
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
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition: thread_management.h:1647
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
GridRefinement::coarsen
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
Definition: grid_refinement.cc:88
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
SolverBase
Definition: solver.h:333
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition: symengine_math.cc:273
Subscriptor
Definition: subscriptor.h:62
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
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
TensorAccessors::extract
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
Definition: tensor_accessors.h:226
GridRefinement::refine_and_coarsen_fixed_fraction
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
Definition: grid_refinement.cc:257
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: dof_handler.cc:1030
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
Solver
SolverBase< VectorType > Solver
Definition: solver.h:476
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
TimeStepping::invalid
@ invalid
Definition: time_stepping.h:71
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
SymmetricTensor::sum
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
KellyErrorEstimator::estimate
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandlerType &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, typename InputVector::value_type > * > &neumann_bc, const InputVector &solution, Vector< float > &error, const ComponentMask &component_mask=ComponentMask(), const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
SparsityPattern
Definition: sparsity_pattern.h:865
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
types
Definition: types.h:31
FiniteElement< dim >
Threads::TaskGroup::join_all
void join_all() const
Definition: thread_management.h:1833
QGauss
Definition: quadrature_lib.h:40
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
SmartPointer
Definition: smartpointer.h:68
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
unsigned int
value
static const bool value
Definition: dof_tools_constraints.cc:433
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints< double >
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
Utilities::MPI::internal::Tags::enumeration
enumeration
Definition: mpi_tags.h:51
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition: fe_update_flags.h:136
GridOut::write_svg
void write_svg(const Triangulation< 2, 2 > &tria, std::ostream &out) const
Definition: grid_out.cc:1544
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
update_hessians
@ update_hessians
Second derivatives of shape functions.
Definition: fe_update_flags.h:90
Threads::TaskGroup
Definition: thread_management.h:1798
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
DataOutInterface::write
void write(std::ostream &out, const DataOutBase::OutputFormat output_format=DataOutBase::default_format) const
Definition: data_out_base.cc:7757
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
Point< dim >
Differentiation::SD::sign
Expression sign(const Expression &x)
Definition: symengine_math.cc:280
PreconditionSSOR
Definition: precondition.h:665
FullMatrix< double >
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
Function< dim >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
SolverControl
Definition: solver_control.h:67
FEFaceValues
Definition: fe.h:40
SubCellData
Definition: tria_description.h:199
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
Threads::Task::join
void join() const
Definition: thread_management.h:1525
Quadrature
Definition: quadrature.h:85
first
Point< 2 > first
Definition: grid_out.cc:4352
GridOut
Definition: grid_out.h:1020
DataOut< dim >
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
AffineConstraints::close
void close()
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
parallel
Definition: distributed.h:416
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
DataOut_DoFData::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090