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