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