676 *
const unsigned int = 0) const override
685 * hessian_list(
const std::vector<
Point<dim>> & points,
687 *
const unsigned int = 0) const override
689 *
for (
unsigned i = 0; i < points.size(); ++i)
691 *
const double x = points[i][0];
692 *
const double y = points[i][1];
703 *
class RightHandSide :
public Function<dim>
706 *
static_assert(dim == 2,
"Only dim==2 is implemented");
709 *
const unsigned int = 0) const override
723 * <a name=
"Themainclass"></a>
724 * <h3>The main
class</h3>
728 * The following is the principal
class of this tutorial program. It has
729 * the structure of many of the other tutorial programs and there should
730 * really be
nothing particularly surprising about its contents or
731 * the constructor that follows it.
735 *
class BiharmonicProblem
738 * BiharmonicProblem(
const unsigned int fe_degree);
744 *
void setup_system();
745 *
void assemble_system();
747 *
void compute_errors();
748 *
void output_results(
const unsigned int iteration)
const;
768 * BiharmonicProblem<dim>::BiharmonicProblem(
const unsigned int fe_degree)
779 * unit square) and
set up the constraints, vectors, and matrices on
780 * each mesh. Again, both of these are essentially unchanged from many
781 * previous tutorial programs.
785 *
void BiharmonicProblem<dim>::make_grid()
790 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
799 *
void BiharmonicProblem<dim>::setup_system()
801 * dof_handler.distribute_dofs(fe);
803 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
806 * constraints.clear();
811 * ExactSolution::Solution<dim>(),
813 * constraints.close();
818 * sparsity_pattern.copy_from(dsp);
819 * system_matrix.reinit(sparsity_pattern);
821 * solution.reinit(dof_handler.n_dofs());
822 * system_rhs.reinit(dof_handler.n_dofs());
830 * <a name=
"Assemblingthelinearsystem"></a>
831 * <h4>Assembling the linear system</h4>
835 * The following pieces of code are more interesting. They all relate to the
836 * assembly of the linear system. While assembling the cell-interior terms
837 * is not of great difficulty -- that works in essence like the assembly
838 * of the corresponding terms of the Laplace equation, and you have seen
839 * how
this works in @ref step_4
"step-4" or @ref step_6
"step-6",
for example -- the difficulty
840 * is with the penalty terms in the formulation. These require the evaluation
843 * two sides is adaptively refined, then
one actually needs an
FEFaceValues
845 * shape
functions live where, and
finally we need to ensure that every
846 * face is visited only once. All of
this is a substantial overhead to the
847 * logic we really want to implement (namely the penalty terms in the
848 * bilinear form). As a consequence, we will make use of the
851 * directly access what we really care about: jumps, averages, etc.
855 * But
this doesn
't yet solve our problem of having to keep track of
856 * which faces we have already visited when we loop over all cells and
857 * all of their faces. To make this process simpler, we use the
858 * MeshWorker::mesh_loop() function that provides a simple interface
859 * for this task: Based on the ideas outlined in the WorkStream
860 * namespace documentation, MeshWorker::mesh_loop() requires three
861 * functions that do work on cells, interior faces, and boundary
862 * faces. These functions work on scratch objects for intermediate
863 * results, and then copy the result of their computations into
864 * copy data objects from where a copier function copies them into
865 * the global matrix and right hand side objects.
869 * The following structures then provide the scratch and copy objects
870 * that are necessary for this approach. You may look up the WorkStream
871 * namespace as well as the
872 * @ref threads "Parallel computing with multiple processors"
873 * module for more information on how they typically work.
879 * ScratchData(const Mapping<dim> & mapping,
880 * const FiniteElement<dim> &fe,
881 * const unsigned int quadrature_degree,
882 * const UpdateFlags update_flags,
883 * const UpdateFlags interface_update_flags)
884 * : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
885 * , fe_interface_values(mapping,
887 * QGauss<dim - 1>(quadrature_degree),
888 * interface_update_flags)
892 * ScratchData(const ScratchData<dim> &scratch_data)
893 * : fe_values(scratch_data.fe_values.get_mapping(),
894 * scratch_data.fe_values.get_fe(),
895 * scratch_data.fe_values.get_quadrature(),
896 * scratch_data.fe_values.get_update_flags())
897 * , fe_interface_values(scratch_data.fe_values.get_mapping(),
898 * scratch_data.fe_values.get_fe(),
899 * scratch_data.fe_interface_values.get_quadrature(),
900 * scratch_data.fe_interface_values.get_update_flags())
903 * FEValues<dim> fe_values;
904 * FEInterfaceValues<dim> fe_interface_values;
911 * CopyData(const unsigned int dofs_per_cell)
912 * : cell_matrix(dofs_per_cell, dofs_per_cell)
913 * , cell_rhs(dofs_per_cell)
914 * , local_dof_indices(dofs_per_cell)
918 * CopyData(const CopyData &) = default;
921 * CopyData(CopyData &&) = default;
924 * ~CopyData() = default;
927 * CopyData &operator=(const CopyData &) = default;
930 * CopyData &operator=(CopyData &&) = default;
935 * FullMatrix<double> cell_matrix;
936 * std::vector<types::global_dof_index> joint_dof_indices;
939 * FullMatrix<double> cell_matrix;
940 * Vector<double> cell_rhs;
941 * std::vector<types::global_dof_index> local_dof_indices;
942 * std::vector<FaceData> face_data;
949 * The more interesting part is where we actually assemble the linear system.
950 * Fundamentally, this function has five parts:
951 * - The definition of the `cell_worker` lambda function, a small
952 * function that is defined within the `assemble_system()`
953 * function and that will be responsible for computing the local
954 * integrals on an individual cell. It will work on a copy of the
955 * `ScratchData` class and put its results into the corresponding
957 * - The definition of the `face_worker` lambda function that does
958 * the integration of all terms that live on the interfaces between
960 * - The definition of the `boundary_worker` function that does the
961 * same but for cell faces located on the boundary of the domain.
962 * - The definition of the `copier` function that is responsible
963 * for copying all of the data the previous three functions have
964 * put into copy objects for a single cell, into the global matrix
965 * and right hand side.
969 * The fifth part is the one where we bring all of this together.
973 * Let us go through each of these pieces necessary for the assembly
978 * void BiharmonicProblem<dim>::assemble_system()
980 * using Iterator = typename DoFHandler<dim>::active_cell_iterator;
984 * The first piece is the `cell_worker` that does the assembly
985 * on the cell interiors. It is a (lambda) function that takes
986 * a cell (input), a scratch object, and a copy object (output)
987 * as arguments. It looks like the assembly functions of many
988 * other of the tutorial programs, or at least the body of the
989 * loop over all cells.
993 * The terms we integrate here are the cell contribution
995 * A^K_{ij} = \int_K \nabla^2\varphi_i(x) : \nabla^2\varphi_j(x) dx
997 * to the global matrix, and
999 * f^K_i = \int_K \varphi_i(x) f(x) dx
1001 * to the right hand side vector.
1005 * We use the same technique as used in the assembly of @ref step_22 "step-22"
1006 * to accelerate the function: Instead of calling
1007 * `fe_values.shape_hessian(i, qpoint)` in the innermost loop,
1008 * we create a variable `hessian_i` that evaluates this
1009 * value once in the loop over `i` and re-use the so-evaluated
1010 * value in the loop over `j`. For symmetry, we do the same with a
1011 * variable `hessian_j`, although it is indeed only used once and
1012 * we could have left the call to `fe_values.shape_hessian(j,qpoint)`
1013 * in the instruction that computes the scalar product between
1017 * auto cell_worker = [&](const Iterator & cell,
1018 * ScratchData<dim> &scratch_data,
1019 * CopyData & copy_data) {
1020 * copy_data.cell_matrix = 0;
1021 * copy_data.cell_rhs = 0;
1023 * FEValues<dim> &fe_values = scratch_data.fe_values;
1024 * fe_values.reinit(cell);
1026 * cell->get_dof_indices(copy_data.local_dof_indices);
1028 * const ExactSolution::RightHandSide<dim> right_hand_side;
1030 * const unsigned int dofs_per_cell =
1031 * scratch_data.fe_values.get_fe().n_dofs_per_cell();
1033 * for (unsigned int qpoint = 0; qpoint < fe_values.n_quadrature_points;
1036 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1038 * const Tensor<2, dim> &hessian_i =
1039 * fe_values.shape_hessian(i, qpoint);
1041 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1043 * const Tensor<2, dim> &hessian_j =
1044 * fe_values.shape_hessian(j, qpoint);
1046 * copy_data.cell_matrix(i, j) +=
1047 * scalar_product(hessian_i, // nabla^2 phi_i(x)
1048 * hessian_j) * // nabla^2 phi_j(x)
1049 * fe_values.JxW(qpoint); // dx
1052 * copy_data.cell_rhs(i) +=
1053 * fe_values.shape_value(i, qpoint) * // phi_i(x)
1054 * right_hand_side.value(
1055 * fe_values.quadrature_point(qpoint)) * // f(x)
1056 * fe_values.JxW(qpoint); // dx
1064 * The next building block is the one that assembles penalty terms on each
1065 * of the interior faces of the mesh. As described in the documentation of
1066 * MeshWorker::mesh_loop(), this function receives arguments that denote
1067 * a cell and its neighboring cell, as well as (for each of the two
1068 * cells) the face (and potentially sub-face) we have to integrate
1069 * over. Again, we also get a scratch object, and a copy object
1070 * for putting the results in.
1074 * The function has three parts itself. At the top, we initialize
1075 * the FEInterfaceValues object and create a new `CopyData::FaceData`
1076 * object to store our input in. This gets pushed to the end of the
1077 * `copy_data.face_data` variable. We need to do this because
1078 * the number of faces (or subfaces) over which we integrate for a
1079 * given cell differs from cell to cell, and the sizes of these
1080 * matrices also differ, depending on what degrees of freedom
1081 * are adjacent to the face or subface. As discussed in the documentation
1082 * of MeshWorker::mesh_loop(), the copy object is reset every time a new
1083 * cell is visited, so that what we push to the end of
1084 * `copy_data.face_data()` is really all that the later `copier` function
1085 * gets to see when it copies the contributions of each cell to the global
1086 * matrix and right hand side objects.
1089 * auto face_worker = [&](const Iterator & cell,
1090 * const unsigned int &f,
1091 * const unsigned int &sf,
1092 * const Iterator & ncell,
1093 * const unsigned int &nf,
1094 * const unsigned int &nsf,
1095 * ScratchData<dim> & scratch_data,
1096 * CopyData & copy_data) {
1097 * FEInterfaceValues<dim> &fe_interface_values =
1098 * scratch_data.fe_interface_values;
1099 * fe_interface_values.reinit(cell, f, sf, ncell, nf, nsf);
1101 * copy_data.face_data.emplace_back();
1102 * CopyData::FaceData ©_data_face = copy_data.face_data.back();
1104 * copy_data_face.joint_dof_indices =
1105 * fe_interface_values.get_interface_dof_indices();
1107 * const unsigned int n_interface_dofs =
1108 * fe_interface_values.n_current_interface_dofs();
1109 * copy_data_face.cell_matrix.reinit(n_interface_dofs, n_interface_dofs);
1113 * The second part deals with determining what the penalty
1114 * parameter should be. By looking at the units of the various
1115 * terms in the bilinear form, it is clear that the penalty has
1116 * to have the form @f$\frac{\gamma}{h_K}@f$ (i.e., one over length
1117 * scale), but it is not a priori obvious how one should choose
1118 * the dimension-less number @f$\gamma@f$. From the discontinuous
1119 * Galerkin theory for the Laplace equation, one might
1120 * conjecture that the right choice is @f$\gamma=p(p+1)@f$ is the
1121 * right choice, where @f$p@f$ is the polynomial degree of the
1122 * finite element used. We will discuss this choice in a bit
1123 * more detail in the results section of this program.
1127 * In the formula above, @f$h_K@f$ is the size of cell @f$K@f$. But this
1128 * is not quite so straightforward either: If one uses highly
1129 * stretched cells, then a more involved theory says that @f$h@f$
1130 * should be replaced by the diameter of cell @f$K@f$ normal to the
1131 * direction of the edge in question. It turns out that there
1132 * is a function in deal.II for that. Secondly, @f$h_K@f$ may be
1133 * different when viewed from the two different sides of a face.
1137 * To stay on the safe side, we take the maximum of the two values.
1138 * We will note that it is possible that this computation has to be
1139 * further adjusted if one were to use hanging nodes resulting from
1140 * adaptive mesh refinement.
1143 * const unsigned int p = fe.degree;
1144 * const double gamma_over_h =
1145 * std::max((1.0 * p * (p + 1) /
1146 * cell->extent_in_direction(
1147 * GeometryInfo<dim>::unit_normal_direction[f])),
1148 * (1.0 * p * (p + 1) /
1149 * ncell->extent_in_direction(
1150 * GeometryInfo<dim>::unit_normal_direction[nf])));
1154 * Finally, and as usual, we loop over the quadrature points and
1155 * indices `i` and `j` to add up the contributions of this face
1156 * or sub-face. These are then stored in the
1157 * `copy_data.face_data` object created above. As for the cell
1158 * worker, we pull the evaluation of averages and jumps out of
1159 * the loops if possible, introducing local variables that store
1160 * these results. The assembly then only needs to use these
1161 * local variables in the innermost loop. Regarding the concrete
1162 * formula this code implements, recall that the interface terms
1163 * of the bilinear form were as follows:
1165 * -\sum_{e \in \mathbb{F}} \int_{e}
1166 * \jump{ \frac{\partial v_h}{\partial \mathbf n}}
1167 * \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds
1168 * -\sum_{e \in \mathbb{F}} \int_{e}
1169 * \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}}
1170 * \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds
1171 * + \sum_{e \in \mathbb{F}}
1172 * \frac{\gamma}{h_e}
1174 * \jump{\frac{\partial v_h}{\partial \mathbf n}}
1175 * \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds.
1179 * for (unsigned int qpoint = 0;
1180 * qpoint < fe_interface_values.n_quadrature_points;
1183 * const auto &n = fe_interface_values.normal(qpoint);
1185 * for (unsigned int i = 0; i < n_interface_dofs; ++i)
1187 * const double av_hessian_i_dot_n_dot_n =
1188 * (fe_interface_values.average_hessian(i, qpoint) * n * n);
1189 * const double jump_grad_i_dot_n =
1190 * (fe_interface_values.jump_gradient(i, qpoint) * n);
1192 * for (unsigned int j = 0; j < n_interface_dofs; ++j)
1194 * const double av_hessian_j_dot_n_dot_n =
1195 * (fe_interface_values.average_hessian(j, qpoint) * n * n);
1196 * const double jump_grad_j_dot_n =
1197 * (fe_interface_values.jump_gradient(j, qpoint) * n);
1199 * copy_data_face.cell_matrix(i, j) +=
1200 * (-av_hessian_i_dot_n_dot_n // - {grad^2 v n n }
1201 * * jump_grad_j_dot_n // [grad u n]
1202 * - av_hessian_j_dot_n_dot_n // - {grad^2 u n n }
1203 * * jump_grad_i_dot_n // [grad v n]
1205 * gamma_over_h * // gamma/h
1206 * jump_grad_i_dot_n * // [grad v n]
1207 * jump_grad_j_dot_n) * // [grad u n]
1208 * fe_interface_values.JxW(qpoint); // dx
1217 * The third piece is to do the same kind of assembly for faces that
1218 * are at the boundary. The idea is the same as above, of course,
1219 * with only the difference that there are now penalty terms that
1220 * also go into the right hand side.
1224 * As before, the first part of the function simply sets up some
1228 * auto boundary_worker = [&](const Iterator & cell,
1229 * const unsigned int &face_no,
1230 * ScratchData<dim> & scratch_data,
1231 * CopyData & copy_data) {
1232 * FEInterfaceValues<dim> &fe_interface_values =
1233 * scratch_data.fe_interface_values;
1234 * fe_interface_values.reinit(cell, face_no);
1235 * const auto &q_points = fe_interface_values.get_quadrature_points();
1237 * copy_data.face_data.emplace_back();
1238 * CopyData::FaceData ©_data_face = copy_data.face_data.back();
1240 * const unsigned int n_dofs =
1241 * fe_interface_values.n_current_interface_dofs();
1242 * copy_data_face.joint_dof_indices =
1243 * fe_interface_values.get_interface_dof_indices();
1245 * copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
1247 * const std::vector<double> &JxW = fe_interface_values.get_JxW_values();
1248 * const std::vector<Tensor<1, dim>> &normals =
1249 * fe_interface_values.get_normal_vectors();
1252 * const ExactSolution::Solution<dim> exact_solution;
1253 * std::vector<Tensor<1, dim>> exact_gradients(q_points.size());
1254 * exact_solution.gradient_list(q_points, exact_gradients);
1259 * Positively, because we now only deal with one cell adjacent to the
1260 * face (as we are on the boundary), the computation of the penalty
1261 * factor @f$\gamma@f$ is substantially simpler:
1264 * const unsigned int p = fe.degree;
1265 * const double gamma_over_h =
1266 * (1.0 * p * (p + 1) /
1267 * cell->extent_in_direction(
1268 * GeometryInfo<dim>::unit_normal_direction[face_no]));
1272 * The third piece is the assembly of terms. This is now
1273 * slightly more involved since these contains both terms for
1274 * the matrix and for the right hand side. The former is exactly
1275 * the same as for the interior faces stated above if one just
1276 * defines the jump and average appropriately (which is what the
1277 * FEInterfaceValues class does). The latter requires us to
1278 * evaluate the boundary conditions @f$j(\mathbf x)@f$, which in the
1279 * current case (where we know the exact solution) we compute
1280 * from @f$j(\mathbf x) = \frac{\partial u(\mathbf x)}{\partial
1281 * {\mathbf n}}@f$. The term to be added to the right hand side
1283 * @f$\frac{\gamma}{h_e}\int_e
1284 * \jump{\frac{\partial v_h}{\partial \mathbf n}} j \ ds@f$.
1287 * for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
1289 * const auto &n = normals[qpoint];
1291 * for (unsigned int i = 0; i < n_dofs; ++i)
1293 * const double av_hessian_i_dot_n_dot_n =
1294 * (fe_interface_values.average_hessian(i, qpoint) * n * n);
1295 * const double jump_grad_i_dot_n =
1296 * (fe_interface_values.jump_gradient(i, qpoint) * n);
1298 * for (unsigned int j = 0; j < n_dofs; ++j)
1300 * const double av_hessian_j_dot_n_dot_n =
1301 * (fe_interface_values.average_hessian(j, qpoint) * n * n);
1302 * const double jump_grad_j_dot_n =
1303 * (fe_interface_values.jump_gradient(j, qpoint) * n);
1305 * copy_data_face.cell_matrix(i, j) +=
1306 * (-av_hessian_i_dot_n_dot_n // - {grad^2 v n n}
1307 * * jump_grad_j_dot_n // [grad u n]
1309 * - av_hessian_j_dot_n_dot_n // - {grad^2 u n n}
1310 * * jump_grad_i_dot_n // [grad v n]
1312 * + gamma_over_h // gamma/h
1313 * * jump_grad_i_dot_n // [grad v n]
1314 * * jump_grad_j_dot_n // [grad u n]
1316 * JxW[qpoint]; // dx
1319 * copy_data.cell_rhs(i) +=
1320 * (-av_hessian_i_dot_n_dot_n * // - {grad^2 v n n }
1321 * (exact_gradients[qpoint] * n) // (grad u_exact . n)
1323 * gamma_over_h // gamma/h
1324 * * jump_grad_i_dot_n // [grad v n]
1325 * * (exact_gradients[qpoint] * n) // (grad u_exact . n)
1327 * JxW[qpoint]; // dx
1334 * Part 4 is a small function that copies the data produced by the
1335 * cell, interior, and boundary face assemblers above into the
1336 * global matrix and right hand side vector. There really is not
1337 * very much to do here: We distribute the cell matrix and right
1338 * hand side contributions as we have done in almost all of the
1339 * other tutorial programs using the constraints objects. We then
1340 * also have to do the same for the face matrix contributions
1341 * that have gained content for the faces (interior and boundary)
1342 * and that the `face_worker` and `boundary_worker` have added
1343 * to the `copy_data.face_data` array.
1346 * auto copier = [&](const CopyData ©_data) {
1347 * constraints.distribute_local_to_global(copy_data.cell_matrix,
1348 * copy_data.cell_rhs,
1349 * copy_data.local_dof_indices,
1353 * for (auto &cdf : copy_data.face_data)
1355 * constraints.distribute_local_to_global(cdf.cell_matrix,
1356 * cdf.joint_dof_indices,
1364 * Having set all of this up, what remains is to just create a scratch
1365 * and copy data object and call the MeshWorker::mesh_loop() function
1366 * that then goes over all cells and faces, calls the respective workers
1367 * on them, and then the copier function that puts things into the
1368 * global matrix and right hand side. As an additional benefit,
1369 * MeshWorker::mesh_loop() does all of this in parallel, using
1370 * as many processor cores as your machine happens to have.
1373 * const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
1374 * ScratchData<dim> scratch_data(mapping,
1377 * update_values | update_gradients |
1378 * update_hessians | update_quadrature_points |
1379 * update_JxW_values,
1380 * update_values | update_gradients |
1381 * update_hessians | update_quadrature_points |
1382 * update_JxW_values | update_normal_vectors);
1383 * CopyData copy_data(dof_handler.get_fe().n_dofs_per_cell());
1384 * MeshWorker::mesh_loop(dof_handler.begin_active(),
1385 * dof_handler.end(),
1390 * MeshWorker::assemble_own_cells |
1391 * MeshWorker::assemble_boundary_faces |
1392 * MeshWorker::assemble_own_interior_faces_once,
1402 * <a name="Solvingthelinearsystemandpostprocessing"></a>
1403 * <h4>Solving the linear system and postprocessing</h4>
1407 * The show is essentially over at this point: The remaining functions are
1408 * not overly interesting or novel. The first one simply uses a direct
1409 * solver to solve the linear system (see also @ref step_29 "step-29"):
1412 * template <int dim>
1413 * void BiharmonicProblem<dim>::solve()
1415 * std::cout << " Solving system..." << std::endl;
1417 * SparseDirectUMFPACK A_direct;
1418 * A_direct.initialize(system_matrix);
1419 * A_direct.vmult(solution, system_rhs);
1421 * constraints.distribute(solution);
1428 * The next function evaluates the error between the computed solution
1429 * and the exact solution (which is known here because we have chosen
1430 * the right hand side and boundary values in a way so that we know
1431 * the corresponding solution). In the first two code blocks below,
1432 * we compute the error in the @f$L_2@f$ norm and the @f$H^1@f$ semi-norm.
1435 * template <int dim>
1436 * void BiharmonicProblem<dim>::compute_errors()
1439 * Vector<float> norm_per_cell(triangulation.n_active_cells());
1440 * VectorTools::integrate_difference(mapping,
1443 * ExactSolution::Solution<dim>(),
1445 * QGauss<dim>(fe.degree + 2),
1446 * VectorTools::L2_norm);
1447 * const double error_norm =
1448 * VectorTools::compute_global_error(triangulation,
1450 * VectorTools::L2_norm);
1451 * std::cout << " Error in the L2 norm : " << error_norm
1456 * Vector<float> norm_per_cell(triangulation.n_active_cells());
1457 * VectorTools::integrate_difference(mapping,
1460 * ExactSolution::Solution<dim>(),
1462 * QGauss<dim>(fe.degree + 2),
1463 * VectorTools::H1_seminorm);
1464 * const double error_norm =
1465 * VectorTools::compute_global_error(triangulation,
1467 * VectorTools::H1_seminorm);
1468 * std::cout << " Error in the H1 seminorm : " << error_norm
1474 * Now also compute an approximation to the @f$H^2@f$ seminorm error. The actual
1475 * @f$H^2@f$ seminorm would require us to integrate second derivatives of the
1476 * solution @f$u_h@f$, but given the Lagrange shape functions we use, @f$u_h@f$ of
1477 * course has kinks at the interfaces between cells, and consequently second
1478 * derivatives are singular at interfaces. As a consequence, we really only
1479 * integrate over the interior of cells and ignore the interface
1480 * contributions. This is *not* an equivalent norm to the energy norm for
1481 * the problem, but still gives us an idea of how fast the error converges.
1485 * We note that one could address this issue by defining a norm that
1486 * is equivalent to the energy norm. This would involve adding up not
1487 * only the integrals over cell interiors as we do below, but also adding
1488 * penalty terms for the jump of the derivative of @f$u_h@f$ across interfaces,
1489 * with an appropriate scaling of the two kinds of terms. We will leave
1490 * this for later work.
1494 * const QGauss<dim> quadrature_formula(fe.degree + 2);
1495 * ExactSolution::Solution<dim> exact_solution;
1496 * Vector<double> error_per_cell(triangulation.n_active_cells());
1498 * FEValues<dim> fe_values(mapping,
1500 * quadrature_formula,
1501 * update_values | update_hessians |
1502 * update_quadrature_points | update_JxW_values);
1504 * FEValuesExtractors::Scalar scalar(0);
1505 * const unsigned int n_q_points = quadrature_formula.size();
1507 * std::vector<SymmetricTensor<2, dim>> exact_hessians(n_q_points);
1508 * std::vector<Tensor<2, dim>> hessians(n_q_points);
1509 * for (auto &cell : dof_handler.active_cell_iterators())
1511 * fe_values.reinit(cell);
1512 * fe_values[scalar].get_function_hessians(solution, hessians);
1513 * exact_solution.hessian_list(fe_values.get_quadrature_points(),
1516 * double local_error = 0;
1517 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1520 * ((exact_hessians[q_point] - hessians[q_point]).norm_square() *
1521 * fe_values.JxW(q_point));
1523 * error_per_cell[cell->active_cell_index()] = std::sqrt(local_error);
1526 * const double error_norm = error_per_cell.l2_norm();
1527 * std::cout << " Error in the broken H2 seminorm: " << error_norm
1536 * Equally uninteresting is the function that generates graphical output.
1537 * It looks exactly like the one in @ref step_6 "step-6", for example.
1540 * template <int dim>
1542 * BiharmonicProblem<dim>::output_results(const unsigned int iteration) const
1544 * std::cout << " Writing graphical output..." << std::endl;
1546 * DataOut<dim> data_out;
1548 * data_out.attach_dof_handler(dof_handler);
1549 * data_out.add_data_vector(solution, "solution");
1550 * data_out.build_patches();
1552 * const std::string filename =
1553 * ("output_" + Utilities::int_to_string(iteration, 6) + ".vtu");
1554 * std::ofstream output_vtu(filename);
1555 * data_out.write_vtu(output_vtu);
1562 * The same is true for the `run()` function: Just like in previous
1566 * template <int dim>
1567 * void BiharmonicProblem<dim>::run()
1571 * const unsigned int n_cycles = 4;
1572 * for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1574 * std::cout << "Cycle " << cycle << " of " << n_cycles << std::endl;
1576 * triangulation.refine_global(1);
1579 * assemble_system();
1582 * output_results(cycle);
1585 * std::cout << std::endl;
1588 * } // namespace Step47
1595 * <a name="Themainfunction"></a>
1596 * <h3>The main() function</h3>
1600 * Finally for the `main()` function. There is, again, not very much to see
1601 * here: It looks like the ones in previous tutorial programs. There
1602 * is a variable that allows selecting the polynomial degree of the element
1603 * we want to use for solving the equation. Because the C0IP formulation
1604 * we use requires the element degree to be at least two, we check with
1605 * an assertion that whatever one sets for the polynomial degree actually
1613 * using namespace dealii;
1614 * using namespace Step47;
1616 * const unsigned int fe_degree = 2;
1617 * Assert(fe_degree >= 2,
1618 * ExcMessage("The C0IP formulation for the biharmonic problem "
1619 * "only works if one uses elements of polynomial "
1620 * "degree at least 2."));
1622 * BiharmonicProblem<2> biharmonic_problem(fe_degree);
1623 * biharmonic_problem.run();
1625 * catch (std::exception &exc)
1627 * std::cerr << std::endl
1629 * << "----------------------------------------------------"
1631 * std::cerr << "Exception on processing: " << std::endl
1632 * << exc.what() << std::endl
1633 * << "Aborting!" << std::endl
1634 * << "----------------------------------------------------"
1641 * std::cerr << std::endl
1643 * << "----------------------------------------------------"
1645 * std::cerr << "Unknown exception!" << std::endl
1646 * << "Aborting!" << std::endl
1647 * << "----------------------------------------------------"
1655<a name="Results"></a><h1>Results</h1>
1658We run the program with right hand side and boundary values as
1659discussed in the introduction. These will produce the
1660solution @f$u = \sin(\pi x) \sin(\pi y)@f$ on the domain @f$\Omega = (0,1)^2@f$.
1661We test this setup using @f$Q_2@f$, @f$Q_3@f$, and @f$Q_4@f$ elements, which one can
1662change via the `fe_degree` variable in the `main()` function. With mesh
1663refinement, the @f$L_2@f$ convergence rates, @f$H^1@f$-seminorm rate,
1664and @f$H^2@f$-seminorm convergence of @f$u@f$
1665should then be around 2, 2, 1 for @f$Q_2@f$ (with the @f$L_2@f$ norm
1666sub-optimal as discussed in the introduction); 4, 3, 2 for
1667@f$Q_3@f$; and 5, 4, 3 for @f$Q_4@f$, respectively.
1669From the literature, it is not immediately clear what
1670the penalty parameter @f$\gamma@f$ should be. For example,
1671@cite Brenner2009 state that it needs to be larger than one, and
1672choose @f$\gamma=5@f$. The FEniCS/Dolphin tutorial chooses it as
1674https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html
1675. @cite Wells2007 uses a value for @f$\gamma@f$ larger than the
1676number of edges belonging to an element for Kirchhoff plates (see
1677their Section 4.2). This suggests that maybe
1678@f$\gamma = 1@f$, @f$2@f$, are too small; on the other hand, a value
1679@f$p(p+1)@f$ would be reasonable,
1680where @f$p@f$ is the degree of polynomials. The last of these choices is
1681the one one would expect to work by comparing
1682to the discontinuous Galerkin formulations for the Laplace equation
1683(see, for example, the discussions in @ref step_39 "step-39" and @ref step_74 "step-74"),
1684and it will turn out to also work here.
1685But we should check what value of @f$\gamma@f$ is right, and we will do so
1686below; changing @f$\gamma@f$ is easy in the two `face_worker` and
1687`boundary_worker` functions defined in `assemble_system()`.
1690<a name="TestresultsoniQsub2subiwithigammapp1i"></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>γ = p(p+1)</i> </h3>
1693We run the code with differently refined meshes
1694and get the following convergence rates.
1696<table align="center" class="doxtable">
1698 <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>
1701 <td> 2 </td><td> 8.780e-03 </td><td> </td><td> 7.095e-02 </td><td> </td><td> 1.645 </td><td> </td>
1704 <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>
1707 <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>
1710 <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>
1713We can see that the @f$L_2@f$ convergence rates are around 2,
1714@f$H^1@f$-seminorm convergence rates are around 2,
1715and @f$H^2@f$-seminorm convergence rates are around 1. The latter two
1716match the theoretically expected rates; for the former, we have no
1717theorem but are not surprised that it is sub-optimal given the remark
1721<a name="TestresultsoniQsub3subiwithigammapp1i"></a><h3>Test results on <i>Q<sub>3</sub></i> with <i>γ = p(p+1)</i> </h3>
1725<table align="center" class="doxtable">
1727 <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>
1730 <td> 2 </td><td> 2.045e-04 </td><td> </td><td> 4.402e-03 </td><td> </td><td> 1.641e-01 </td><td> </td>
1733 <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>
1736 <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>
1739 <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>
1742We can see that the @f$L_2@f$ convergence rates are around 4,
1743@f$H^1@f$-seminorm convergence rates are around 3,
1744and @f$H^2@f$-seminorm convergence rates are around 2.
1745This, of course, matches our theoretical expectations.
1748<a name="TestresultsoniQsub4subiwithigammapp1i"></a><h3>Test results on <i>Q<sub>4</sub></i> with <i>γ = p(p+1)</i> </h3>
1751<table align="center" class="doxtable">
1753 <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>
1756 <td> 2 </td><td> 6.510e-06 </td><td> </td><td> 2.215e-04 </td><td> </td><td> 1.275e-02 </td><td> </td>
1759 <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>
1762 <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>
1765 <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>
1768We can see that the @f$L_2@f$ norm convergence rates are around 5,
1769@f$H^1@f$-seminorm convergence rates are around 4,
1770and @f$H^2@f$-seminorm convergence rates are around 3.
1771On the finest mesh, the @f$L_2@f$ norm convergence rate
1772is much smaller than our theoretical expectations
1773because the linear solver becomes the limiting factor due
1774to round-off. Of course the @f$L_2@f$ error is also very small already in
1778<a name="TestresultsoniQsub2subiwithigamma1i"></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>γ = 1</i> </h3>
1781For comparison with the results above, let us now also consider the
1782case where we simply choose @f$\gamma=1@f$:
1784<table align="center" class="doxtable">
1786 <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>
1789 <td> 2 </td><td> 7.350e-02 </td><td> </td><td> 7.323e-01 </td><td> </td><td> 10.343 </td><td> </td>
1792 <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>
1795 <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>
1798 <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>
1801Although @f$L_2@f$ norm convergence rates of @f$u@f$ more or less
1802follows the theoretical expectations,
1803the @f$H^1@f$-seminorm and @f$H^2@f$-seminorm do not seem to converge as expected.
1804Comparing results from @f$\gamma = 1@f$ and @f$\gamma = p(p+1)@f$, it is clear that
1805@f$\gamma = p(p+1)@f$ is a better penalty.
1806Given 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
1807experiment with a @f$Q_3@f$ element, the results are even more disappointing: One again only obtains convergence
1808rates of 2, 1, zero -- i.e., no better than for the @f$Q_2@f$ element (although the errors are smaller in magnitude).
1809Maybe surprisingly, however, one obtains more or less the expected convergence orders when using @f$Q_4@f$
1810elements. Regardless, this uncertainty suggests that @f$\gamma=1@f$ is at best a risky choice, and at worst an
1811unreliable one and that we should choose @f$\gamma@f$ larger.
1814<a name="TestresultsoniQsub2subiwithigamma2i"></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>γ = 2</i> </h3>
1817Since @f$\gamma=1@f$ is clearly too small, one might conjecture that
1818@f$\gamma=2@f$ might actually work better. Here is what one obtains in
1821<table align="center" class="doxtable">
1823 <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>
1826 <td> 2 </td><td> 4.133e-02 </td><td> </td><td> 2.517e-01 </td><td> </td><td> 3.056 </td><td> </td>
1829 <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>
1832 <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>
1835 <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>
1838In this case, the convergence rates more or less follow the
1839theoretical expectations, but, compared to the results from @f$\gamma =
1840p(p+1)@f$, are more variable.
1841Again, 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
1842obtain roughly the expected convergence rates. Of more interest may then be to compare the absolute
1843size of the errors. While in the table above, for the @f$Q_2@f$ case, the errors on the finest grid are comparable between
1844the @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
1845@f$\gamma=p(p+1)@f$. The same is true for the @f$Q_4@f$ case.
1848<a name="Conclusionsforthechoiceofthepenaltyparameter"></a><h3> Conclusions for the choice of the penalty parameter </h3>
1851The conclusions for which of the "reasonable" choices one should use for the penalty parameter
1852is that @f$\gamma=p(p+1)@f$ yields the expected results. It is, consequently, what the code
1853uses as currently written.
1856<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1859There are a number of obvious extensions to this program that would
1862- The program uses a square domain and a uniform mesh. Real problems
1863 don't come
this way, and
one should verify convergence also on
1864 domains with other shapes and, in particular, curved boundaries. One
1865 may also be interested in resolving areas of less regularity by
1866 using adaptive mesh refinement.
1868- From a more theoretical perspective, the convergence results above
1869 only used the
"broken" @f$H^2@f$ seminorm @f$|\cdot|^\circ_{H^2}@f$ instead
1870 of the
"equivalent" norm @f$|\cdot|_h@f$. This is good enough to
1871 convince ourselves that the program isn
't fundamentally
1872 broken. However, it might be interesting to measure the error in the
1873 actual norm for which we have theoretical results. Implementing this
1874 addition should not be overly difficult using, for example, the
1875 FEInterfaceValues class combined with MeshWorker::mesh_loop() in the
1876 same spirit as we used for the assembly of the linear system.
1879<a name="Derivationforthesimplysupportedplates"></a> <h4> Derivation for the simply supported plates </h4>
1882 Similar to the "clamped" boundary condition addressed in the implementation,
1883 we will derive the @f$C^0@f$ IP finite element scheme for simply supported plates:
1885 \Delta^2 u(\mathbf x) &= f(\mathbf x)
1886 \qquad \qquad &&\forall \mathbf x \in \Omega,
1887 u(\mathbf x) &= g(\mathbf x) \qquad \qquad
1888 &&\forall \mathbf x \in \partial\Omega, \\
1889 \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad
1890 &&\forall \mathbf x \in \partial\Omega.
1892 We multiply the biharmonic equation by the test function @f$v_h@f$ and integrate over @f$ K @f$ and get:
1894 \int_K v_h (\Delta^2 u_h)
1895 &= \int_K (D^2 v_h) : (D^2 u_h)
1896 + \int_{\partial K} v_h \frac{\partial (\Delta u_h)}{\partial \mathbf{n}}
1897 -\int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}).
1900 Summing up over all cells @f$K \in \mathbb{T}@f$,since normal directions of @f$\Delta u_h@f$ are pointing at
1901 opposite directions on each interior edge shared by two cells and @f$v_h = 0@f$ on @f$\partial \Omega@f$,
1903 \sum_{K \in \mathbb{T}} \int_{\partial K} v_h \frac{\partial (\Delta u_h)}{\partial \mathbf{n}} = 0,
1905 and by the definition of jump over cell interfaces,
1907 -\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}}).
1909 We separate interior faces and boundary faces of the domain,
1911 -\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}})
1912 - \sum_{e \in \partial \Omega} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} h,
1914 where @f$\mathbb{F}^i@f$ is the set of interior faces.
1917 \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
1918 = \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.
1921 In order to symmetrize and stabilize the discrete problem,
1922 we add symmetrization and stabilization term.
1923 We finally get the @f$C^0@f$ IP finite element scheme for the biharmonic equation:
1924 find @f$u_h@f$ such that @f$u_h =g@f$ on @f$\partial \Omega@f$ and
1926 \mathcal{A}(v_h,u_h)&=\mathcal{F}(v_h) \quad \text{holds for all test functions } v_h,
1930 \mathcal{A}(v_h,u_h):=&\sum_{K \in \mathbb{T}}\int_K D^2v_h:D^2u_h \ dx
1933 -\sum_{e \in \mathbb{F}^i} \int_{e}
1934 \jump{\frac{\partial v_h}{\partial \mathbf n}}
1935 \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds
1936 -\sum_{e \in \mathbb{F}^i} \int_{e}
1937 \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}}
1938 \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds
1940 &+ \sum_{e \in \mathbb{F}^i}
1943 \jump{\frac{\partial v_h}{\partial \mathbf n}}
1944 \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds,
1948 \mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx
1950 \sum_{e\subset\partial\Omega}
1951 \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} h \ ds.
1953 The implementation of this boundary case is similar to the "clamped" version
1954 except that `boundary_worker` is no longer needed for system assembling
1955 and the right hand side is changed according to the formulation.
1958<a name="PlainProg"></a>
1959<h1> The plain program</h1>
1960@include "step-47.cc"
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
__global__ void set(Number *val, const Number s, const size_type N)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
static const types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
static constexpr double PI
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation