Reference documentation for deal.II version GIT relicensing-478-g3275795167 2024-04-24 07:10: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-61.h
Go to the documentation of this file.
1,
734 *   const unsigned int /*component*/) const
735 *   {
736 *   return 0;
737 *   }
738 *  
739 *  
740 *  
741 *   template <int dim>
742 *   class RightHandSide : public Function<dim>
743 *   {
744 *   public:
745 *   virtual double value(const Point<dim> &p,
746 *   const unsigned int component = 0) const override;
747 *   };
748 *  
749 *  
750 *  
751 *   template <int dim>
752 *   double RightHandSide<dim>::value(const Point<dim> &p,
753 *   const unsigned int /*component*/) const
754 *   {
755 *   return (2 * numbers::PI * numbers::PI * std::sin(numbers::PI * p[0]) *
756 *   std::sin(numbers::PI * p[1]));
757 *   }
758 *  
759 *  
760 *  
761 * @endcode
762 *
763 * The class that implements the exact pressure solution has an
764 * oddity in that we implement it as a vector-valued one with two
765 * components. (We say that it has two components in the constructor
766 * where we call the constructor of the base Function class.) In the
767 * `value()` function, we do not test for the value of the
768 * `component` argument, which implies that we return the same value
769 * for both components of the vector-valued function. We do this
770 * because we describe the finite element in use in this program as
771 * a vector-valued system that contains the interior and the
772 * interface pressures, and when we compute errors, we will want to
773 * use the same pressure solution to test both of these components.
774 *
775 * @code
776 *   template <int dim>
777 *   class ExactPressure : public Function<dim>
778 *   {
779 *   public:
780 *   ExactPressure()
781 *   : Function<dim>(2)
782 *   {}
783 *  
784 *   virtual double value(const Point<dim> &p,
785 *   const unsigned int component) const override;
786 *   };
787 *  
788 *  
789 *  
790 *   template <int dim>
791 *   double ExactPressure<dim>::value(const Point<dim> &p,
792 *   const unsigned int /*component*/) const
793 *   {
794 *   return std::sin(numbers::PI * p[0]) * std::sin(numbers::PI * p[1]);
795 *   }
796 *  
797 *  
798 *  
799 *   template <int dim>
800 *   class ExactVelocity : public TensorFunction<1, dim>
801 *   {
802 *   public:
803 *   ExactVelocity()
805 *   {}
806 *  
807 *   virtual Tensor<1, dim> value(const Point<dim> &p) const override;
808 *   };
809 *  
810 *  
811 *  
812 *   template <int dim>
813 *   Tensor<1, dim> ExactVelocity<dim>::value(const Point<dim> &p) const
814 *   {
815 *   Tensor<1, dim> return_value;
816 *   return_value[0] = -numbers::PI * std::cos(numbers::PI * p[0]) *
817 *   std::sin(numbers::PI * p[1]);
818 *   return_value[1] = -numbers::PI * std::sin(numbers::PI * p[0]) *
819 *   std::cos(numbers::PI * p[1]);
820 *   return return_value;
821 *   }
822 *  
823 *  
824 *  
825 * @endcode
826 *
827 *
828 * <a name="step_61-WGDarcyEquationclassimplementation"></a>
829 * <h3>WGDarcyEquation class implementation</h3>
830 *
831
832 *
833 *
834 * <a name="step_61-WGDarcyEquationWGDarcyEquation"></a>
835 * <h4>WGDarcyEquation::WGDarcyEquation</h4>
836 *
837
838 *
839 * In this constructor, we create a finite element space for vector valued
840 * functions, which will here include the ones used for the interior and
841 * interface pressures, @f$p^\circ@f$ and @f$p^\partial@f$.
842 *
843 * @code
844 *   template <int dim>
845 *   WGDarcyEquation<dim>::WGDarcyEquation(const unsigned int degree)
846 *   : fe(FE_DGQ<dim>(degree), FE_FaceQ<dim>(degree))
847 *   , dof_handler(triangulation)
848 *   , fe_dgrt(degree)
849 *   , dof_handler_dgrt(triangulation)
850 *   {}
851 *  
852 *  
853 *  
854 * @endcode
855 *
856 *
857 * <a name="step_61-WGDarcyEquationmake_grid"></a>
858 * <h4>WGDarcyEquation::make_grid</h4>
859 *
860
861 *
862 * We generate a mesh on the unit square domain and refine it.
863 *
864 * @code
865 *   template <int dim>
866 *   void WGDarcyEquation<dim>::make_grid()
867 *   {
869 *   triangulation.refine_global(5);
870 *  
871 *   std::cout << " Number of active cells: " << triangulation.n_active_cells()
872 *   << std::endl
873 *   << " Total number of cells: " << triangulation.n_cells()
874 *   << std::endl;
875 *   }
876 *  
877 *  
878 *  
879 * @endcode
880 *
881 *
882 * <a name="step_61-WGDarcyEquationsetup_system"></a>
883 * <h4>WGDarcyEquation::setup_system</h4>
884 *
885
886 *
887 * After we have created the mesh above, we distribute degrees of
888 * freedom and resize matrices and vectors. The only piece of
889 * interest in this function is how we interpolate the boundary
890 * values for the pressure. Since the pressure consists of interior
891 * and interface components, we need to make sure that we only
892 * interpolate onto that component of the vector-valued solution
893 * space that corresponds to the interface pressures (as these are
894 * the only ones that are defined on the boundary of the domain). We
895 * do this via a component mask object for only the interface
896 * pressures.
897 *
898 * @code
899 *   template <int dim>
900 *   void WGDarcyEquation<dim>::setup_system()
901 *   {
902 *   dof_handler.distribute_dofs(fe);
903 *   dof_handler_dgrt.distribute_dofs(fe_dgrt);
904 *  
905 *   std::cout << " Number of pressure degrees of freedom: "
906 *   << dof_handler.n_dofs() << std::endl;
907 *  
908 *   solution.reinit(dof_handler.n_dofs());
909 *   system_rhs.reinit(dof_handler.n_dofs());
910 *  
911 *  
912 *   {
913 *   constraints.clear();
914 *   const FEValuesExtractors::Scalar interface_pressure(1);
915 *   const ComponentMask interface_pressure_mask =
916 *   fe.component_mask(interface_pressure);
918 *   0,
919 *   BoundaryValues<dim>(),
920 *   constraints,
921 *   interface_pressure_mask);
922 *   constraints.close();
923 *   }
924 *  
925 *  
926 * @endcode
927 *
928 * In the bilinear form, there is no integration term over faces
929 * between two neighboring cells, so we can just use
930 * <code>DoFTools::make_sparsity_pattern</code> to calculate the sparse
931 * matrix.
932 *
933 * @code
934 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
935 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
936 *   sparsity_pattern.copy_from(dsp);
937 *  
938 *   system_matrix.reinit(sparsity_pattern);
939 *   }
940 *  
941 *  
942 *  
943 * @endcode
944 *
945 *
946 * <a name="step_61-WGDarcyEquationassemble_system"></a>
947 * <h4>WGDarcyEquation::assemble_system</h4>
948 *
949
950 *
951 * This function is more interesting. As detailed in the
952 * introduction, the assembly of the linear system requires us to
953 * evaluate the weak gradient of the shape functions, which is an
954 * element in the Raviart-Thomas space. As a consequence, we need to
955 * define a Raviart-Thomas finite element object, and have FEValues
956 * objects that evaluate it at quadrature points. We then need to
957 * compute the matrix @f$C^K@f$ on every cell @f$K@f$, for which we need the
958 * matrices @f$M^K@f$ and @f$G^K@f$ mentioned in the introduction.
959 *
960
961 *
962 * A point that may not be obvious is that in all previous tutorial
963 * programs, we have always called FEValues::reinit() with a cell
964 * iterator from a DoFHandler. This is so that one can call
965 * functions such as FEValuesBase::get_function_values() that
966 * extract the values of a finite element function (represented by a
967 * vector of DoF values) on the quadrature points of a cell. For
968 * this operation to work, one needs to know which vector elements
969 * correspond to the degrees of freedom on a given cell -- i.e.,
970 * exactly the kind of information and operation provided by the
971 * DoFHandler class.
972 *
973
974 *
975 * We could create a DoFHandler object for the "broken" Raviart-Thomas space
976 * (using the FE_DGRaviartThomas class), but we really don't want to here: At
977 * least in the current function, we have no need for any globally defined
978 * degrees of freedom associated with this broken space, but really only
979 * need to reference the shape functions of such a space on the current
980 * cell. As a consequence, we use the fact that one can call
981 * FEValues::reinit() also with cell iterators into Triangulation
982 * objects (rather than DoFHandler objects). In this case, FEValues
983 * can of course only provide us with information that only
984 * references information about cells, rather than degrees of freedom
985 * enumerated on these cells. So we can't use
986 * FEValuesBase::get_function_values(), but we can use
987 * FEValues::shape_value() to obtain the values of shape functions
988 * at quadrature points on the current cell. It is this kind of
989 * functionality we will make use of below. The variable that will
990 * give us this information about the Raviart-Thomas functions below
991 * is then the `fe_values_rt` (and corresponding `fe_face_values_rt`)
992 * object.
993 *
994
995 *
996 * Given this introduction, the following declarations should be
997 * pretty obvious:
998 *
999 * @code
1000 *   template <int dim>
1001 *   void WGDarcyEquation<dim>::assemble_system()
1002 *   {
1003 *   const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1004 *   const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1005 *  
1006 *   FEValues<dim> fe_values(fe,
1007 *   quadrature_formula,
1009 *   update_JxW_values);
1010 *   FEFaceValues<dim> fe_face_values(fe,
1011 *   face_quadrature_formula,
1014 *   update_JxW_values);
1015 *  
1016 *   FEValues<dim> fe_values_dgrt(fe_dgrt,
1017 *   quadrature_formula,
1020 *   update_JxW_values);
1021 *   FEFaceValues<dim> fe_face_values_dgrt(fe_dgrt,
1022 *   face_quadrature_formula,
1023 *   update_values |
1026 *   update_JxW_values);
1027 *  
1028 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1029 *   const unsigned int dofs_per_cell_dgrt = fe_dgrt.n_dofs_per_cell();
1030 *  
1031 *   const unsigned int n_q_points = fe_values.get_quadrature().size();
1032 *   const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1033 *  
1034 *   const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
1035 *  
1036 *   RightHandSide<dim> right_hand_side;
1037 *   std::vector<double> right_hand_side_values(n_q_points);
1038 *  
1039 *   const Coefficient<dim> coefficient;
1040 *   std::vector<Tensor<2, dim>> coefficient_values(n_q_points);
1041 *  
1042 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1043 *  
1044 *  
1045 * @endcode
1046 *
1047 * Next, let us declare the various cell matrices discussed in the
1048 * introduction:
1049 *
1050 * @code
1051 *   FullMatrix<double> cell_matrix_M(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
1052 *   FullMatrix<double> cell_matrix_G(dofs_per_cell_dgrt, dofs_per_cell);
1053 *   FullMatrix<double> cell_matrix_C(dofs_per_cell, dofs_per_cell_dgrt);
1055 *   Vector<double> cell_rhs(dofs_per_cell);
1056 *   Vector<double> cell_solution(dofs_per_cell);
1057 *  
1058 * @endcode
1059 *
1060 * We need <code>FEValuesExtractors</code> to access the @p interior and
1061 * @p face component of the shape functions.
1062 *
1063 * @code
1064 *   const FEValuesExtractors::Vector velocities(0);
1065 *   const FEValuesExtractors::Scalar pressure_interior(0);
1066 *   const FEValuesExtractors::Scalar pressure_face(1);
1067 *  
1068 * @endcode
1069 *
1070 * This finally gets us in position to loop over all cells. On
1071 * each cell, we will first calculate the various cell matrices
1072 * used to construct the local matrix -- as they depend on the
1073 * cell in question, they need to be re-computed on each cell. We
1074 * need shape functions for the Raviart-Thomas space as well, for
1075 * which we need to create first an iterator to the cell of the
1076 * triangulation, which we can obtain by assignment from the cell
1077 * pointing into the DoFHandler.
1078 *
1079 * @code
1080 *   for (const auto &cell : dof_handler.active_cell_iterators())
1081 *   {
1082 *   fe_values.reinit(cell);
1083 *  
1084 *   const typename Triangulation<dim>::active_cell_iterator cell_dgrt =
1085 *   cell;
1086 *   fe_values_dgrt.reinit(cell_dgrt);
1087 *  
1088 *   right_hand_side.value_list(fe_values.get_quadrature_points(),
1089 *   right_hand_side_values);
1090 *   coefficient.value_list(fe_values.get_quadrature_points(),
1091 *   coefficient_values);
1092 *  
1093 * @endcode
1094 *
1095 * The first cell matrix we will compute is the @ref GlossMassMatrix "mass matrix"
1096 * for the Raviart-Thomas space. Hence, we need to loop over
1097 * all the quadrature points for the velocity FEValues object.
1098 *
1099 * @code
1100 *   cell_matrix_M = 0;
1101 *   for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
1102 *   for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1103 *   {
1104 *   const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
1105 *   for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1106 *   {
1107 *   const Tensor<1, dim> v_k =
1108 *   fe_values_dgrt[velocities].value(k, q);
1109 *   cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
1110 *   }
1111 *   }
1112 * @endcode
1113 *
1114 * Next we take the inverse of this matrix by using
1115 * FullMatrix::gauss_jordan(). It will be used to calculate
1116 * the coefficient matrix @f$C^K@f$ later. It is worth recalling
1117 * later that `cell_matrix_M` actually contains the *inverse*
1118 * of @f$M^K@f$ after this call.
1119 *
1120 * @code
1121 *   cell_matrix_M.gauss_jordan();
1122 *  
1123 * @endcode
1124 *
1125 * From the introduction, we know that the right hand side
1126 * @f$G^K@f$ of the equation that defines @f$C^K@f$ is the difference
1127 * between a face integral and a cell integral. Here, we
1128 * approximate the negative of the contribution in the
1129 * interior. Each component of this matrix is the integral of
1130 * a product between a basis function of the polynomial space
1131 * and the divergence of a basis function of the
1132 * Raviart-Thomas space. These basis functions are defined in
1133 * the interior.
1134 *
1135 * @code
1136 *   cell_matrix_G = 0;
1137 *   for (unsigned int q = 0; q < n_q_points; ++q)
1138 *   for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1139 *   {
1140 *   const double div_v_i =
1141 *   fe_values_dgrt[velocities].divergence(i, q);
1142 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1143 *   {
1144 *   const double phi_j_interior =
1145 *   fe_values[pressure_interior].value(j, q);
1146 *  
1147 *   cell_matrix_G(i, j) -=
1148 *   (div_v_i * phi_j_interior * fe_values.JxW(q));
1149 *   }
1150 *   }
1151 *  
1152 *  
1153 * @endcode
1154 *
1155 * Next, we approximate the integral on faces by quadrature.
1156 * Each component is the integral of a product between a basis function
1157 * of the polynomial space and the dot product of a basis function of
1158 * the Raviart-Thomas space and the normal vector. So we loop over all
1159 * the faces of the element and obtain the normal vector.
1160 *
1161 * @code
1162 *   for (const auto &face : cell->face_iterators())
1163 *   {
1164 *   fe_face_values.reinit(cell, face);
1165 *   fe_face_values_dgrt.reinit(cell_dgrt, face);
1166 *  
1167 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1168 *   {
1169 *   const Tensor<1, dim> &normal = fe_face_values.normal_vector(q);
1170 *  
1171 *   for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1172 *   {
1173 *   const Tensor<1, dim> v_i =
1174 *   fe_face_values_dgrt[velocities].value(i, q);
1175 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1176 *   {
1177 *   const double phi_j_face =
1178 *   fe_face_values[pressure_face].value(j, q);
1179 *  
1180 *   cell_matrix_G(i, j) +=
1181 *   ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
1182 *   }
1183 *   }
1184 *   }
1185 *   }
1186 *  
1187 * @endcode
1188 *
1189 * @p cell_matrix_C is then the matrix product between the
1190 * transpose of @f$G^K@f$ and the inverse of the mass matrix
1191 * (where this inverse is stored in @p cell_matrix_M):
1192 *
1193 * @code
1194 *   cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
1195 *  
1196 * @endcode
1197 *
1198 * Finally we can compute the local matrix @f$A^K@f$. Element
1199 * @f$A^K_{ij}@f$ is given by @f$\int_{E} \sum_{k,l} C_{ik} C_{jl}
1200 * (\mathbf{K} \mathbf{v}_k) \cdot \mathbf{v}_l
1201 * \mathrm{d}x@f$. We have calculated the coefficients @f$C@f$ in
1202 * the previous step, and so obtain the following after
1203 * suitably re-arranging the loops:
1204 *
1205 * @code
1206 *   local_matrix = 0;
1207 *   for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
1208 *   {
1209 *   for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1210 *   {
1211 *   const Tensor<1, dim> v_k =
1212 *   fe_values_dgrt[velocities].value(k, q);
1213 *   for (unsigned int l = 0; l < dofs_per_cell_dgrt; ++l)
1214 *   {
1215 *   const Tensor<1, dim> v_l =
1216 *   fe_values_dgrt[velocities].value(l, q);
1217 *  
1218 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1219 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1220 *   local_matrix(i, j) +=
1221 *   (coefficient_values[q] * cell_matrix_C[i][k] * v_k) *
1222 *   cell_matrix_C[j][l] * v_l * fe_values_dgrt.JxW(q);
1223 *   }
1224 *   }
1225 *   }
1226 *  
1227 * @endcode
1228 *
1229 * Next, we calculate the right hand side, @f$\int_{K} f q \mathrm{d}x@f$:
1230 *
1231 * @code
1232 *   cell_rhs = 0;
1233 *   for (unsigned int q = 0; q < n_q_points; ++q)
1234 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1235 *   {
1236 *   cell_rhs(i) += (fe_values[pressure_interior].value(i, q) *
1237 *   right_hand_side_values[q] * fe_values.JxW(q));
1238 *   }
1239 *  
1240 * @endcode
1241 *
1242 * The last step is to distribute components of the local
1243 * matrix into the system matrix and transfer components of
1244 * the cell right hand side into the system right hand side:
1245 *
1246 * @code
1247 *   cell->get_dof_indices(local_dof_indices);
1248 *   constraints.distribute_local_to_global(
1249 *   local_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1250 *   }
1251 *   }
1252 *  
1253 *  
1254 *  
1255 * @endcode
1256 *
1257 *
1258 * <a name="step_61-WGDarcyEquationdimsolve"></a>
1259 * <h4>WGDarcyEquation<dim>::solve</h4>
1260 *
1261
1262 *
1263 * This step is rather trivial and the same as in many previous
1264 * tutorial programs:
1265 *
1266 * @code
1267 *   template <int dim>
1268 *   void WGDarcyEquation<dim>::solve()
1269 *   {
1270 *   SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
1271 *   SolverCG<Vector<double>> solver(solver_control);
1272 *   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
1273 *   constraints.distribute(solution);
1274 *   }
1275 *  
1276 *  
1277 * @endcode
1278 *
1279 *
1280 * <a name="step_61-WGDarcyEquationdimcompute_postprocessed_velocity"></a>
1281 * <h4>WGDarcyEquation<dim>::compute_postprocessed_velocity</h4>
1282 *
1283
1284 *
1285 * In this function, compute the velocity field from the pressure
1286 * solution previously computed. The
1287 * velocity is defined as @f$\mathbf{u}_h = \mathbf{Q}_h \left(
1288 * -\mathbf{K}\nabla_{w,d}p_h \right)@f$, which requires us to compute
1289 * many of the same terms as in the assembly of the system matrix.
1290 * There are also the matrices @f$E^K,D^K@f$ we need to assemble (see
1291 * the introduction) but they really just follow the same kind of
1292 * pattern.
1293 *
1294
1295 *
1296 * Computing the same matrices here as we have already done in the
1297 * `assemble_system()` function is of course wasteful in terms of
1298 * CPU time. Likewise, we copy some of the code from there to this
1299 * function, and this is also generally a poor idea. A better
1300 * implementation might provide for a function that encapsulates
1301 * this duplicated code. One could also think of using the classic
1302 * trade-off between computing efficiency and memory efficiency to
1303 * only compute the @f$C^K@f$ matrices once per cell during the
1304 * assembly, storing them somewhere on the side, and re-using them
1305 * here. (This is what @ref step_51 "step-51" does, for example, where the
1306 * `assemble_system()` function takes an argument that determines
1307 * whether the local matrices are recomputed, and a similar approach
1308 * -- maybe with storing local matrices elsewhere -- could be
1309 * adapted for the current program.)
1310 *
1311 * @code
1312 *   template <int dim>
1313 *   void WGDarcyEquation<dim>::compute_postprocessed_velocity()
1314 *   {
1315 *   darcy_velocity.reinit(dof_handler_dgrt.n_dofs());
1316 *  
1317 *   const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1318 *   const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1319 *  
1320 *   FEValues<dim> fe_values(fe,
1321 *   quadrature_formula,
1323 *   update_JxW_values);
1324 *  
1325 *   FEFaceValues<dim> fe_face_values(fe,
1326 *   face_quadrature_formula,
1329 *   update_JxW_values);
1330 *  
1331 *   FEValues<dim> fe_values_dgrt(fe_dgrt,
1332 *   quadrature_formula,
1335 *   update_JxW_values);
1336 *  
1337 *   FEFaceValues<dim> fe_face_values_dgrt(fe_dgrt,
1338 *   face_quadrature_formula,
1339 *   update_values |
1342 *   update_JxW_values);
1343 *  
1344 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1345 *   const unsigned int dofs_per_cell_dgrt = fe_dgrt.n_dofs_per_cell();
1346 *  
1347 *   const unsigned int n_q_points = fe_values.get_quadrature().size();
1348 *   const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1349 *  
1350 *   const unsigned int n_face_q_points = fe_face_values.get_quadrature().size();
1351 *  
1352 *  
1353 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1354 *   std::vector<types::global_dof_index> local_dof_indices_dgrt(
1355 *   dofs_per_cell_dgrt);
1356 *  
1357 *   FullMatrix<double> cell_matrix_M(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
1358 *   FullMatrix<double> cell_matrix_G(dofs_per_cell_dgrt, dofs_per_cell);
1359 *   FullMatrix<double> cell_matrix_C(dofs_per_cell, dofs_per_cell_dgrt);
1360 *   FullMatrix<double> cell_matrix_D(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
1361 *   FullMatrix<double> cell_matrix_E(dofs_per_cell_dgrt, dofs_per_cell_dgrt);
1362 *  
1363 *   Vector<double> cell_solution(dofs_per_cell);
1364 *   Vector<double> cell_velocity(dofs_per_cell_dgrt);
1365 *  
1366 *   const Coefficient<dim> coefficient;
1367 *   std::vector<Tensor<2, dim>> coefficient_values(n_q_points_dgrt);
1368 *  
1369 *   const FEValuesExtractors::Vector velocities(0);
1370 *   const FEValuesExtractors::Scalar pressure_interior(0);
1371 *   const FEValuesExtractors::Scalar pressure_face(1);
1372 *  
1373 * @endcode
1374 *
1375 * In the introduction, we explained how to calculate the numerical velocity
1376 * on the cell. We need the pressure solution values on each cell,
1377 * coefficients of the Gram matrix and coefficients of the @f$L_2@f$ projection.
1378 * We have already calculated the global solution, so we will extract the
1379 * cell solution from the global solution. The coefficients of the Gram
1380 * matrix have been calculated when we assembled the system matrix for the
1381 * pressures. We will do the same way here. For the coefficients of the
1382 * projection, we do matrix multiplication, i.e., the inverse of the Gram
1383 * matrix times the matrix with @f$(\mathbf{K} \mathbf{w}, \mathbf{w})@f$ as
1384 * components. Then, we multiply all these coefficients and call them beta.
1385 * The numerical velocity is the product of beta and the basis functions of
1386 * the Raviart-Thomas space.
1387 *
1388 * @code
1390 *   cell = dof_handler.begin_active(),
1391 *   endc = dof_handler.end(), cell_dgrt = dof_handler_dgrt.begin_active();
1392 *   for (; cell != endc; ++cell, ++cell_dgrt)
1393 *   {
1394 *   fe_values.reinit(cell);
1395 *   fe_values_dgrt.reinit(cell_dgrt);
1396 *  
1397 *   coefficient.value_list(fe_values_dgrt.get_quadrature_points(),
1398 *   coefficient_values);
1399 *  
1400 * @endcode
1401 *
1402 * The component of this <code>cell_matrix_E</code> is the integral of
1403 * @f$(\mathbf{K} \mathbf{w}, \mathbf{w})@f$. <code>cell_matrix_M</code> is
1404 * the Gram matrix.
1405 *
1406 * @code
1407 *   cell_matrix_M = 0;
1408 *   cell_matrix_E = 0;
1409 *   for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
1410 *   for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1411 *   {
1412 *   const Tensor<1, dim> v_i = fe_values_dgrt[velocities].value(i, q);
1413 *   for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1414 *   {
1415 *   const Tensor<1, dim> v_k =
1416 *   fe_values_dgrt[velocities].value(k, q);
1417 *  
1418 *   cell_matrix_E(i, k) +=
1419 *   (coefficient_values[q] * v_i * v_k * fe_values_dgrt.JxW(q));
1420 *  
1421 *   cell_matrix_M(i, k) += (v_i * v_k * fe_values_dgrt.JxW(q));
1422 *   }
1423 *   }
1424 *  
1425 * @endcode
1426 *
1427 * To compute the matrix @f$D@f$ mentioned in the introduction, we
1428 * then need to evaluate @f$D=M^{-1}E@f$ as explained in the
1429 * introduction:
1430 *
1431 * @code
1432 *   cell_matrix_M.gauss_jordan();
1433 *   cell_matrix_M.mmult(cell_matrix_D, cell_matrix_E);
1434 *  
1435 * @endcode
1436 *
1437 * Then we also need, again, to compute the matrix @f$C@f$ that is
1438 * used to evaluate the weak discrete gradient. This is the
1439 * exact same code as used in the assembly of the system
1440 * matrix, so we just copy it from there:
1441 *
1442 * @code
1443 *   cell_matrix_G = 0;
1444 *   for (unsigned int q = 0; q < n_q_points; ++q)
1445 *   for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1446 *   {
1447 *   const double div_v_i =
1448 *   fe_values_dgrt[velocities].divergence(i, q);
1449 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1450 *   {
1451 *   const double phi_j_interior =
1452 *   fe_values[pressure_interior].value(j, q);
1453 *  
1454 *   cell_matrix_G(i, j) -=
1455 *   (div_v_i * phi_j_interior * fe_values.JxW(q));
1456 *   }
1457 *   }
1458 *  
1459 *   for (const auto &face : cell->face_iterators())
1460 *   {
1461 *   fe_face_values.reinit(cell, face);
1462 *   fe_face_values_dgrt.reinit(cell_dgrt, face);
1463 *  
1464 *   for (unsigned int q = 0; q < n_face_q_points; ++q)
1465 *   {
1466 *   const Tensor<1, dim> &normal = fe_face_values.normal_vector(q);
1467 *  
1468 *   for (unsigned int i = 0; i < dofs_per_cell_dgrt; ++i)
1469 *   {
1470 *   const Tensor<1, dim> v_i =
1471 *   fe_face_values_dgrt[velocities].value(i, q);
1472 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1473 *   {
1474 *   const double phi_j_face =
1475 *   fe_face_values[pressure_face].value(j, q);
1476 *  
1477 *   cell_matrix_G(i, j) +=
1478 *   ((v_i * normal) * phi_j_face * fe_face_values.JxW(q));
1479 *   }
1480 *   }
1481 *   }
1482 *   }
1483 *   cell_matrix_G.Tmmult(cell_matrix_C, cell_matrix_M);
1484 *  
1485 * @endcode
1486 *
1487 * Finally, we need to extract the pressure unknowns that
1488 * correspond to the current cell:
1489 *
1490 * @code
1491 *   cell->get_dof_values(solution, cell_solution);
1492 *  
1493 * @endcode
1494 *
1495 * We are now in a position to compute the local velocity
1496 * unknowns (with respect to the Raviart-Thomas space we are
1497 * projecting the term @f$-\mathbf K \nabla_{w,d} p_h@f$ into):
1498 *
1499 * @code
1500 *   cell_velocity = 0;
1501 *   for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1502 *   for (unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
1503 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1504 *   cell_velocity(k) +=
1505 *   -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
1506 *  
1507 * @endcode
1508 *
1509 * We compute Darcy velocity.
1510 * This is same as cell_velocity but used to graph Darcy velocity.
1511 *
1512 * @code
1513 *   cell_dgrt->get_dof_indices(local_dof_indices_dgrt);
1514 *   for (unsigned int k = 0; k < dofs_per_cell_dgrt; ++k)
1515 *   for (unsigned int j = 0; j < dofs_per_cell_dgrt; ++j)
1516 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1517 *   darcy_velocity(local_dof_indices_dgrt[k]) +=
1518 *   -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
1519 *   }
1520 *   }
1521 *  
1522 *  
1523 *  
1524 * @endcode
1525 *
1526 *
1527 * <a name="step_61-WGDarcyEquationdimcompute_pressure_error"></a>
1528 * <h4>WGDarcyEquation<dim>::compute_pressure_error</h4>
1529 *
1530
1531 *
1532 * This part is to calculate the @f$L_2@f$ error of the pressure. We
1533 * define a vector that holds the norm of the error on each cell.
1534 * Next, we use VectorTool::integrate_difference() to compute the
1535 * error in the @f$L_2@f$ norm on each cell. However, we really only
1536 * care about the error in the interior component of the solution
1537 * vector (we can't even evaluate the interface pressures at the
1538 * quadrature points because these are all located in the interior
1539 * of cells) and consequently have to use a weight function that
1540 * ensures that the interface component of the solution variable is
1541 * ignored. This is done by using the ComponentSelectFunction whose
1542 * arguments indicate which component we want to select (component
1543 * zero, i.e., the interior pressures) and how many components there
1544 * are in total (two).
1545 *
1546 * @code
1547 *   template <int dim>
1548 *   void WGDarcyEquation<dim>::compute_pressure_error()
1549 *   {
1550 *   Vector<float> difference_per_cell(triangulation.n_active_cells());
1551 *   const ComponentSelectFunction<dim> select_interior_pressure(0, 2);
1552 *   VectorTools::integrate_difference(dof_handler,
1553 *   solution,
1554 *   ExactPressure<dim>(),
1555 *   difference_per_cell,
1556 *   QGauss<dim>(fe.degree + 2),
1558 *   &select_interior_pressure);
1559 *  
1560 *   const double L2_error = difference_per_cell.l2_norm();
1561 *   std::cout << "L2_error_pressure " << L2_error << std::endl;
1562 *   }
1563 *  
1564 *  
1565 *  
1566 * @endcode
1567 *
1568 *
1569 * <a name="step_61-WGDarcyEquationdimcompute_velocity_error"></a>
1570 * <h4>WGDarcyEquation<dim>::compute_velocity_error</h4>
1571 *
1572
1573 *
1574 * In this function, we evaluate @f$L_2@f$ errors for the velocity on
1575 * each cell, and @f$L_2@f$ errors for the flux on faces. The function
1576 * relies on the `compute_postprocessed_velocity()` function having
1577 * previous computed, which computes the velocity field based on the
1578 * pressure solution that has previously been computed.
1579 *
1580
1581 *
1582 * We are going to evaluate velocities on each cell and calculate
1583 * the difference between numerical and exact velocities.
1584 *
1585 * @code
1586 *   template <int dim>
1587 *   void WGDarcyEquation<dim>::compute_velocity_errors()
1588 *   {
1589 *   const QGauss<dim> quadrature_formula(fe_dgrt.degree + 1);
1590 *   const QGauss<dim - 1> face_quadrature_formula(fe_dgrt.degree + 1);
1591 *  
1592 *   FEValues<dim> fe_values_dgrt(fe_dgrt,
1593 *   quadrature_formula,
1596 *   update_JxW_values);
1597 *  
1598 *   FEFaceValues<dim> fe_face_values_dgrt(fe_dgrt,
1599 *   face_quadrature_formula,
1600 *   update_values |
1603 *   update_JxW_values);
1604 *  
1605 *   const unsigned int n_q_points_dgrt = fe_values_dgrt.get_quadrature().size();
1606 *   const unsigned int n_face_q_points_dgrt =
1607 *   fe_face_values_dgrt.get_quadrature().size();
1608 *  
1609 *   std::vector<Tensor<1, dim>> velocity_values(n_q_points_dgrt);
1610 *   std::vector<Tensor<1, dim>> velocity_face_values(n_face_q_points_dgrt);
1611 *  
1612 *   const FEValuesExtractors::Vector velocities(0);
1613 *  
1614 *   const ExactVelocity<dim> exact_velocity;
1615 *  
1616 *   double L2_err_velocity_cell_sqr_global = 0;
1617 *   double L2_err_flux_sqr = 0;
1618 *  
1619 * @endcode
1620 *
1621 * Having previously computed the postprocessed velocity, we here
1622 * only have to extract the corresponding values on each cell and
1623 * face and compare it to the exact values.
1624 *
1625 * @code
1626 *   for (const auto &cell_dgrt : dof_handler_dgrt.active_cell_iterators())
1627 *   {
1628 *   fe_values_dgrt.reinit(cell_dgrt);
1629 *  
1630 * @endcode
1631 *
1632 * First compute the @f$L_2@f$ error between the postprocessed velocity
1633 * field and the exact one:
1634 *
1635 * @code
1636 *   fe_values_dgrt[velocities].get_function_values(darcy_velocity,
1637 *   velocity_values);
1638 *   double L2_err_velocity_cell_sqr_local = 0;
1639 *   for (unsigned int q = 0; q < n_q_points_dgrt; ++q)
1640 *   {
1641 *   const Tensor<1, dim> velocity = velocity_values[q];
1642 *   const Tensor<1, dim> true_velocity =
1643 *   exact_velocity.value(fe_values_dgrt.quadrature_point(q));
1644 *  
1645 *   L2_err_velocity_cell_sqr_local +=
1646 *   ((velocity - true_velocity) * (velocity - true_velocity) *
1647 *   fe_values_dgrt.JxW(q));
1648 *   }
1649 *   L2_err_velocity_cell_sqr_global += L2_err_velocity_cell_sqr_local;
1650 *  
1651 * @endcode
1652 *
1653 * For reconstructing the flux we need the size of cells and
1654 * faces. Since fluxes are calculated on faces, we have the
1655 * loop over all four faces of each cell. To calculate the
1656 * face velocity, we extract values at the quadrature points from the
1657 * `darcy_velocity` which we have computed previously. Then, we
1658 * calculate the squared velocity error in normal direction. Finally, we
1659 * calculate the @f$L_2@f$ flux error on the cell by appropriately scaling
1660 * with face and cell areas and add it to the global error.
1661 *
1662 * @code
1663 *   const double cell_area = cell_dgrt->measure();
1664 *   for (const auto &face_dgrt : cell_dgrt->face_iterators())
1665 *   {
1666 *   const double face_length = face_dgrt->measure();
1667 *   fe_face_values_dgrt.reinit(cell_dgrt, face_dgrt);
1668 *   fe_face_values_dgrt[velocities].get_function_values(
1669 *   darcy_velocity, velocity_face_values);
1670 *  
1671 *   double L2_err_flux_face_sqr_local = 0;
1672 *   for (unsigned int q = 0; q < n_face_q_points_dgrt; ++q)
1673 *   {
1674 *   const Tensor<1, dim> velocity = velocity_face_values[q];
1675 *   const Tensor<1, dim> true_velocity =
1676 *   exact_velocity.value(fe_face_values_dgrt.quadrature_point(q));
1677 *  
1678 *   const Tensor<1, dim> &normal =
1679 *   fe_face_values_dgrt.normal_vector(q);
1680 *  
1681 *   L2_err_flux_face_sqr_local +=
1682 *   ((velocity * normal - true_velocity * normal) *
1683 *   (velocity * normal - true_velocity * normal) *
1684 *   fe_face_values_dgrt.JxW(q));
1685 *   }
1686 *   const double err_flux_each_face =
1687 *   L2_err_flux_face_sqr_local / face_length * cell_area;
1688 *   L2_err_flux_sqr += err_flux_each_face;
1689 *   }
1690 *   }
1691 *  
1692 * @endcode
1693 *
1694 * After adding up errors over all cells and faces, we take the
1695 * square root and get the @f$L_2@f$ errors of velocity and
1696 * flux. These we output to screen.
1697 *
1698 * @code
1699 *   const double L2_err_velocity_cell =
1700 *   std::sqrt(L2_err_velocity_cell_sqr_global);
1701 *   const double L2_err_flux_face = std::sqrt(L2_err_flux_sqr);
1702 *  
1703 *   std::cout << "L2_error_vel: " << L2_err_velocity_cell << std::endl
1704 *   << "L2_error_flux: " << L2_err_flux_face << std::endl;
1705 *   }
1706 *  
1707 *  
1708 * @endcode
1709 *
1710 *
1711 * <a name="step_61-WGDarcyEquationoutput_results"></a>
1712 * <h4>WGDarcyEquation::output_results</h4>
1713 *
1714
1715 *
1716 * We have two sets of results to output: the interior solution and
1717 * the skeleton solution. We use <code>DataOut</code> to visualize
1718 * interior results. The graphical output for the skeleton results
1719 * is done by using the DataOutFaces class.
1720 *
1721
1722 *
1723 * In both of the output files, both the interior and the face
1724 * variables are stored. For the interface output, the output file
1725 * simply contains the interpolation of the interior pressures onto
1726 * the faces, but because it is undefined which of the two interior
1727 * pressure variables you get from the two adjacent cells, it is
1728 * best to ignore the interior pressure in the interface output
1729 * file. Conversely, for the cell interior output file, it is of
1730 * course impossible to show any interface pressures @f$p^\partial@f$,
1731 * because these are only available on interfaces and not cell
1732 * interiors. Consequently, you will see them shown as an invalid
1733 * value (such as an infinity).
1734 *
1735
1736 *
1737 * For the cell interior output, we also want to output the velocity
1738 * variables. This is a bit tricky since it lives on the same mesh
1739 * but uses a different DoFHandler object (the pressure variables live
1740 * on the `dof_handler` object, the Darcy velocity on the `dof_handler_dgrt`
1741 * object). Fortunately, there are variations of the
1742 * DataOut::add_data_vector() function that allow specifying which
1743 * DoFHandler a vector corresponds to, and consequently we can visualize
1744 * the data from both DoFHandler objects within the same file.
1745 *
1746 * @code
1747 *   template <int dim>
1748 *   void WGDarcyEquation<dim>::output_results() const
1749 *   {
1750 *   {
1751 *   DataOut<dim> data_out;
1752 *  
1753 * @endcode
1754 *
1755 * First attach the pressure solution to the DataOut object:
1756 *
1757 * @code
1758 *   const std::vector<std::string> solution_names = {"interior_pressure",
1759 *   "interface_pressure"};
1760 *   data_out.add_data_vector(dof_handler, solution, solution_names);
1761 *  
1762 * @endcode
1763 *
1764 * Then do the same with the Darcy velocity field, and continue
1765 * with writing everything out into a file.
1766 *
1767 * @code
1768 *   const std::vector<std::string> velocity_names(dim, "velocity");
1769 *   const std::vector<
1771 *   velocity_component_interpretation(
1773 *   data_out.add_data_vector(dof_handler_dgrt,
1774 *   darcy_velocity,
1775 *   velocity_names,
1776 *   velocity_component_interpretation);
1777 *  
1778 *   data_out.build_patches(fe.degree);
1779 *   std::ofstream output("solution_interior.vtu");
1780 *   data_out.write_vtu(output);
1781 *   }
1782 *  
1783 *   {
1784 *   DataOutFaces<dim> data_out_faces(false);
1785 *   data_out_faces.attach_dof_handler(dof_handler);
1786 *   data_out_faces.add_data_vector(solution, "Pressure_Face");
1787 *   data_out_faces.build_patches(fe.degree);
1788 *   std::ofstream face_output("solution_interface.vtu");
1789 *   data_out_faces.write_vtu(face_output);
1790 *   }
1791 *   }
1792 *  
1793 *  
1794 * @endcode
1795 *
1796 *
1797 * <a name="step_61-WGDarcyEquationrun"></a>
1798 * <h4>WGDarcyEquation::run</h4>
1799 *
1800
1801 *
1802 * This is the final function of the main class. It calls the other functions
1803 * of our class.
1804 *
1805 * @code
1806 *   template <int dim>
1807 *   void WGDarcyEquation<dim>::run()
1808 *   {
1809 *   std::cout << "Solving problem in " << dim << " space dimensions."
1810 *   << std::endl;
1811 *   make_grid();
1812 *   setup_system();
1813 *   assemble_system();
1814 *   solve();
1815 *   compute_postprocessed_velocity();
1816 *   compute_pressure_error();
1817 *   compute_velocity_errors();
1818 *   output_results();
1819 *   }
1820 *  
1821 *   } // namespace Step61
1822 *  
1823 *  
1824 * @endcode
1825 *
1826 *
1827 * <a name="step_61-Thecodemaincodefunction"></a>
1828 * <h3>The <code>main</code> function</h3>
1829 *
1830
1831 *
1832 * This is the main function. We can change the dimension here to run in 3d.
1833 *
1834 * @code
1835 *   int main()
1836 *   {
1837 *   try
1838 *   {
1839 *   Step61::WGDarcyEquation<2> wg_darcy(0);
1840 *   wg_darcy.run();
1841 *   }
1842 *   catch (std::exception &exc)
1843 *   {
1844 *   std::cerr << std::endl
1845 *   << std::endl
1846 *   << "----------------------------------------------------"
1847 *   << std::endl;
1848 *   std::cerr << "Exception on processing: " << std::endl
1849 *   << exc.what() << std::endl
1850 *   << "Aborting!" << std::endl
1851 *   << "----------------------------------------------------"
1852 *   << std::endl;
1853 *   return 1;
1854 *   }
1855 *   catch (...)
1856 *   {
1857 *   std::cerr << std::endl
1858 *   << std::endl
1859 *   << "----------------------------------------------------"
1860 *   << std::endl;
1861 *   std::cerr << "Unknown exception!" << std::endl
1862 *   << "Aborting!" << std::endl
1863 *   << "----------------------------------------------------"
1864 *   << std::endl;
1865 *   return 1;
1866 *   }
1867 *  
1868 *   return 0;
1869 *   }
1870 * @endcode
1871<a name="step_61-Results"></a><h1>Results</h1>
1872
1873
1874We run the program with a right hand side that will produce the
1875solution @f$p = \sin(\pi x) \sin(\pi y)@f$ and with homogeneous Dirichlet
1876boundary conditions in the domain @f$\Omega = (0,1)^2@f$. In addition, we
1877choose the coefficient matrix in the differential operator
1878@f$\mathbf{K}@f$ as the identity matrix. We test this setup using
1879@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ and
1880@f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ element combinations, which one can
1881select by using the appropriate constructor argument for the
1882`WGDarcyEquation` object in `main()`. We will then visualize pressure
1883values in interiors of cells and on faces. We want to see that the
1884pressure maximum is around 1 and the minimum is around 0. With mesh
1885refinement, the convergence rates of pressure, velocity and flux
1886should then be around 1 for @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ , 2 for
1887@f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$, and 3 for
1888@f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$.
1889
1890
1891<a name="step_61-TestresultsoniWGQsub0subQsub0subRTsub0subi"></a><h3>Test results on <i>WG(Q<sub>0</sub>,Q<sub>0</sub>;RT<sub>[0]</sub>)</i></h3>
1892
1893
1894The following figures show interior pressures and face pressures using the
1895@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ element. The mesh is refined 2 times (top)
1896and 4 times (bottom), respectively. (This number can be adjusted in the
1897`make_grid()` function.) When the mesh is coarse, one can see
1898the face pressures @f$p^\partial@f$ neatly between the values of the interior
1899pressures @f$p^\circ@f$ on the two adjacent cells.
1900
1901<table align="center">
1902 <tr>
1903 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_2d_2.png" alt=""></td>
1904 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_3d_2.png" alt=""></td>
1905 </tr>
1906 <tr>
1907 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_2d_4.png" alt=""></td>
1908 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg000_3d_4.png" alt=""></td>
1909 </tr>
1910</table>
1911
1912From the figures, we can see that with the mesh refinement, the maximum and
1913minimum pressure values are approaching the values we expect.
1914Since the mesh is a rectangular mesh and numbers of cells in each direction is even, we
1915have symmetric solutions. From the 3d figures on the right,
1916we can see that on @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, the pressure is a constant
1917in the interior of the cell, as expected.
1918
1919<a name="step_61-Convergencetableforik0i"></a><h4>Convergence table for <i>k=0</i></h4>
1920
1921
1922We run the code with differently refined meshes (chosen in the `make_grid()` function)
1923and get the following convergence rates of pressure,
1924velocity, and flux (as defined in the introduction).
1925
1926<table align="center" class="doxtable">
1927 <tr>
1928 <th>number of refinements </th><th> @f$\|p-p_h^\circ\|@f$ </th><th> @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
1929 </tr>
1930 <tr>
1931 <td> 2 </td><td> 1.587e-01 </td><td> 5.113e-01 </td><td> 7.062e-01 </td>
1932 </tr>
1933 <tr>
1934 <td> 3 </td><td> 8.000e-02 </td><td> 2.529e-01 </td><td> 3.554e-01 </td>
1935 </tr>
1936 <tr>
1937 <td> 4 </td><td> 4.006e-02 </td><td> 1.260e-01 </td><td> 1.780e-01 </td>
1938 </tr>
1939 <tr>
1940 <td> 5 </td><td> 2.004e-02 </td><td> 6.297e-02 </td><td> 8.902e-02 </td>
1941 </tr>
1942 <tr>
1943 <th>Conv.rate </th><th> 1.00 </th><th> 1.00 </th><th> 1.00 </th>
1944 </tr>
1945</table>
1946
1947We can see that the convergence rates of @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$ are around 1.
1948This, of course, matches our theoretical expectations.
1949
1950
1951<a name="step_61-TestresultsoniWGQsub1subQsub1subRTsub1subi"></a><h3>Test results on <i>WG(Q<sub>1</sub>,Q<sub>1</sub>;RT<sub>[1]</sub>)</i></h3>
1952
1953
1954We can repeat the experiment from above using the next higher polynomial
1955degree:
1956The following figures are interior pressures and face pressures implemented using
1957@f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$. The mesh is refined 4 times. Compared to the
1958previous figures using
1959@f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$, on each cell, the solution is no longer constant
1960on each cell, as we now use bilinear polynomials to do the approximation.
1961Consequently, there are 4 pressure values in one interior, 2 pressure values on
1962each face.
1963
1964<table align="center">
1965 <tr>
1966 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg111_2d_4.png" alt=""></td>
1967 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg111_3d_4.png" alt=""></td>
1968 </tr>
1969</table>
1970
1971Compared to the corresponding image for the @f$\mbox{WG}(Q_0,Q_0;RT_{[0]})@f$
1972combination, the solution is now substantially more accurate and, in
1973particular so close to being continuous at the interfaces that we can
1974no longer distinguish the interface pressures @f$p^\partial@f$ from the
1975interior pressures @f$p^\circ@f$ on the adjacent cells.
1976
1977<a name="step_61-Convergencetableforik1i"></a><h4>Convergence table for <i>k=1</i></h4>
1978
1979
1980The following are the convergence rates of pressure, velocity, and flux
1981we obtain from using the @f$\mbox{WG}(Q_1,Q_1;RT_{[1]})@f$ element combination:
1982
1983<table align="center" class="doxtable">
1984 <tr>
1985 <th>number of refinements </th><th> @f$\|p-p_h^\circ\|@f$ </th><th> @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
1986 </tr>
1987 <tr>
1988 <td> 2 </td><td> 1.613e-02 </td><td> 5.093e-02 </td><td> 7.167e-02 </td>
1989 </tr>
1990 <tr>
1991 <td> 3 </td><td> 4.056e-03 </td><td> 1.276e-02 </td><td> 1.802e-02 </td>
1992 </tr>
1993 <tr>
1994 <td> 4 </td><td> 1.015e-03 </td><td> 3.191e-03 </td><td> 4.512e-03 </td>
1995 </tr>
1996 <tr>
1997 <td> 5 </td><td> 2.540e-04 </td><td> 7.979e-04 </td><td> 1.128e-03 </td>
1998 </tr>
1999 <tr>
2000 <th>Conv.rate </th><th> 2.00 </th><th> 2.00 </th><th> 2.00 </th>
2001 </tr>
2002</table>
2003
2004The convergence rates of @f$WG(Q_1,Q_1;RT_{[1]})@f$ are around 2, as expected.
2005
2006
2007
2008<a name="step_61-TestresultsoniWGQsub2subQsub2subRTsub2subi"></a><h3>Test results on <i>WG(Q<sub>2</sub>,Q<sub>2</sub>;RT<sub>[2]</sub>)</i></h3>
2009
2010
2011Let us go one polynomial degree higher.
2012The following are interior pressures and face pressures implemented using
2013@f$WG(Q_2,Q_2;RT_{[2]})@f$, with mesh size @f$h = 1/32@f$ (i.e., 5 global mesh
2014refinement steps). In the program, we use
2015`data_out_face.build_patches(fe.degree)` when generating graphical output
2016(see the documentation of DataOut::build_patches()), which here implies that
2017we divide each 2d cell interior into 4 subcells in order to provide a better
2018visualization of the quadratic polynomials.
2019<table align="center">
2020 <tr>
2021 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg222_2d_5.png" alt=""></td>
2022 <td><img src="https://www.dealii.org/images/steps/developer/step-61.wg222_3d_5.png" alt=""></td>
2023 </tr>
2024</table>
2025
2026
2027<a name="step_61-Convergencetableforik2i"></a><h4>Convergence table for <i>k=2</i></h4>
2028
2029
2030As before, we can generate convergence data for the
2031@f$L_2@f$ errors of pressure, velocity, and flux
2032using the @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ combination:
2033
2034<table align="center" class="doxtable">
2035 <tr>
2036 <th>number of refinements </th><th> @f$\|p-p_h^\circ\|@f$ </th><th> @f$\|\mathbf{u}-\mathbf{u}_h\|@f$ </th><th> @f$\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|@f$ </th>
2037 </tr>
2038 <tr>
2039 <td> 2 </td><td> 1.072e-03 </td><td> 3.375e-03 </td><td> 4.762e-03 </td>
2040 </tr>
2041 <tr>
2042 <td> 3 </td><td> 1.347e-04 </td><td> 4.233e-04 </td><td> 5.982e-04 </td>
2043 </tr>
2044 <tr>
2045 <td> 4 </td><td> 1.685e-05 </td><td> 5.295e-05 </td><td> 7.487e-05 </td>
2046 </tr>
2047 <tr>
2048 <td> 5 </td><td> 2.107e-06 </td><td> 6.620e-06 </td><td> 9.362e-06 </td>
2049 </tr>
2050 <tr>
2051 <th>Conv.rate </th><th> 3.00 </th><th> 3.00 </th><th> 3.00 </th>
2052 </tr>
2053</table>
2054
2055Once more, the convergence rates of @f$\mbox{WG}(Q_2,Q_2;RT_{[2]})@f$ is
2056as expected, with values around 3.
2057 *
2058 *
2059<a name="step_61-PlainProg"></a>
2060<h1> The plain program</h1>
2061@include "step-61.cc"
2062*/
std::vector< bool > component_mask
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
const unsigned int dofs_per_cell
const SmartPointer< const FiniteElement< dim, spacedim >, FEValuesBase< dim, spacedim > > fe
const Quadrature< dim > quadrature
Definition fe_values.h:172
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
void gauss_jordan()
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111
virtual value_type value(const Point< dim > &p) const
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > first
Definition grid_out.cc:4623
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
const std::vector< bool > & used
spacedim & mesh
Definition grid_tools.h:989
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
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 integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ReadVector< Number > &fe_function, const Function< spacedim, Number > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
spacedim const DoFHandler< dim, spacedim > const FullMatrix< double > & transfer
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)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static constexpr double PI
Definition numbers.h:259
const InputIterator OutputIterator out
Definition parallel.h:167
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation