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-47.h
Go to the documentation of this file.
1  = 0) const override
673  * {
674  * return std::sin(PI * p[0]) * std::sin(PI * p[1]);
675  * }
676  *
677  * virtual Tensor<1, dim>
678  * gradient(const Point<dim> &p,
679  * const unsigned int /*component*/ = 0) const override
680  * {
681  * Tensor<1, dim> r;
682  * r[0] = PI * std::cos(PI * p[0]) * std::sin(PI * p[1]);
683  * r[1] = PI * std::cos(PI * p[1]) * std::sin(PI * p[0]);
684  * return r;
685  * }
686  *
687  * virtual void
688  * hessian_list(const std::vector<Point<dim>> & points,
689  * std::vector<SymmetricTensor<2, dim>> &hessians,
690  * const unsigned int /*component*/ = 0) const override
691  * {
692  * for (unsigned i = 0; i < points.size(); ++i)
693  * {
694  * const double x = points[i][0];
695  * const double y = points[i][1];
696  *
697  * hessians[i][0][0] = -PI * PI * std::sin(PI * x) * std::sin(PI * y);
698  * hessians[i][0][1] = PI * PI * std::cos(PI * x) * std::cos(PI * y);
699  * hessians[i][1][1] = -PI * PI * std::sin(PI * x) * std::sin(PI * y);
700  * }
701  * }
702  * };
703  *
704  *
705  * template <int dim>
706  * class RightHandSide : public Function<dim>
707  * {
708  * public:
709  * static_assert(dim == 2, "Only dim==2 is implemented");
710  *
711  * virtual double value(const Point<dim> &p,
712  * const unsigned int /*component*/ = 0) const override
713  *
714  * {
715  * return 4 * std::pow(PI, 4.0) * std::sin(PI * p[0]) *
716  * std::sin(PI * p[1]);
717  * }
718  * };
719  * } // namespace ExactSolution
720  *
721  *
722  *
723  * @endcode
724  *
725  *
726  * <a name="Themainclass"></a>
727  * <h3>The main class</h3>
728  *
729 
730  *
731  * The following is the principal class of this tutorial program. It has
732  * the structure of many of the other tutorial programs and there should
733  * really be nothing particularly surprising about its contents or
734  * the constructor that follows it.
735  *
736  * @code
737  * template <int dim>
738  * class BiharmonicProblem
739  * {
740  * public:
741  * BiharmonicProblem(const unsigned int fe_degree);
742  *
743  * void run();
744  *
745  * private:
746  * void make_grid();
747  * void setup_system();
748  * void assemble_system();
749  * void solve();
750  * void compute_errors();
751  * void output_results(const unsigned int iteration) const;
752  *
754  *
755  * MappingQ<dim> mapping;
756  *
757  * FE_Q<dim> fe;
758  * DoFHandler<dim> dof_handler;
759  * AffineConstraints<double> constraints;
760  *
761  * SparsityPattern sparsity_pattern;
762  * SparseMatrix<double> system_matrix;
763  *
764  * Vector<double> solution;
765  * Vector<double> system_rhs;
766  * };
767  *
768  *
769  *
770  * template <int dim>
771  * BiharmonicProblem<dim>::BiharmonicProblem(const unsigned int fe_degree)
772  * : mapping(1)
773  * , fe(fe_degree)
774  * , dof_handler(triangulation)
775  * {}
776  *
777  *
778  *
779  * @endcode
780  *
781  * Next up are the functions that create the initial mesh (a once refined
782  * unit square) and set up the constraints, vectors, and matrices on
783  * each mesh. Again, both of these are essentially unchanged from many
784  * previous tutorial programs.
785  *
786  * @code
787  * template <int dim>
788  * void BiharmonicProblem<dim>::make_grid()
789  * {
791  * triangulation.refine_global(1);
792  *
793  * std::cout << "Number of active cells: " << triangulation.n_active_cells()
794  * << std::endl
795  * << "Total number of cells: " << triangulation.n_cells()
796  * << std::endl;
797  * }
798  *
799  *
800  *
801  * template <int dim>
802  * void BiharmonicProblem<dim>::setup_system()
803  * {
804  * dof_handler.distribute_dofs(fe);
805  *
806  * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
807  * << std::endl;
808  *
809  * constraints.clear();
810  * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
811  *
813  * 0,
814  * ExactSolution::Solution<dim>(),
815  * constraints);
816  * constraints.close();
817  *
818  *
819  * DynamicSparsityPattern dsp(dof_handler.n_dofs());
820  * DoFTools::make_flux_sparsity_pattern(dof_handler, dsp, constraints, true);
821  * sparsity_pattern.copy_from(dsp);
822  * system_matrix.reinit(sparsity_pattern);
823  *
824  * solution.reinit(dof_handler.n_dofs());
825  * system_rhs.reinit(dof_handler.n_dofs());
826  * }
827  *
828  *
829  *
830  * @endcode
831  *
832  *
833  * <a name="Assemblingthelinearsystem"></a>
834  * <h4>Assembling the linear system</h4>
835  *
836 
837  *
838  * The following pieces of code are more interesting. They all relate to the
839  * assembly of the linear system. While assembling the cell-interior terms
840  * is not of great difficulty -- that works in essence like the assembly
841  * of the corresponding terms of the Laplace equation, and you have seen
842  * how this works in @ref step_4 "step-4" or @ref step_6 "step-6", for example -- the difficulty
843  * is with the penalty terms in the formulation. These require the evaluation
844  * of gradients of shape functions at interfaces of cells. At the least,
845  * one would therefore need to use two FEFaceValues objects, but if one of the
846  * two sides is adaptively refined, then one actually needs an FEFaceValues
847  * and one FESubfaceValues objects; one also needs to keep track which
848  * shape functions live where, and finally we need to ensure that every
849  * face is visited only once. All of this is a substantial overhead to the
850  * logic we really want to implement (namely the penalty terms in the
851  * bilinear form). As a consequence, we will make use of the
852  * FEInterfaceValues class -- a helper class in deal.II that allows us
853  * to abstract away the two FEFaceValues or FESubfaceValues objects and
854  * directly access what we really care about: jumps, averages, etc.
855  *
856 
857  *
858  * But this doesn't yet solve our problem of having to keep track of
859  * which faces we have already visited when we loop over all cells and
860  * all of their faces. To make this process simpler, we use the
861  * MeshWorker::mesh_loop() function that provides a simple interface
862  * for this task: Based on the ideas outlined in the WorkStream
863  * namespace documentation, MeshWorker::mesh_loop() requires three
864  * functions that do work on cells, interior faces, and boundary
865  * faces. These functions work on scratch objects for intermediate
866  * results, and then copy the result of their computations into
867  * copy data objects from where a copier function copies them into
868  * the global matrix and right hand side objects.
869  *
870 
871  *
872  * The following structures then provide the scratch and copy objects
873  * that are necessary for this approach. You may look up the WorkStream
874  * namespace as well as the
875  * @ref threads "Parallel computing with multiple processors"
876  * module for more information on how they typically work.
877  *
878  * @code
879  * template <int dim>
880  * struct ScratchData
881  * {
882  * ScratchData(const Mapping<dim> & mapping,
883  * const FiniteElement<dim> &fe,
884  * const unsigned int quadrature_degree,
885  * const UpdateFlags update_flags,
886  * const UpdateFlags interface_update_flags)
887  * : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
888  * , fe_interface_values(mapping,
889  * fe,
890  * QGauss<dim - 1>(quadrature_degree),
891  * interface_update_flags)
892  * {}
893  *
894  *
895  * ScratchData(const ScratchData<dim> &scratch_data)
896  * : fe_values(scratch_data.fe_values.get_mapping(),
897  * scratch_data.fe_values.get_fe(),
898  * scratch_data.fe_values.get_quadrature(),
899  * scratch_data.fe_values.get_update_flags())
900  * , fe_interface_values(scratch_data.fe_values.get_mapping(),
901  * scratch_data.fe_values.get_fe(),
902  * scratch_data.fe_interface_values.get_quadrature(),
903  * scratch_data.fe_interface_values.get_update_flags())
904  * {}
905  *
906  * FEValues<dim> fe_values;
907  * FEInterfaceValues<dim> fe_interface_values;
908  * };
909  *
910  *
911  *
912  * struct CopyData
913  * {
914  * CopyData(const unsigned int dofs_per_cell)
915  * : cell_matrix(dofs_per_cell, dofs_per_cell)
916  * , cell_rhs(dofs_per_cell)
917  * , local_dof_indices(dofs_per_cell)
918  * {}
919  *
920  *
921  * CopyData(const CopyData &) = default;
922  *
923  *
924  * struct FaceData
925  * {
926  * FullMatrix<double> cell_matrix;
927  * std::vector<types::global_dof_index> joint_dof_indices;
928  * };
929  *
930  * FullMatrix<double> cell_matrix;
931  * Vector<double> cell_rhs;
932  * std::vector<types::global_dof_index> local_dof_indices;
933  * std::vector<FaceData> face_data;
934  * };
935  *
936  *
937  *
938  * @endcode
939  *
940  * The more interesting part is where we actually assemble the linear system.
941  * Fundamentally, this function has five parts:
942  * - The definition of the `cell_worker` lambda function, a small
943  * function that is defined within the `assemble_system()`
944  * function and that will be responsible for computing the local
945  * integrals on an individual cell. It will work on a copy of the
946  * `ScratchData` class and put its results into the corresponding
947  * `CopyData` object.
948  * - The definition of the `face_worker` lambda function that does
949  * the integration of all terms that live on the interfaces between
950  * cells.
951  * - The definition of the `boundary_worker` function that does the
952  * same but for cell faces located on the boundary of the domain.
953  * - The definition of the `copier` function that is responsible
954  * for copying all of the data the previous three functions have
955  * put into copy objects for a single cell, into the global matrix
956  * and right hand side.
957  *
958 
959  *
960  * The fifth part is the one where we bring all of this together.
961  *
962 
963  *
964  * Let us go through each of these pieces necessary for the assembly
965  * in turns.
966  *
967  * @code
968  * template <int dim>
969  * void BiharmonicProblem<dim>::assemble_system()
970  * {
971  * using Iterator = typename DoFHandler<dim>::active_cell_iterator;
972  *
973  * @endcode
974  *
975  * The first piece is the `cell_worker` that does the assembly
976  * on the cell interiors. It is a (lambda) function that takes
977  * a cell (input), a scratch object, and a copy object (output)
978  * as arguments. It looks like the assembly functions of many
979  * other of the tutorial programs, or at least the body of the
980  * loop over all cells.
981  *
982 
983  *
984  * The terms we integrate here are the cell contribution
985  * @f{align*}{
986  * A^K_{ij} = \int_K \nabla^2\varphi_i(x) : \nabla^2\varphi_j(x) dx
987  * @f}
988  * to the global matrix, and
989  * @f{align*}{
990  * f^K_i = \int_K \varphi_i(x) f(x) dx
991  * @f}
992  * to the right hand side vector.
993  *
994 
995  *
996  * We use the same technique as used in the assembly of @ref step_22 "step-22"
997  * to accelerate the function: Instead of calling
998  * `fe_values.shape_hessian(i, qpoint)` in the innermost loop,
999  * we instead create a variable `hessian_i` that evaluates this
1000  * value once in the loop over `i` and re-use the so-evaluated
1001  * value in the loop over `j`. For symmetry, we do the same with a
1002  * variable `hessian_j`, although it is indeed only used once and
1003  * we could have left the call to `fe_values.shape_hessian(j,qpoint)`
1004  * in the instruction that computes the scalar product between
1005  * the two terms.
1006  *
1007  * @code
1008  * auto cell_worker = [&](const Iterator & cell,
1009  * ScratchData<dim> &scratch_data,
1010  * CopyData & copy_data) {
1011  * copy_data.cell_matrix = 0;
1012  * copy_data.cell_rhs = 0;
1013  *
1014  * FEValues<dim> &fe_values = scratch_data.fe_values;
1015  * fe_values.reinit(cell);
1016  *
1017  * cell->get_dof_indices(copy_data.local_dof_indices);
1018  *
1019  * const ExactSolution::RightHandSide<dim> right_hand_side;
1020  *
1021  * const unsigned int dofs_per_cell =
1022  * scratch_data.fe_values.get_fe().dofs_per_cell;
1023  *
1024  * for (unsigned int qpoint = 0; qpoint < fe_values.n_quadrature_points;
1025  * ++qpoint)
1026  * {
1027  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1028  * {
1029  * const Tensor<2, dim> hessian_i =
1030  * fe_values.shape_hessian(i, qpoint);
1031  *
1032  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1033  * {
1034  * const Tensor<2, dim> hessian_j =
1035  * fe_values.shape_hessian(j, qpoint);
1036  *
1037  * copy_data.cell_matrix(i, j) +=
1038  * scalar_product(hessian_i, // nabla^2 phi_i(x)
1039  * hessian_j) * // nabla^2 phi_j(x)
1040  * fe_values.JxW(qpoint); // dx
1041  * }
1042  *
1043  * copy_data.cell_rhs(i) +=
1044  * fe_values.shape_value(i, qpoint) * // phi_i(x)
1045  * right_hand_side.value(
1046  * fe_values.quadrature_point(qpoint)) * // f(x)
1047  * fe_values.JxW(qpoint); // dx
1048  * }
1049  * }
1050  * };
1051  *
1052  *
1053  * @endcode
1054  *
1055  * The next building block is the one that assembles penalty terms on each
1056  * of the interior faces of the mesh. As described in the documentation of
1057  * MeshWorker::mesh_loop(), this function receives arguments that denote
1058  * a cell and its neighboring cell, as well as (for each of the two
1059  * cells) the face (and potentially sub-face) we have to integrate
1060  * over. Again, we also get a scratch object, and a copy object
1061  * for putting the results in.
1062  *
1063 
1064  *
1065  * The function has three parts itself. At the top, we initialize
1066  * the FEInterfaceValues object and create a new `CopyData::FaceData`
1067  * object to store our input in. This gets pushed to the end of the
1068  * `copy_data.face_data` variable. We need to do this because
1069  * the number of faces (or subfaces) over which we integrate for a
1070  * given cell differs from cell to cell, and the sizes of these
1071  * matrices also differ, depending on what degrees of freedom
1072  * are adjacent to the face or subface. As discussed in the documentation
1073  * of MeshWorker::mesh_loop(), the copy object is reset every time a new
1074  * cell is visited, so that what we push to the end of
1075  * `copy_data.face_data()` is really all that the later `copier` function
1076  * gets to see when it copies the contributions of each cell to the global
1077  * matrix and right hand side objects.
1078  *
1079  * @code
1080  * auto face_worker = [&](const Iterator & cell,
1081  * const unsigned int &f,
1082  * const unsigned int &sf,
1083  * const Iterator & ncell,
1084  * const unsigned int &nf,
1085  * const unsigned int &nsf,
1086  * ScratchData<dim> & scratch_data,
1087  * CopyData & copy_data) {
1088  * FEInterfaceValues<dim> &fe_interface_values =
1089  * scratch_data.fe_interface_values;
1090  * fe_interface_values.reinit(cell, f, sf, ncell, nf, nsf);
1091  *
1092  * copy_data.face_data.emplace_back();
1093  * CopyData::FaceData &copy_data_face = copy_data.face_data.back();
1094  *
1095  * copy_data_face.joint_dof_indices =
1096  * fe_interface_values.get_interface_dof_indices();
1097  *
1098  * const unsigned int n_interface_dofs =
1099  * fe_interface_values.n_current_interface_dofs();
1100  * copy_data_face.cell_matrix.reinit(n_interface_dofs, n_interface_dofs);
1101  *
1102  * @endcode
1103  *
1104  * The second part deals with determining what the penalty
1105  * parameter should be. By looking at the units of the various
1106  * terms in the bilinear form, it is clear that the penalty has
1107  * to have the form @f$\frac{\gamma}{h_K}@f$ (i.e., one over length
1108  * scale), but it is not a priori obvious how one should choose
1109  * the dimension-less number @f$\gamma@f$. From the discontinuous
1110  * Galerkin theory for the Laplace equation, one might
1111  * conjecture that the right choice is @f$\gamma=p(p+1)@f$ is the
1112  * right choice, where @f$p@f$ is the polynomial degree of the
1113  * finite element used. We will discuss this choice in a bit
1114  * more detail in the results section of this program.
1115  *
1116 
1117  *
1118  * In the formula above, @f$h_K@f$ is the size of cell @f$K@f$. But this
1119  * is not quite so straightforward either: If one uses highly
1120  * stretched cells, then a more involved theory says that @f$h@f$
1121  * should be replaced by the diameter of cell @f$K@f$ normal to the
1122  * direction of the edge in question. It turns out that there
1123  * is a function in deal.II for that. Secondly, @f$h_K@f$ may be
1124  * different when viewed from the two different sides of a face.
1125  *
1126 
1127  *
1128  * To stay on the safe side, we take the maximum of the two values.
1129  * We will note that it is possible that this computation has to be
1130  * further adjusted if one were to use hanging nodes resulting from
1131  * adaptive mesh refinement.
1132  *
1133  * @code
1134  * const unsigned int p = fe.degree;
1135  * const double gamma_over_h =
1136  * std::max((1.0 * p * (p + 1) /
1137  * cell->extent_in_direction(
1138  * GeometryInfo<dim>::unit_normal_direction[f])),
1139  * (1.0 * p * (p + 1) /
1140  * ncell->extent_in_direction(
1141  * GeometryInfo<dim>::unit_normal_direction[nf])));
1142  *
1143  * @endcode
1144  *
1145  * Finally, and as usual, we loop over the quadrature points and
1146  * indices `i` and `j` to add up the contributions of this face
1147  * or sub-face. These are then stored in the
1148  * `copy_data.face_data` object created above. As for the cell
1149  * worker, we pull the evaluation of averages and jumps out of
1150  * the loops if possible, introducing local variables that store
1151  * these results. The assembly then only needs to use these
1152  * local variables in the innermost loop. Regarding the concrete
1153  * formula this code implements, recall that the interface terms
1154  * of the bilinear form were as follows:
1155  * @f{align*}{
1156  * -\sum_{e \in \mathbb{F}} \int_{e}
1157  * \jump{ \frac{\partial v_h}{\partial \mathbf n}}
1158  * \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds
1159  * -\sum_{e \in \mathbb{F}} \int_{e}
1160  * \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}}
1161  * \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds
1162  * + \sum_{e \in \mathbb{F}}
1163  * \frac{\gamma}{h_e}
1164  * \int_e
1165  * \jump{\frac{\partial v_h}{\partial \mathbf n}}
1166  * \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds.
1167  * @f}
1168  *
1169  * @code
1170  * for (unsigned int qpoint = 0;
1171  * qpoint < fe_interface_values.n_quadrature_points;
1172  * ++qpoint)
1173  * {
1174  * const auto &n = fe_interface_values.normal(qpoint);
1175  *
1176  * for (unsigned int i = 0; i < n_interface_dofs; ++i)
1177  * {
1178  * const double av_hessian_i_dot_n_dot_n =
1179  * (fe_interface_values.average_hessian(i, qpoint) * n * n);
1180  * const double jump_grad_i_dot_n =
1181  * (fe_interface_values.jump_gradient(i, qpoint) * n);
1182  *
1183  * for (unsigned int j = 0; j < n_interface_dofs; ++j)
1184  * {
1185  * const double av_hessian_j_dot_n_dot_n =
1186  * (fe_interface_values.average_hessian(j, qpoint) * n * n);
1187  * const double jump_grad_j_dot_n =
1188  * (fe_interface_values.jump_gradient(j, qpoint) * n);
1189  *
1190  * copy_data_face.cell_matrix(i, j) +=
1191  * (-av_hessian_i_dot_n_dot_n // - {grad^2 v n n }
1192  * * jump_grad_j_dot_n // [grad u n]
1193  * - av_hessian_j_dot_n_dot_n // - {grad^2 u n n }
1194  * * jump_grad_i_dot_n // [grad v n]
1195  * + // +
1196  * gamma_over_h * // gamma/h
1197  * jump_grad_i_dot_n * // [grad v n]
1198  * jump_grad_j_dot_n) * // [grad u n]
1199  * fe_interface_values.JxW(qpoint); // dx
1200  * }
1201  * }
1202  * }
1203  * };
1204  *
1205  *
1206  * @endcode
1207  *
1208  * The third piece is to do the same kind of assembly for faces that
1209  * are at the boundary. The idea is the same as above, of course,
1210  * with only the difference that there are now penalty terms that
1211  * also go into the right hand side.
1212  *
1213 
1214  *
1215  * As before, the first part of the function simply sets up some
1216  * helper objects:
1217  *
1218  * @code
1219  * auto boundary_worker = [&](const Iterator & cell,
1220  * const unsigned int &face_no,
1221  * ScratchData<dim> & scratch_data,
1222  * CopyData & copy_data) {
1223  * FEInterfaceValues<dim> &fe_interface_values =
1224  * scratch_data.fe_interface_values;
1225  * fe_interface_values.reinit(cell, face_no);
1226  * const auto &q_points = fe_interface_values.get_quadrature_points();
1227  *
1228  * copy_data.face_data.emplace_back();
1229  * CopyData::FaceData &copy_data_face = copy_data.face_data.back();
1230  *
1231  * const unsigned int n_dofs =
1232  * fe_interface_values.n_current_interface_dofs();
1233  * copy_data_face.joint_dof_indices =
1234  * fe_interface_values.get_interface_dof_indices();
1235  *
1236  * copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
1237  *
1238  * const std::vector<double> &JxW = fe_interface_values.get_JxW_values();
1239  * const std::vector<Tensor<1, dim>> &normals =
1240  * fe_interface_values.get_normal_vectors();
1241  *
1242  *
1243  * const ExactSolution::Solution<dim> exact_solution;
1244  * std::vector<Tensor<1, dim>> exact_gradients(q_points.size());
1245  * exact_solution.gradient_list(q_points, exact_gradients);
1246  *
1247  *
1248  * @endcode
1249  *
1250  * Positively, because we now only deal with one cell adjacent to the
1251  * face (as we are on the boundary), the computation of the penalty
1252  * factor @f$\gamma@f$ is substantially simpler:
1253  *
1254  * @code
1255  * const unsigned int p = fe.degree;
1256  * const double gamma_over_h =
1257  * (1.0 * p * (p + 1) /
1258  * cell->extent_in_direction(
1259  * GeometryInfo<dim>::unit_normal_direction[face_no]));
1260  *
1261  * @endcode
1262  *
1263  * The third piece is the assembly of terms. This is now
1264  * slightly more involved since these contains both terms for
1265  * the matrix and for the right hand side. The former is exactly
1266  * the same as for the interior faces stated above if one just
1267  * defines the jump and average appropriately (which is what the
1268  * FEInterfaceValues class does). The latter requires us to
1269  * evaluate the boundary conditions @f$j(\mathbf x)@f$, which in the
1270  * current case (where we know the exact solution) we compute
1271  * from @f$j(\mathbf x) = \frac{\partial u(\mathbf x)}{\partial
1272  * {\mathbf n}}@f$. The term to be added to the right hand side
1273  * vector is then
1274  * @f$\frac{\gamma}{h_e}\int_e
1275  * \jump{\frac{\partial v_h}{\partial \mathbf n}} j \ ds@f$.
1276  *
1277  * @code
1278  * for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
1279  * {
1280  * const auto &n = normals[qpoint];
1281  *
1282  * for (unsigned int i = 0; i < n_dofs; ++i)
1283  * {
1284  * const double av_hessian_i_dot_n_dot_n =
1285  * (fe_interface_values.average_hessian(i, qpoint) * n * n);
1286  * const double jump_grad_i_dot_n =
1287  * (fe_interface_values.jump_gradient(i, qpoint) * n);
1288  *
1289  * for (unsigned int j = 0; j < n_dofs; ++j)
1290  * {
1291  * const double av_hessian_j_dot_n_dot_n =
1292  * (fe_interface_values.average_hessian(j, qpoint) * n * n);
1293  * const double jump_grad_j_dot_n =
1294  * (fe_interface_values.jump_gradient(j, qpoint) * n);
1295  *
1296  * copy_data_face.cell_matrix(i, j) +=
1297  * (-av_hessian_i_dot_n_dot_n // - {grad^2 v n n}
1298  * * jump_grad_j_dot_n // [grad u n]
1299  *
1300  * - av_hessian_j_dot_n_dot_n // - {grad^2 u n n}
1301  * * jump_grad_i_dot_n // [grad v n]
1302  *
1303  * + gamma_over_h // gamma/h
1304  * * jump_grad_i_dot_n // [grad v n]
1305  * * jump_grad_j_dot_n // [grad u n]
1306  * ) *
1307  * JxW[qpoint]; // dx
1308  * }
1309  *
1310  * copy_data.cell_rhs(i) +=
1311  * (-av_hessian_i_dot_n_dot_n * // - {grad^2 v n n }
1312  * (exact_gradients[qpoint] * n) // (grad u_exact . n)
1313  * + // +
1314  * gamma_over_h // gamma/h
1315  * * jump_grad_i_dot_n // [grad v n]
1316  * * (exact_gradients[qpoint] * n) // (grad u_exact . n)
1317  * ) *
1318  * JxW[qpoint]; // dx
1319  * }
1320  * }
1321  * };
1322  *
1323  * @endcode
1324  *
1325  * Part 4 is a small function that copies the data produced by the
1326  * cell, interior, and boundary face assemblers above into the
1327  * global matrix and right hand side vector. There really is not
1328  * very much to do here: We distribute the cell matrix and right
1329  * hand side contributions as we have done in almost all of the
1330  * other tutorial programs using the constraints objects. We then
1331  * also have to do the same for the face matrix contributions
1332  * that have gained content for the faces (interior and boundary)
1333  * and that the `face_worker` and `boundary_worker` have added
1334  * to the `copy_data.face_data` array.
1335  *
1336  * @code
1337  * auto copier = [&](const CopyData &copy_data) {
1338  * constraints.distribute_local_to_global(copy_data.cell_matrix,
1339  * copy_data.cell_rhs,
1340  * copy_data.local_dof_indices,
1341  * system_matrix,
1342  * system_rhs);
1343  *
1344  * for (auto &cdf : copy_data.face_data)
1345  * {
1346  * constraints.distribute_local_to_global(cdf.cell_matrix,
1347  * cdf.joint_dof_indices,
1348  * system_matrix);
1349  * }
1350  * };
1351  *
1352  *
1353  * @endcode
1354  *
1355  * Having set all of this up, what remains is to just create a scratch
1356  * and copy data object and call the MeshWorker::mesh_loop() function
1357  * that then goes over all cells and faces, calls the respective workers
1358  * on them, and then the copier function that puts things into the
1359  * global matrix and right hand side. As an additional benefit,
1360  * MeshWorker::mesh_loop() does all of this in parallel, using
1361  * as many processor cores as your machine happens to have.
1362  *
1363  * @code
1364  * const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
1365  * ScratchData<dim> scratch_data(mapping,
1366  * fe,
1367  * n_gauss_points,
1368  * update_values | update_gradients |
1369  * update_hessians | update_quadrature_points |
1370  * update_JxW_values,
1371  * update_values | update_gradients |
1372  * update_hessians | update_quadrature_points |
1373  * update_JxW_values | update_normal_vectors);
1374  * CopyData copy_data(dof_handler.get_fe().dofs_per_cell);
1375  * MeshWorker::mesh_loop(dof_handler.begin_active(),
1376  * dof_handler.end(),
1377  * cell_worker,
1378  * copier,
1379  * scratch_data,
1380  * copy_data,
1381  * MeshWorker::assemble_own_cells |
1382  * MeshWorker::assemble_boundary_faces |
1383  * MeshWorker::assemble_own_interior_faces_once,
1384  * boundary_worker,
1385  * face_worker);
1386  * }
1387  *
1388  *
1389  *
1390  * @endcode
1391  *
1392  *
1393  * <a name="Solvingthelinearsystemandpostprocessing"></a>
1394  * <h4>Solving the linear system and postprocessing</h4>
1395  *
1396 
1397  *
1398  * The show is essentially over at this point: The remaining functions are
1399  * not overly interesting or novel. The first one simply uses a direct
1400  * solver to solve the linear system (see also @ref step_29 "step-29"):
1401  *
1402  * @code
1403  * template <int dim>
1404  * void BiharmonicProblem<dim>::solve()
1405  * {
1406  * std::cout << " Solving system..." << std::endl;
1407  *
1408  * SparseDirectUMFPACK A_direct;
1409  * A_direct.initialize(system_matrix);
1410  * A_direct.vmult(solution, system_rhs);
1411  *
1412  * constraints.distribute(solution);
1413  * }
1414  *
1415  *
1416  *
1417  * @endcode
1418  *
1419  * The next function evaluates the error between the computed solution
1420  * and the exact solution (which is known here because we have chosen
1421  * the right hand side and boundary values in a way so that we know
1422  * the corresponding solution). In the first two code blocks below,
1423  * we compute the error in the @f$L_2@f$ norm and the @f$H^1@f$ semi-norm.
1424  *
1425  * @code
1426  * template <int dim>
1427  * void BiharmonicProblem<dim>::compute_errors()
1428  * {
1429  * {
1430  * Vector<float> norm_per_cell(triangulation.n_active_cells());
1431  * VectorTools::integrate_difference(mapping,
1432  * dof_handler,
1433  * solution,
1434  * ExactSolution::Solution<dim>(),
1435  * norm_per_cell,
1436  * QGauss<dim>(fe.degree + 2),
1437  * VectorTools::L2_norm);
1438  * const double error_norm =
1439  * VectorTools::compute_global_error(triangulation,
1440  * norm_per_cell,
1441  * VectorTools::L2_norm);
1442  * std::cout << " Error in the L2 norm : " << error_norm
1443  * << std::endl;
1444  * }
1445  *
1446  * {
1447  * Vector<float> norm_per_cell(triangulation.n_active_cells());
1448  * VectorTools::integrate_difference(mapping,
1449  * dof_handler,
1450  * solution,
1451  * ExactSolution::Solution<dim>(),
1452  * norm_per_cell,
1453  * QGauss<dim>(fe.degree + 2),
1454  * VectorTools::H1_seminorm);
1455  * const double error_norm =
1456  * VectorTools::compute_global_error(triangulation,
1457  * norm_per_cell,
1458  * VectorTools::H1_seminorm);
1459  * std::cout << " Error in the H1 seminorm : " << error_norm
1460  * << std::endl;
1461  * }
1462  *
1463  * @endcode
1464  *
1465  * Now also compute an approximation to the @f$H^2@f$ seminorm error. The actual
1466  * @f$H^2@f$ seminorm would require us to integrate second derivatives of the
1467  * solution @f$u_h@f$, but given the Lagrange shape functions we use, @f$u_h@f$ of
1468  * course has kinks at the interfaces between cells, and consequently second
1469  * derivatives are singular at interfaces. As a consequence, we really only
1470  * integrate over the interior of cells and ignore the interface
1471  * contributions. This is *not* an equivalent norm to the energy norm for
1472  * the problem, but still gives us an idea of how fast the error converges.
1473  *
1474 
1475  *
1476  * We note that one could address this issue by defining a norm that
1477  * is equivalent to the energy norm. This would involve adding up not
1478  * only the integrals over cell interiors as we do below, but also adding
1479  * penalty terms for the jump of the derivative of @f$u_h@f$ across interfaces,
1480  * with an appropriate scaling of the two kinds of terms. We will leave
1481  * this for later work.
1482  *
1483  * @code
1484  * {
1485  * const QGauss<dim> quadrature_formula(fe.degree + 2);
1486  * ExactSolution::Solution<dim> exact_solution;
1487  * Vector<double> error_per_cell(triangulation.n_active_cells());
1488  *
1489  * FEValues<dim> fe_values(mapping,
1490  * fe,
1491  * quadrature_formula,
1492  * update_values | update_hessians |
1493  * update_quadrature_points | update_JxW_values);
1494  *
1495  * FEValuesExtractors::Scalar scalar(0);
1496  * const unsigned int n_q_points = quadrature_formula.size();
1497  *
1498  * std::vector<SymmetricTensor<2, dim>> exact_hessians(n_q_points);
1499  * std::vector<Tensor<2, dim>> hessians(n_q_points);
1500  * for (auto &cell : dof_handler.active_cell_iterators())
1501  * {
1502  * fe_values.reinit(cell);
1503  * fe_values[scalar].get_function_hessians(solution, hessians);
1504  * exact_solution.hessian_list(fe_values.get_quadrature_points(),
1505  * exact_hessians);
1506  *
1507  * double local_error = 0;
1508  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1509  * {
1510  * local_error +=
1511  * ((exact_hessians[q_point] - hessians[q_point]).norm_square() *
1512  * fe_values.JxW(q_point));
1513  * }
1514  * error_per_cell[cell->active_cell_index()] = std::sqrt(local_error);
1515  * }
1516  *
1517  * const double error_norm = error_per_cell.l2_norm();
1518  * std::cout << " Error in the broken H2 seminorm: " << error_norm
1519  * << std::endl;
1520  * }
1521  * }
1522  *
1523  *
1524  *
1525  * @endcode
1526  *
1527  * Equally uninteresting is the function that generates graphical output.
1528  * It looks exactly like the one in @ref step_6 "step-6", for example.
1529  *
1530  * @code
1531  * template <int dim>
1532  * void
1533  * BiharmonicProblem<dim>::output_results(const unsigned int iteration) const
1534  * {
1535  * std::cout << " Writing graphical output..." << std::endl;
1536  *
1537  * DataOut<dim> data_out;
1538  *
1539  * data_out.attach_dof_handler(dof_handler);
1540  * data_out.add_data_vector(solution, "solution");
1541  * data_out.build_patches();
1542  *
1543  * std::ofstream output_vtu(
1544  * ("output_" + Utilities::int_to_string(iteration, 6) + ".vtu").c_str());
1545  * data_out.write_vtu(output_vtu);
1546  * }
1547  *
1548  *
1549  *
1550  * @endcode
1551  *
1552  * The same is true for the `run()` function: Just like in previous
1553  * programs.
1554  *
1555  * @code
1556  * template <int dim>
1557  * void BiharmonicProblem<dim>::run()
1558  * {
1559  * make_grid();
1560  *
1561  * const unsigned int n_cycles = 4;
1562  * for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1563  * {
1564  * std::cout << "Cycle: " << cycle << " of " << n_cycles << std::endl;
1565  *
1566  * triangulation.refine_global(1);
1567  * setup_system();
1568  *
1569  * assemble_system();
1570  * solve();
1571  *
1572  * output_results(cycle);
1573  *
1574  * compute_errors();
1575  * std::cout << std::endl;
1576  * }
1577  * }
1578  * } // namespace Step47
1579  *
1580  *
1581  *
1582  * @endcode
1583  *
1584  *
1585  * <a name="Themainfunction"></a>
1586  * <h3>The main() function</h3>
1587  *
1588 
1589  *
1590  * Finally for the `main()` function. There is, again, not very much to see
1591  * here: It looks like the ones in previous tutorial programs. There
1592  * is a variable that allows selecting the polynomial degree of the element
1593  * we want to use for solving the equation. Because the C0IP formulation
1594  * we use requires the element degree to be at least two, we check with
1595  * an assertion that whatever one sets for the polynomial degree actually
1596  * makes sense.
1597  *
1598  * @code
1599  * int main()
1600  * {
1601  * try
1602  * {
1603  * using namespace dealii;
1604  * using namespace Step47;
1605  *
1606  * const unsigned int fe_degree = 2;
1607  * Assert(fe_degree >= 2,
1608  * ExcMessage("The C0IP formulation for the biharmonic problem "
1609  * "only works if one uses elements of polynomial "
1610  * "degree at least 2."));
1611  *
1612  * BiharmonicProblem<2> biharmonic_problem(fe_degree);
1613  * biharmonic_problem.run();
1614  * }
1615  * catch (std::exception &exc)
1616  * {
1617  * std::cerr << std::endl
1618  * << std::endl
1619  * << "----------------------------------------------------"
1620  * << std::endl;
1621  * std::cerr << "Exception on processing: " << std::endl
1622  * << exc.what() << std::endl
1623  * << "Aborting!" << std::endl
1624  * << "----------------------------------------------------"
1625  * << std::endl;
1626  *
1627  * return 1;
1628  * }
1629  * catch (...)
1630  * {
1631  * std::cerr << std::endl
1632  * << std::endl
1633  * << "----------------------------------------------------"
1634  * << std::endl;
1635  * std::cerr << "Unknown exception!" << std::endl
1636  * << "Aborting!" << std::endl
1637  * << "----------------------------------------------------"
1638  * << std::endl;
1639  * return 1;
1640  * }
1641  *
1642  * return 0;
1643  * }
1644  * @endcode
1645 <a name="Results"></a><h1>Results</h1>
1646 
1647 
1648 We run the program with right hand side and boundary values as
1649 discussed in the introduction. These will produce the
1650 solution @f$u = \sin(\pi x) \sin(\pi y)@f$ on the domain @f$\Omega = (0,1)^2@f$.
1651 We test this setup using @f$Q_2@f$, @f$Q_3@f$, and @f$Q_4@f$ elements, which one can
1652 change via the `fe_degree` variable in the `main()` function. With mesh
1653 refinement, the @f$L_2@f$ convergence rates, @f$H^1@f$-seminorm rate,
1654 and @f$H^2@f$-seminorm convergence of @f$u@f$
1655 should then be around 2, 2, 1 for @f$Q_2@f$ (with the @f$L_2@f$ norm
1656 sub-optimal as discussed in the introduction); 4, 3, 2 for
1657 @f$Q_3@f$; and 5, 4, 3 for @f$Q_4@f$, respectively.
1658 
1659 From the literature, it is not immediately clear what
1660 the penalty parameter @f$\gamma@f$ should be. For example,
1661 @cite Brenner2009 state that it needs to be larger than one, and
1662 choose @f$\gamma=5@f$. The FEniCS/Dolphin tutorial chooses it as
1663 @f$\gamma=8@f$, see
1664 https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html
1665 . @cite Wells2007 uses a value for @f$\gamma@f$ larger than the
1666 number of edges belonging to an element for Kirchhoff plates (see
1667 their Section 4.2). This suggests that maybe
1668 @f$\gamma = 1@f$, @f$2@f$, are too small; on the other hand, a value
1669 @f$p(p+1)@f$ would be reasonable,
1670 where @f$p@f$ is the degree of polynomials. The last of these choices is
1671 the one one would expect to work by comparing
1672 to the discontinuous Galerkin formulations for the Laplace equation,
1673 and it will turn out to also work here.
1674 But we should check what value of @f$\gamma@f$ is right, and we will do so
1675 below; changing @f$\gamma@f$ is easy in the two `face_worker` and
1676 `boundary_worker` functions defined in `assemble_system()`.
1677 
1678 
1679 <a name="TestresultsoniQsub2subiwithigammapp1i"></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>&gamma; = p(p+1)</i> </h3>
1680 
1681 
1682 We run the code with differently refined meshes
1683 and get the following convergence rates.
1684 
1685 <table align="center" class="doxtable">
1686  <tr>
1687  <th>Number of refinements </th><th> @f$\|u-u_h^\circ\|_{L_2}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^1}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^2}@f$ </th><th> Conv. rates </th>
1688  </tr>
1689  <tr>
1690  <td> 2 </td><td> 8.780e-03 </td><td> </td><td> 7.095e-02 </td><td> </td><td> 1.645 </td><td> </td>
1691  </tr>
1692  <tr>
1693  <td> 3 </td><td> 3.515e-03 </td><td> 1.32 </td><td> 2.174e-02 </td><td> 1.70 </td><td> 8.121e-01 </td><td> 1.018 </td>
1694  </tr>
1695  <tr>
1696  <td> 4 </td><td> 1.103e-03 </td><td> 1.67 </td><td> 6.106e-03 </td><td> 1.83 </td><td> 4.015e-01 </td><td> 1.016 </td>
1697  </tr>
1698  <tr>
1699  <td> 5 </td><td> 3.084e-04 </td><td> 1.83 </td><td> 1.622e-03 </td><td> 1.91 </td><td> 1.993e-01 </td><td> 1.010 </td>
1700  </tr>
1701 </table>
1702 We can see that the @f$L_2@f$ convergence rates are around 2,
1703 @f$H^1@f$-seminorm convergence rates are around 2,
1704 and @f$H^2@f$-seminorm convergence rates are around 1. The latter two
1705 match the theoretically expected rates; for the former, we have no
1706 theorem but are not surprised that it is sub-optimal given the remark
1707 in the introduction.
1708 
1709 
1710 <a name="TestresultsoniQsub3subiwithigammapp1i"></a><h3>Test results on <i>Q<sub>3</sub></i> with <i>&gamma; = p(p+1)</i> </h3>
1711 
1712 
1713 
1714 <table align="center" class="doxtable">
1715  <tr>
1716  <th>Number of refinements </th><th> @f$\|u-u_h^\circ\|_{L_2}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^1}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^2}@f$ </th><th> Conv. rates </th>
1717  </tr>
1718  <tr>
1719  <td> 2 </td><td> 2.045e-04 </td><td> </td><td> 4.402e-03 </td><td> </td><td> 1.641e-01 </td><td> </td>
1720  </tr>
1721  <tr>
1722  <td> 3 </td><td> 1.312e-05 </td><td> 3.96 </td><td> 5.537e-04 </td><td> 2.99 </td><td> 4.096e-02 </td><td> 2.00 </td>
1723  </tr>
1724  <tr>
1725  <td> 4 </td><td> 8.239e-07 </td><td> 3.99 </td><td> 6.904e-05 </td><td> 3.00 </td><td> 1.023e-02 </td><td> 2.00 </td>
1726  </tr>
1727  <tr>
1728  <td> 5 </td><td> 5.158e-08 </td><td> 3.99 </td><td> 8.621e-06 </td><td> 3.00 </td><td> 2.558e-03 </td><td> 2.00 </td>
1729  </tr>
1730 </table>
1731 We can see that the @f$L_2@f$ convergence rates are around 4,
1732 @f$H^1@f$-seminorm convergence rates are around 3,
1733 and @f$H^2@f$-seminorm convergence rates are around 2.
1734 This, of course, matches our theoretical expectations.
1735 
1736 
1737 <a name="TestresultsoniQsub4subiwithigammapp1i"></a><h3>Test results on <i>Q<sub>4</sub></i> with <i>&gamma; = p(p+1)</i> </h3>
1738 
1739 
1740 <table align="center" class="doxtable">
1741  <tr>
1742  <th>Number of refinements </th><th> @f$\|u-u_h^\circ\|_{L_2}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^1}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^2}@f$ </th><th> Conv. rates </th>
1743  </tr>
1744  <tr>
1745  <td> 2 </td><td> 6.510e-06 </td><td> </td><td> 2.215e-04 </td><td> </td><td> 1.275e-02 </td><td> </td>
1746  </tr>
1747  <tr>
1748  <td> 3 </td><td> 2.679e-07 </td><td> 4.60 </td><td> 1.569e-05 </td><td> 3.81 </td><td> 1.496e-03 </td><td> 3.09 </td>
1749  </tr>
1750  <tr>
1751  <td> 4 </td><td> 9.404e-09 </td><td> 4.83 </td><td> 1.040e-06 </td><td> 3.91 </td><td> 1.774e-04 </td><td> 3.07 </td>
1752  </tr>
1753  <tr>
1754  <td> 5 </td><td> 7.943e-10 </td><td> 3.56 </td><td> 6.693e-08 </td><td> 3.95 </td><td> 2.150e-05 </td><td> 3.04 </td>
1755  </tr>
1756 </table>
1757 We can see that the @f$L_2@f$ norm convergence rates are around 5,
1758 @f$H^1@f$-seminorm convergence rates are around 4,
1759 and @f$H^2@f$-seminorm convergence rates are around 3.
1760 On the finest mesh, the @f$L_2@f$ norm convergence rate
1761 is much smaller than our theoretical expectations
1762 because the linear solver becomes the limiting factor due
1763 to round-off. Of course the @f$L_2@f$ error is also very small already in
1764 that case.
1765 
1766 
1767 <a name="TestresultsoniQsub2subiwithigamma1i"></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>&gamma; = 1</i> </h3>
1768 
1769 
1770 For comparison with the results above, let us now also consider the
1771 case where we simply choose @f$\gamma=1@f$:
1772 
1773 <table align="center" class="doxtable">
1774  <tr>
1775  <th>Number of refinements </th><th> @f$\|u-u_h^\circ\|_{L_2}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^1}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^2}@f$ </th><th> Conv. rates </th>
1776  </tr>
1777  <tr>
1778  <td> 2 </td><td> 7.350e-02 </td><td> </td><td> 7.323e-01 </td><td> </td><td> 10.343 </td><td> </td>
1779  </tr>
1780  <tr>
1781  <td> 3 </td><td> 6.798e-03 </td><td> 3.43 </td><td> 1.716e-01 </td><td> 2.09 </td><td>4.836 </td><td> 1.09 </td>
1782  </tr>
1783  <tr>
1784  <td> 4 </td><td> 9.669e-04 </td><td> 2.81 </td><td> 6.436e-02 </td><td> 1.41 </td><td> 3.590 </td><td> 0.430 </td>
1785  </tr>
1786  <tr>
1787  <td> 5 </td><td> 1.755e-04 </td><td> 2.46 </td><td> 2.831e-02 </td><td> 1.18 </td><td>3.144 </td><td> 0.19 </td>
1788  </tr>
1789 </table>
1790 Although @f$L_2@f$ norm convergence rates of @f$u@f$ more or less
1791 follows the theoretical expectations,
1792 the @f$H^1@f$-seminorm and @f$H^2@f$-seminorm do not seem to converge as expected.
1793 Comparing results from @f$\gamma = 1@f$ and @f$\gamma = p(p+1)@f$, it is clear that
1794 @f$\gamma = p(p+1)@f$ is a better penalty.
1795 Given that @f$\gamma=1@f$ is already too small for @f$Q_2@f$ elements, it may not be surprising that if one repeated the
1796 experiment with a @f$Q_3@f$ element, the results are even more disappointing: One again only obtains convergence
1797 rates of 2, 1, zero -- i.e., no better than for the @f$Q_2@f$ element (although the errors are smaller in magnitude).
1798 Maybe surprisingly, however, one obtains more or less the expected convergence orders when using @f$Q_4@f$
1799 elements. Regardless, this uncertainty suggests that @f$\gamma=1@f$ is at best a risky choice, and at worst an
1800 unreliable one and that we should choose @f$\gamma@f$ larger.
1801 
1802 
1803 <a name="TestresultsoniQsub2subiwithigamma2i"></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>&gamma; = 2</i> </h3>
1804 
1805 
1806 Since @f$\gamma=1@f$ is clearly too small, one might conjecture that
1807 @f$\gamma=2@f$ might actually work better. Here is what one obtains in
1808 that case:
1809 
1810 <table align="center" class="doxtable">
1811  <tr>
1812  <th>Number of refinements </th><th> @f$\|u-u_h^\circ\|_{L_2}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^1}@f$ </th><th> Conv. rates </th><th> @f$|u-u_h|_{H^2}@f$ </th><th> Conv. rates </th>
1813  </tr>
1814  <tr>
1815  <td> 2 </td><td> 4.133e-02 </td><td> </td><td> 2.517e-01 </td><td> </td><td> 3.056 </td><td> </td>
1816  </tr>
1817  <tr>
1818  <td> 3 </td><td> 6.500e-03 </td><td>2.66 </td><td> 5.916e-02 </td><td> 2.08 </td><td>1.444 </td><td> 1.08 </td>
1819  </tr>
1820  <tr>
1821  <td> 4 </td><td> 6.780e-04 </td><td> 3.26 </td><td> 1.203e-02 </td><td> 2.296 </td><td> 6.151e-01 </td><td> 1.231 </td>
1822  </tr>
1823  <tr>
1824  <td> 5 </td><td> 1.622e-04 </td><td> 2.06 </td><td> 2.448e-03 </td><td> 2.297 </td><td> 2.618e-01 </td><td> 1.232 </td>
1825  </tr>
1826 </table>
1827 In this case, the convergence rates more or less follow the
1828 theoretical expectations, but, compared to the results from @f$\gamma =
1829 p(p+1)@f$, are more variable.
1830 Again, we could repeat this kind of experiment for @f$Q_3@f$ and @f$Q_4@f$ elements. In both cases, we will find that we
1831 obtain roughly the expected convergence rates. Of more interest may then be to compare the absolute
1832 size of the errors. While in the table above, for the @f$Q_2@f$ case, the errors on the finest grid are comparable between
1833 the @f$\gamma=p(p+1)@f$ and @f$\gamma=2@f$ case, for @f$Q_3@f$ the errors are substantially larger for @f$\gamma=2@f$ than for
1834 @f$\gamma=p(p+1)@f$. The same is true for the @f$Q_4@f$ case.
1835 
1836 
1837 <a name="Conclusionsforthechoiceofthepenaltyparameter"></a><h3> Conclusions for the choice of the penalty parameter </h3>
1838 
1839 
1840 The conclusions for which of the "reasonable" choices one should use for the penalty parameter
1841 is that @f$\gamma=p(p+1)@f$ yields the expected results. It is, consequently, what the code
1842 uses as currently written.
1843 
1844 
1845 <a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1846 
1847 
1848 There are a number of obvious extensions to this program that would
1849 make sense:
1850 
1851 - The program uses a square domain and a uniform mesh. Real problems
1852  don't come this way, and one should verify convergence also on
1853  domains with other shapes and, in particular, curved boundaries. One
1854  may also be interested in resolving areas of less regularity by
1855  using adaptive mesh refinement.
1856 
1857 - From a more theoretical perspective, the convergence results above
1858  only used the "broken" @f$H^2@f$ seminorm @f$|\cdot|^\circ_{H^2}@f$ instead
1859  of the "equivalent" norm @f$|\cdot|_h@f$. This is good enough to
1860  convince ourselves that the program isn't fundamentally
1861  broken. However, it might be interesting to measure the error in the
1862  actual norm for which we have theoretical results. Implementing this
1863  addition should not be overly difficult using, for example, the
1864  FEInterfaceValues class combined with MeshWorker::mesh_loop() in the
1865  same spirit as we used for the assembly of the linear system.
1866 
1867 
1868 <a name="Derivationforthesimplysupportedplates"></a> <h4> Derivation for the simply supported plates </h4>
1869 
1870 
1871  Similar to the "clamped" boundary condition addressed in the implementation,
1872  we will derive the @f$C^0@f$ IP finite element scheme for simply supported plates:
1873  @f{align*}{
1874  \Delta^2 u(\mathbf x) &= f(\mathbf x)
1875  \qquad \qquad &&\forall \mathbf x \in \Omega,
1876  u(\mathbf x) &= g(\mathbf x) \qquad \qquad
1877  &&\forall \mathbf x \in \partial\Omega, \\
1878  \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad
1879  &&\forall \mathbf x \in \partial\Omega.
1880  @f}
1881  We multiply the biharmonic equation by the test function @f$v_h@f$ and integrate over @f$ K @f$ and get:
1882  @f{align*}{
1883  \int_K v_h (\Delta^2 u_h)
1884  &= \int_K (D^2 v_h) : (D^2 u_h)
1885  + \int_{\partial K} v_h \frac{\partial (\Delta u_h)}{\partial \mathbf{n}}
1886  -\int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}).
1887  @f}
1888 
1889  Summing up over all cells @f$K \in \mathbb{T}@f$,since normal directions of @f$\Delta u_h@f$ are pointing at
1890  opposite directions on each interior edge shared by two cells and @f$v_h = 0@f$ on @f$\partial \Omega@f$,
1891  @f{align*}{
1892  \sum_{K \in \mathbb{T}} \int_{\partial K} v_h \frac{\partial (\Delta u_h)}{\partial \mathbf{n}} = 0,
1893  @f}
1894  and by the definition of jump over cell interfaces,
1895  @f{align*}{
1896  -\sum_{K \in \mathbb{T}} \int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}) = -\sum_{e \in \mathbb{F}} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} (\frac{\partial^2 u_h}{\partial \mathbf{n^2}}).
1897  @f}
1898  We separate interior faces and boundary faces of the domain,
1899  @f{align*}{
1900  -\sum_{K \in \mathbb{T}} \int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}) = -\sum_{e \in \mathbb{F}^i} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} (\frac{\partial^2 u_h}{\partial \mathbf{n^2}})
1901  - \sum_{e \in \partial \Omega} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} h,
1902  @f}
1903  where @f$\mathbb{F}^i@f$ is the set of interior faces.
1904  This leads us to
1905  @f{align*}{
1906  \sum_{K \in \mathbb{T}} \int_K (D^2 v_h) : (D^2 u_h) \ dx - \sum_{e \in \mathbb{F}^i} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} (\frac{\partial^2 u_h}{\partial \mathbf{n^2}}) \ ds
1907  = \sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx + \sum_{e\subset\partial\Omega} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} h \ ds.
1908  @f}
1909 
1910  In order to symmetrize and stabilize the discrete problem,
1911  we add symmetrization and stabilization term.
1912  We finally get the @f$C^0@f$ IP finite element scheme for the biharmonic equation:
1913  find @f$u_h@f$ such that @f$u_h =g@f$ on @f$\partial \Omega@f$ and
1914  @f{align*}{
1915  \mathcal{A}(v_h,u_h)&=\mathcal{F}(v_h) \quad \text{holds for all test functions } v_h,
1916  @f}
1917  where
1918  @f{align*}{
1919  \mathcal{A}(v_h,u_h)&=\mathcal{F}(v_h) \quad \text{holds for all test functions } v_h,
1920  @f}
1921  where
1922  @f{align*}{
1923  \mathcal{A}(v_h,u_h):=&\sum_{K \in \mathbb{T}}\int_K D^2v_h:D^2u_h \ dx
1924  \\
1925  &
1926  -\sum_{e \in \mathbb{F}^i} \int_{e}
1927  \jump{\frac{\partial v_h}{\partial \mathbf n}}
1928  \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds
1929  -\sum_{e \in \mathbb{F}^i} \int_{e}
1930  \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}}
1931  \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds
1932  \\
1933  &+ \sum_{e \in \mathbb{F}^i}
1934  \frac{\gamma}{h_e}
1935  \int_e
1936  \jump{\frac{\partial v_h}{\partial \mathbf n}}
1937  \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds,
1938  @f}
1939  and
1940  @f{align*}{
1941  \mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx
1942  +
1943  \sum_{e\subset\partial\Omega}
1944  \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} h \ ds.
1945  @f}
1946  The implementation of this boundary case is similar to "clamped" version
1947  except for `boundary_worker` is no longer needed for system assembling
1948  and the right hand side is changed according to the formulation.
1949  *
1950  *
1951 <a name="PlainProg"></a>
1952 <h1> The plain program</h1>
1953 @include "step-47.cc"
1954 */
FEInterfaceValues
Definition: fe_interface_values.h:53
MappingQ
Definition: mapping_manifold.h:33
FE_Q
Definition: fe_q.h:554
SymmetricTensor< 2, dim >
DoFTools::make_flux_sparsity_pattern
void make_flux_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern)
Definition: dof_tools_sparsity.cc:689
Triangulation< dim >
SparseMatrix< double >
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)
std::cos
inline ::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5324
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
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
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
Tensor< 1, dim >
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
SparsityPattern
Definition: sparsity_pattern.h:865
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
FESubfaceValues
Definition: fe.h:42
AffineConstraints< double >
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
Point< dim >
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
FEFaceValues
Definition: fe.h:40
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector< double >
std::sin
inline ::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5297