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