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