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-3.h
Go to the documentation of this file.
1 1)
767 *   , dof_handler(triangulation)
768 *   {}
769 *  
770 *  
771 * @endcode
772 *
773 *
774 * <a name="step_3-Step3make_grid"></a>
776 *
777
778 *
779 * Now, the first thing we've got to do is to generate the triangulation on
780 * which we would like to do our computation and number each vertex with a
781 * degree of freedom. We have seen these two steps in @ref step_1 "step-1" and @ref step_2 "step-2"
782 * before, respectively.
783 *
784
785 *
786 * This function does the first part, creating the mesh. We create the grid
787 * and refine all cells five times. Since the initial grid (which is the
788 * square @f$[-1,1] \times [-1,1]@f$) consists of only one cell, the final grid
789 * has 32 times 32 cells, for a total of 1024.
790 *
791
792 *
793 * Unsure that 1024 is the correct number? We can check that by outputting the
794 * number of cells using the <code>n_active_cells()</code> function on the
795 * triangulation.
796 *
797 * @code
798 *   void Step3::make_grid()
799 *   {
800 *   GridGenerator::hyper_cube(triangulation, -1, 1);
801 *   triangulation.refine_global(5);
802 *  
803 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
804 *   << std::endl;
805 *   }
806 *  
807 * @endcode
808 *
809 * @note We call the Triangulation::n_active_cells() function, rather than
810 * Triangulation::n_cells(). Here, <i>active</i> means the cells that aren't
812 * cells, namely the parent cells of the finest cells, their parents, etc, up
813 * to the one cell which made up the initial grid. Of course, on the next
814 * coarser level, the number of cells is one quarter that of the cells on the
815 * finest level, i.e. 256, then 64, 16, 4, and 1. If you called
816 * <code>triangulation.n_cells()</code> instead in the code above, you would
817 * consequently get a value of 1365 instead. On the other hand, the number of
818 * cells (as opposed to the number of active cells) is not typically of much
819 * interest, so there is no good reason to print it.
820 *
821
822 *
823 *
824
825 *
826 *
827 * <a name="step_3-Step3setup_system"></a>
829 *
830
831 *
832 * Next we enumerate all the degrees of freedom and set up matrix and vector
833 * objects to hold the system data. Enumerating is done by using
834 * DoFHandler::distribute_dofs(), as we have seen in the @ref step_2 "step-2" example. Since
835 * we use the FE_Q class and have set the polynomial degree to 1 in the
836 * constructor, i.e. bilinear elements, this associates one degree of freedom
837 * with each vertex. While we're at generating output, let us also take a look
838 * at how many degrees of freedom are generated:
839 *
840 * @code
841 *   void Step3::setup_system()
842 *   {
843 *   dof_handler.distribute_dofs(fe);
844 *   std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
845 *   << std::endl;
846 * @endcode
847 *
848 * There should be one DoF for each vertex. Since we have a 32 times 32
849 * grid, the number of DoFs should be 33 times 33, or 1089.
850 *
851
852 *
853 * As we have seen in the previous example, we set up a sparsity pattern by
854 * first creating a temporary structure, tagging those entries that might be
855 * nonzero, and then copying the data over to the SparsityPattern object
856 * that can then be used by the system matrix.
857 *
858 * @code
859 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
860 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
861 *   sparsity_pattern.copy_from(dsp);
862 *  
863 * @endcode
864 *
865 * Note that the SparsityPattern object does not hold the values of the
866 * matrix, it only stores the places where entries are. The entries
867 * themselves are stored in objects of type SparseMatrix, of which our
868 * variable system_matrix is one.
869 *
870
871 *
872 * The distinction between sparsity pattern and matrix was made to allow
873 * several matrices to use the same sparsity pattern. This may not seem
874 * relevant here, but when you consider the size which matrices can have,
875 * and that it may take some time to build the sparsity pattern, this
876 * becomes important in large-scale problems if you have to store several
877 * matrices in your program.
878 *
879 * @code
880 *   system_matrix.reinit(sparsity_pattern);
881 *  
882 * @endcode
883 *
884 * The last thing to do in this function is to set the sizes of the right
885 * hand side vector and the solution vector to the right values:
886 *
887 * @code
888 *   solution.reinit(dof_handler.n_dofs());
889 *   system_rhs.reinit(dof_handler.n_dofs());
890 *   }
891 *  
892 * @endcode
893 *
894 *
895 * <a name="step_3-Step3assemble_system"></a>
896 * <h4>Step3::assemble_system</h4>
897 *
898
899 *
900 *
901
902 *
903 * The next step is to compute the entries of the matrix and right hand side
904 * that form the linear system from which we compute the solution. This is the
905 * central function of each finite element program and we have discussed the
906 * primary steps in the introduction already.
907 *
908
909 *
910 * The general approach to assemble matrices and vectors is to loop over all
911 * cells, and on each cell compute the contribution of that cell to the global
912 * matrix and right hand side by quadrature. The point to realize now is that
913 * we need the values of the shape functions at the locations of quadrature
914 * points on the real cell. However, both the finite element shape functions
915 * as well as the quadrature points are only defined on the reference
916 * cell. They are therefore of little help to us, and we will in fact hardly
917 * ever query information about finite element shape functions or quadrature
918 * points from these objects directly.
919 *
920
921 *
922 * Rather, what is required is a way to map this data from the reference cell
923 * to the real cell. Classes that can do that are derived from the Mapping
924 * class, though one again often does not have to deal with them directly:
925 * many functions in the library can take a mapping object as argument, but
926 * when it is omitted they simply resort to the standard bilinear Q1
927 * mapping. We will go this route, and not bother with it for the moment (we
928 * come back to this in @ref step_10 "step-10", @ref step_11 "step-11", and @ref step_12 "step-12").
929 *
930
931 *
932 * So what we now have is a collection of three classes to deal with: finite
933 * element, quadrature, and mapping objects. That's too much, so there is one
936 * (or two, and an implicit linear mapping), it will be able to provide you
938 * quadrature points on a real cell.
939 *
940
941 *
943 * following function:
944 *
945 * @code
947 *   {
948 * @endcode
949 *
950 * Ok, let's start: we need a quadrature formula for the evaluation of the
951 * integrals on each cell. Let's take a Gauss formula with two quadrature
952 * points in each direction, i.e. a total of four points since we are in
953 * 2d. This quadrature formula integrates polynomials of degrees up to three
954 * exactly (in 1d). It is easy to check that this is sufficient for the
956 *
957 * @code
958 *   const QGauss<2> quadrature_formula(fe.degree + 1);
959 * @endcode
960 *
961 * And we initialize the object which we have briefly talked about above. It
962 * needs to be told which finite element we want to use, and the quadrature
963 * points and their weights (jointly described by a Quadrature object). As
966 * compute on each cell: we need the values of the shape functions at the
967 * quadrature points (for the right hand side @f$(\varphi_i,f)@f$), their
968 * gradients (for the matrix entries @f$(\nabla \varphi_i, \nabla
969 * \varphi_j)@f$), and also the weights of the quadrature points and the
971 * the real cells.
972 *
973
974 *
975 * This list of what kind of information we actually need is given as a
978 * time we go to a new cell, all of these flags start with the prefix
980 * updated. The flag to give if we want the values of the shape functions
986 *
987 * @code
988 *   FEValues<2> fe_values(fe,
991 * @endcode
992 *
993 * The advantage of this approach is that we can specify what kind of
997 * normal vectors to cells, etc are computed on each cell, regardless of
999 *
1000
1001 *
1005 * <code>operator|</code> is the <i>bitwise or operator</i>, i.e.,
1007 * patterns and returns an integer in which every bit is set for
1009 * arguments. For example, consider the operation
1010 * <code>9|10</code>. In binary, <code>9=0b1001</code> (where the
1012 * interpreted as a binary number) and <code>10=0b1010</code>. Going
1013 * through each bit and seeing whether it is set in one of the
1014 * argument, we arrive at <code>0b1001|0b1010=0b1011</code> or, in
1015 * decimal notation, <code>9|10=11</code>. The second piece of
1019 * <code>update_values=0b00001=1</code>,
1020 * <code>update_gradients=0b00010=2</code>,
1021 * <code>update_JxW_values=0b10000=16</code>. Then
1023 * 0b10011 = 19</code>. In other words, we obtain a number that
1026 * one bit in the integer that, if equal to one, means that a
1028 * zero, means that we need not compute it. In other words, even
1034 *
1035
1036 *
1037 * For use further down below, we define a shortcut for a value that will
1038 * be used very frequently. Namely, an abbreviation for the number of degrees
1039 * of freedom on each cell (since we are in 2d and degrees of freedom are
1040 * associated with vertices only, this number is four, but we rather want to
1041 * write the definition of this variable in a way that does not preclude us
1043 * number of degrees of freedom per cell, or work in a different space
1044 * dimension).
1045 *
1046
1047 *
1050 * you may want to change the finite element at some time. Changing the
1051 * element would have to be done in a different function and it is easy
1054 * the right object for the information: Here, we ask the finite element
1055 * to tell us about the number of degrees of freedom per cell and we
1057 * polynomial degree we may have chosen elsewhere in the program.
1058 *
1059
1060 *
1064 * larger programs, and `dofs_per_cell` is one that is more or less the
1065 * conventional name for this kind of object.
1066 *
1067 * @code
1068 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1069 *  
1070 * @endcode
1071 *
1072 * Now, we said that we wanted to assemble the global matrix and vector
1073 * cell-by-cell. We could write the results directly into the global matrix,
1074 * but this is not very efficient since access to the elements of a sparse
1075 * matrix is slow. Rather, we first compute the contribution of each cell in
1076 * a small matrix with the degrees of freedom on the present cell, and only
1077 * transfer them to the global matrix when the computations are finished for
1078 * this cell. We do the same for the right hand side vector. So let's first
1079 * allocate these objects (these being local objects, all degrees of freedom
1080 * are coupling with all others, and we should use a full matrix object
1081 * rather than a sparse one for the local operations; everything will be
1082 * transferred to a global sparse matrix later on):
1083 *
1084 * @code
1085 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1086 *   Vector<double> cell_rhs(dofs_per_cell);
1087 *  
1088 * @endcode
1089 *
1090 * When assembling the contributions of each cell, we do this with the local
1091 * numbering of the degrees of freedom (i.e. the number running from zero
1092 * through dofs_per_cell-1). However, when we transfer the result into the
1093 * global matrix, we have to know the global numbers of the degrees of
1094 * freedom. When we query them, we need a scratch (temporary) array for
1095 * these numbers (see the discussion at the end of the introduction for
1096 * the type, types::global_dof_index, used here):
1097 *
1098 * @code
1099 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1100 *  
1101 * @endcode
1102 *
1103 * Now for the loop over all cells. We have seen before how this works for a
1104 * triangulation. A DoFHandler has cell iterators that are exactly analogous
1105 * to those of a Triangulation, but with extra information about the degrees
1106 * of freedom for the finite element you're using. Looping over the active
1107 * cells of a degree-of-freedom handler works the same as for a triangulation.
1108 *
1109
1110 *
1111 * Note that we declare the type of the cell as `const auto &` instead of
1112 * `auto` this time around. In step 1, we were modifying the cells of the
1114 * examining the cells without modifying them, so it's good practice to
1115 * declare `cell` as `const` in order to enforce this invariant.
1116 *
1117 * @code
1118 *   for (const auto &cell : dof_handler.active_cell_iterators())
1119 *   {
1120 * @endcode
1121 *
1122 * We are now sitting on one cell, and we would like the values and
1124 * determinants of the Jacobian matrices of the mapping between
1125 * reference cell and true cell, at the quadrature points. Since all
1126 * these values depend on the geometry of the cell, we have to have the
1127 * FEValues object re-compute them on each cell:
1128 *
1129 * @code
1130 *   fe_values.reinit(cell);
1131 *  
1132 * @endcode
1133 *
1134 * Next, reset the local cell's contributions to global matrix and
1135 * global right hand side to zero, before we fill them:
1136 *
1137 * @code
1138 *   cell_matrix = 0;
1139 *   cell_rhs = 0;
1140 *  
1141 * @endcode
1142 *
1143 * Now it is time to start integration over the cell, which we
1144 * do by looping over all quadrature points, which we will
1145 * number by q_index.
1146 *
1147 * @code
1148 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
1149 *   {
1150 * @endcode
1151 *
1152 * First assemble the matrix: For the Laplace problem, the
1153 * matrix on each cell is the integral over the gradients of
1154 * shape function i and j. Since we do not integrate, but
1155 * rather use quadrature, this is the sum over all
1156 * quadrature points of the integrands times the determinant
1157 * of the Jacobian matrix at the quadrature point times the
1158 * weight of this quadrature point. You can get the gradient
1159 * of shape function @f$i@f$ at quadrature point with number q_index by
1160 * using <code>fe_values.shape_grad(i,q_index)</code>; this
1161 * gradient is a 2-dimensional vector (in fact it is of type
1162 * Tensor@<1,dim@>, with here dim=2) and the product of two
1163 * such vectors is the scalar product, i.e. the product of
1164 * the two shape_grad function calls is the dot
1165 * product. This is in turn multiplied by the Jacobian
1166 * determinant and the quadrature point weight (that one
1167 * gets together by the call to FEValues::JxW() ). Finally,
1168 * this is repeated for all shape functions @f$i@f$ and @f$j@f$:
1169 *
1170 * @code
1171 *   for (const unsigned int i : fe_values.dof_indices())
1172 *   for (const unsigned int j : fe_values.dof_indices())
1173 *   cell_matrix(i, j) +=
1174 *   (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
1175 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
1176 *   fe_values.JxW(q_index)); // dx
1177 *  
1178 * @endcode
1179 *
1180 * We then do the same thing for the right hand side. Here,
1181 * the integral is over the shape function i times the right
1182 * hand side function, which we choose to be the function
1183 * with constant value one (more interesting examples will
1184 * be considered in the following programs).
1185 *
1186 * @code
1187 *   for (const unsigned int i : fe_values.dof_indices())
1188 *   cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1189 *   1. * // f(x_q)
1190 *   fe_values.JxW(q_index)); // dx
1191 *   }
1192 * @endcode
1193 *
1194 * Now that we have the contribution of this cell, we have to transfer
1195 * it to the global matrix and right hand side. To this end, we first
1196 * have to find out which global numbers the degrees of freedom on this
1197 * cell have. Let's simply ask the cell for that information:
1198 *
1199 * @code
1200 *   cell->get_dof_indices(local_dof_indices);
1201 *  
1202 * @endcode
1203 *
1204 * Then again loop over all shape functions i and j and transfer the
1205 * local elements to the global matrix. The global numbers can be
1206 * obtained using local_dof_indices[i]:
1207 *
1208 * @code
1209 *   for (const unsigned int i : fe_values.dof_indices())
1210 *   for (const unsigned int j : fe_values.dof_indices())
1211 *   system_matrix.add(local_dof_indices[i],
1212 *   local_dof_indices[j],
1213 *   cell_matrix(i, j));
1214 *  
1215 * @endcode
1216 *
1217 * And again, we do the same thing for the right hand side vector.
1218 *
1219 * @code
1220 *   for (const unsigned int i : fe_values.dof_indices())
1221 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
1222 *   }
1223 *  
1224 *  
1225 * @endcode
1226 *
1227 * Now almost everything is set up for the solution of the discrete
1228 * system. However, we have not yet taken care of boundary values (in fact,
1229 * Laplace's equation without Dirichlet boundary values is not even uniquely
1230 * solvable, since you can add an arbitrary constant to the discrete
1231 * solution). We therefore have to do something about the situation.
1232 *
1233
1234 *
1235 * For this, we first obtain a list of the degrees of freedom on the
1236 * boundary and the value the shape function shall have there. For
1237 * simplicity, we only interpolate the boundary value function, rather than
1238 * projecting it onto the boundary. There is a function in the library which
1239 * does exactly this: VectorTools::interpolate_boundary_values(). Its
1240 * parameters are (omitting parameters for which default values exist and
1241 * that we don't care about): the DoFHandler object to get the global
1242 * numbers of the degrees of freedom on the boundary; the component of the
1243 * boundary where the boundary values shall be interpolated; the boundary
1244 * value function itself; and the output object.
1245 *
1246
1247 *
1248 * The component of the boundary is meant as follows: in many cases, you may
1249 * want to impose certain boundary values only on parts of the boundary. For
1254 * function to only compute the boundary values on a certain part of the
1255 * boundary (e.g. the clamped part, or the inflow boundary). By default,
1256 * all boundaries have a 0 boundary indicator, unless otherwise specified.
1258 * If sections of the boundary have different boundary conditions, you have to
1259 * number those parts with different boundary indicators. The function call
1261 * boundary for which the boundary indicator is in fact the zero specified as
1263 *
1264
1265 *
1266 * The function describing the boundary values is an object of type Function
1267 * or of a derived class. One of the derived classes is
1269 * which is zero everywhere. We create such an object in-place and pass it to
1271 *
1272
1273 *
1274 * Finally, the output object is a list of pairs of global degree of freedom
1275 * numbers (i.e. the number of the degrees of freedom on the boundary) and
1276 * their boundary values (which are zero here for all entries). This mapping
1277 * of DoF numbers to boundary values is done by the <code>std::map</code>
1278 * class.
1279 *
1280 * @code
1281 *   std::map<types::global_dof_index, double> boundary_values;
1282 *   VectorTools::interpolate_boundary_values(dof_handler,
1283 *   types::boundary_id(0),
1284 *   Functions::ZeroFunction<2>(),
1286 * @endcode
1287 *
1288 * Now that we got the list of boundary DoFs and their respective boundary
1289 * values, let's use them to modify the system of equations
1290 * accordingly. This is done by the following function call:
1291 *
1292 * @code
1293 *   MatrixTools::apply_boundary_values(boundary_values,
1294 *   system_matrix,
1295 *   solution,
1296 *   system_rhs);
1297 *   }
1298 *  
1299 *  
1300 * @endcode
1301 *
1302 *
1303 * <a name="step_3-Step3solve"></a>
1304 * <h4>Step3::solve</h4>
1305 *
1306
1307 *
1310 * specifically the Conjugate Gradient (CG) method.
1311 *
1312
1313 *
1314 * The way to do this in deal.II is a three-step process:
1315 * - First, we need to have an object that knows how to tell the CG algorithm
1316 * when to stop. This is done by using a SolverControl object, and as
1317 * stopping criterion we say: stop after a maximum of 1000 iterations (which
1318 * is far more than is needed for 1089 variables; see the results section to
1319 * find out how many were really used), and stop if the norm of the residual
1320 * is below @f$\tau=10^{-6}\|\mathbf b\|@f$ where @f$\mathbf b@f$ is the right hand
1321 * side vector. In practice, this latter criterion will be the one
1323 * - Then we need the solver itself. The template parameter to the SolverCG
1324 * class is the type of the vectors we are using.
1325 * - The last step is to actually solve the system of equations. The CG solver
1326 * takes as arguments the components of the linear system @f$Ax=b@f$ (in the
1327 * order in which they appear in this equation), and a preconditioner
1328 * as the fourth argument. We don't feel ready to delve into preconditioners
1329 * yet, so we tell it to use the identity operation as preconditioner. Later
1330 * tutorial programs will spend significant amount of time and space on
1331 * constructing better preconditioners.
1332 *
1333
1334 *
1335 * At the end of this process, the `solution` variable contains the
1336 * nodal values of the solution function. At the end of the function, we
1337 * output how many Conjugate Gradients iterations it took to solve the
1338 * linear system.
1339 *
1340 * @code
1341 *   void Step3::solve()
1342 *   {
1343 *   SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
1344 *   SolverCG<Vector<double>> solver(solver_control);
1345 *   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
1346 *  
1347 *   std::cout << solver_control.last_step()
1348 *   << " CG iterations needed to obtain convergence." << std::endl;
1349 *   }
1350 *  
1351 *  
1352 * @endcode
1353 *
1354 *
1355 * <a name="step_3-Step3output_results"></a>
1356 * <h4>Step3::output_results</h4>
1357 *
1358
1359 *
1360 * The last part of a typical finite element program is to output the results
1361 * and maybe do some postprocessing (for example compute the maximal stress
1362 * values at the boundary, or the average flux across the outflow, etc). We
1363 * have no such postprocessing here, but we would like to write the solution
1364 * to a file.
1365 *
1366 * @code
1367 *   void Step3::output_results() const
1368 *   {
1369 * @endcode
1370 *
1371 * To write the output to a file, we need an object which knows about output
1372 * formats and the like. This is the DataOut class, and we need an object of
1373 * that type:
1374 *
1375 * @code
1376 *   DataOut<2> data_out;
1377 * @endcode
1378 *
1379 * Now we have to tell it where to take the values from which it shall
1380 * write. We tell it which DoFHandler object to use, and the solution vector
1381 * (and the name by which the solution variable shall appear in the output
1382 * file). If we had more than one vector which we would like to look at in
1383 * the output (for example right hand sides, errors per cell, etc) we would
1384 * add them as well:
1385 *
1386 * @code
1387 *   data_out.attach_dof_handler(dof_handler);
1388 *   data_out.add_data_vector(solution, "solution");
1389 * @endcode
1390 *
1391 * After the DataOut object knows which data it is to work on, we have to
1392 * tell it to process them into something the back ends can handle. The
1393 * reason is that we have separated the frontend (which knows about how to
1394 * treat DoFHandler objects and data vectors) from the back end (which knows
1395 * many different output formats) and use an intermediate data format to
1396 * transfer data from the front- to the backend. The data is transformed
1397 * into this intermediate format by the following function:
1398 *
1399 * @code
1400 *   data_out.build_patches();
1401 *  
1402 * @endcode
1403 *
1404 * Now we have everything in place for the actual output. Just open a file
1405 * and write the data into it, using VTK format (there are many other
1406 * functions in the DataOut class we are using here that can write the
1407 * data in postscript, AVS, GMV, Gnuplot, or some other file
1408 * formats):
1409 *
1410 * @code
1411 *   const std::string filename = "solution.vtk";
1412 *   std::ofstream output(filename);
1413 *   data_out.write_vtk(output);
1414 *   std::cout << "Output written to " << filename << std::endl;
1415 *   }
1416 *  
1417 *  
1418 * @endcode
1419 *
1420 *
1421 * <a name="step_3-Step3run"></a>
1422 * <h4>Step3::run</h4>
1423 *
1424
1425 *
1426 * Finally, the last function of this class is the main function which calls
1427 * all the other functions of the <code>Step3</code> class. The order in which
1428 * this is done resembles the order in which most finite element programs
1429 * work. Since the names are mostly self-explanatory, there is not much to
1430 * comment about:
1431 *
1432 * @code
1433 *   void Step3::run()
1434 *   {
1435 *   make_grid();
1436 *   setup_system();
1437 *   assemble_system();
1438 *   solve();
1439 *   output_results();
1440 *   }
1441 *  
1442 *  
1443 * @endcode
1444 *
1445 *
1446 * <a name="step_3-Thecodemaincodefunction"></a>
1447 * <h3>The <code>main</code> function</h3>
1448 *
1449
1450 *
1451 * This is the main function of the program. Since the concept of a
1452 * main function is mostly a remnant from the pre-object oriented era
1453 * before C++ programming, it often does not do much more than
1454 * creating an object of the top-level class and calling its principle
1455 * function.
1456 *
1457 * @code
1458 *   int main()
1459 *   {
1460 *   Step3 laplace_problem;
1461 *   laplace_problem.run();
1462 *  
1463 *   return 0;
1464 *   }
1465 * @endcode
1466<a name="step_3-Results"></a><h1>Results</h1>
1467
1468
1469The output of the program looks as follows:
1470@code
1471Number of active cells: 1024
1472Number of degrees of freedom: 1089
147336 CG iterations needed to obtain convergence.
1474Output written to solution.vtk
1475@endcode
1476
1477The last line is the output we generated at the bottom of the
1478`output_results()` function: The program generated the file
1479<code>solution.vtk</code>, which is in the VTK format that is widely
1480used by many visualization programs today -- including the two
1481heavy-weights <a href="https://www.llnl.gov/visit">VisIt</a> and
1482<a href="https://www.paraview.org">Paraview</a> that are the most
1483commonly used programs for this purpose today.
1484
1485Using VisIt, it is not very difficult to generate a picture of the
1486solution like this:
1487<table width="60%" align="center">
1488 <tr>
1489 <td align="center">
1490 <img src="https://www.dealii.org/images/steps/developer/step-3.solution-3.png" alt="Visualization of the solution of step-3">
1491 </td>
1492 </tr>
1493</table>
1494It shows both the solution and the mesh, elevated above the @f$x@f$-@f$y@f$ plane
1495based on the value of the solution at each point. Of course the solution
1496here is not particularly exciting, but that is a result of both what the
1497Laplace equation represents and the right hand side @f$f(\mathbf x)=1@f$ we
1498have chosen for this program: The Laplace equation describes (among many
1499other uses) the vertical deformation of a membrane subject to an external
1500(also vertical) force. In the current example, the membrane's borders
1504
1507See also <a href="https://www.math.colostate.edu/~bangerth/videos.676.11.html">video lecture 11</a>, <a href="https://www.math.colostate.edu/~bangerth/videos.676.32.html">video lecture 32</a>.
1508
1509
1510
1511<a name="step-3-extensions"></a>
1512<a name="step_3-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1513
1514
1517</p>
1518
1519<ul>
1520 <li>
1521 Change the geometry and mesh: In the program, we have generated a square
1523 function. However, the <code>GridGenerator</code> has a good number of other
1524 functions as well. Try an L-shaped domain, a ring, or other domains you find
1525 there.
1526 </li>
1527
1528 <li>
1529 Change the boundary condition: The code uses the Functions::ZeroFunction
1530 function to generate zero boundary conditions. However, you may want to try
1531 non-zero constant boundary values using
1532 <code>Functions::ConstantFunction&lt;2&gt;(1)</code> instead of
1533 <code>Functions::ZeroFunction&lt;2&gt;()</code> to have unit Dirichlet boundary
1534 values. More exotic functions are described in the documentation of the
1536 values.
1537 </li>
1538
1539 <li> Modify the type of boundary condition: Presently, what happens
1540 is that we use Dirichlet boundary values all around, since the
1541 default is that all boundary parts have boundary indicator zero, and
1542 then we tell the
1543 VectorTools::interpolate_boundary_values() function to
1544 interpolate boundary values to zero on all boundary components with
1545 indicator zero. <p> We can change this behavior if we assign parts
1548 @code
1549 triangulation.begin_active()->face(0)->set_boundary_id(1);
1550 @endcode
1551
1553 return an iterator that points to the first active cell. Of course,
1554 this being the coarse mesh for the triangulation of a square, the
1556 active. Next, we ask the cell to return an iterator to its first
1557 face, and then we ask the face to reset the boundary indicator of
1559 faces of child cells inherit the boundary indicator of their
1560 parents, i.e. even on the finest mesh, the faces on one side of the
1561 square have boundary indicator 1. Later, when we get to
1562 interpolating boundary conditions, the
1563 VectorTools::interpolate_boundary_values() call will only produce boundary
1564 values for those faces that have zero boundary indicator, and leave
1565 those faces alone that have a different boundary indicator. What
1568 normal derivative of the solution, unless one adds additional terms
1571 run the program.
1572
1575 For example, we can label all of the cells along the top and
1577 see if the cell centers' y-coordinates are within a tolerance
1578 (here `1e-12`) of -1 and 1. Try this immediately after calling
1580 @code
1581 for (auto &face : triangulation.active_face_iterators())
1582 if (face->at_boundary())
1583 if (std::fabs(face->center()(1) - (-1.0)) < 1e-12 ||
1584 std::fabs(face->center()(1) - (1.0)) < 1e-12)
1585 face->set_boundary_id(1);
1586 @endcode
1589
1590 <li>
1591 A slight variation of the last point would be to set different boundary
1592 values as above, but then use a different boundary value function for
1593 boundary indicator one. In practice, what you have to do is to add a second
1594 call to <code>interpolate_boundary_values</code> for boundary indicator one:
1595 @code
1596 VectorTools::interpolate_boundary_values(dof_handler,
1597 types::boundary_id(1),
1598 Functions::ConstantFunction<2>(1.),
1600 @endcode
1601 If you have this call immediately after the first one to this function, then
1602 it will interpolate boundary values on faces with boundary indicator 1 to the
1603 unit value, and merge these interpolated values with those previously
1604 computed for boundary indicator 0. The result will be that we will get
1605 discontinuous boundary values, zero on three sides of the square, and one on
1606 the fourth.
1607
1608 <li>
1609 Use triangles: As mentioned in the results section of @ref step_1 "step-1", for
1614
1615 This is *almost* trivial. First, as discussed in @ref step_1 "step-1", we may want
1619 @code
1624 @endcode
1631 of @ref step_4 "step-4", but that goes beyond what we want to demonstrate here.
1632
1635 @code
1636--------------------------------------------------------
1637An error occurred in line <2633> of file </home/bangerth/p/deal.II/1/dealii/include/deal.II/dofs/dof_accessor.templates.h> in function
1642 The reference-cell type used on this cell (Tri) does not match the
1643 reference-cell type of the finite element associated with this cell
1646 a simplex element to a hypercube cell (or the other way around) via
1647 the active_fe_index?
1648 @endcode
1654
1657 you will find that the FE_Q element we use in the main class can
1658 only be used on hypercube meshes; what we *want* to use instead now
1661 have to add `#include <deal.II/fe/fe_simplex_p.h>` at the top of the file.)
1662
1663 The last thing you need to change (which at the time of writing is
1665 we integrate, we need to use a quadrature formula that is
1668
1669 With all of these steps, you then get the following solution:
1670 <img src="https://www.dealii.org/images/steps/developer/step-3.solution-triangles.png" alt="Visualization of the solution of step-3 using triangles">
1671
1672 <li>
1674 @ref step_7 "step-7", but it is easy to check that computations converge
1675 already here. For example, we could evaluate the value of the solution in a
1676 single point and compare the value for different %numbers of global
1677 refinement (the number of global refinement steps is set in
1679 solution at a point, say at @f$(\frac 13, \frac 13)@f$, we could add the
1681 @code
1682 std::cout << "Solution at (1/3,1/3): "
1683 << VectorTools::point_value(dof_handler, solution,
1684 Point<2>(1./3, 1./3))
1685 << std::endl;
1686 @endcode
1687 For 1 through 9 global refinement steps, we then get the following sequence
1688 of point values:
1689 <table align="center" class="doxtable">
1690 <tr> <th># of refinements</th> <th>@f$u_h(\frac 13,\frac13)@f$</th> </tr>
1691 <tr> <td>1</td> <td>0.166667</td> </tr>
1692 <tr> <td>2</td> <td>0.227381</td> </tr>
1693 <tr> <td>3</td> <td>0.237375</td> </tr>
1694 <tr> <td>4</td> <td>0.240435</td> </tr>
1695 <tr> <td>5</td> <td>0.241140</td> </tr>
1696 <tr> <td>6</td> <td>0.241324</td> </tr>
1697 <tr> <td>7</td> <td>0.241369</td> </tr>
1698 <tr> <td>8</td> <td>0.241380</td> </tr>
1699 <tr> <td>9</td> <td>0.241383</td> </tr>
1700 </table>
1702 by about a factor of 4, we can conjecture that the "correct" value may be
1703 @f$u(\frac 13, \frac 13)\approx 0.241384@f$. In fact, if we assumed this to be
1704 the correct value, we could show that the sequence above indeed shows @f${\cal
1708
1709 A slight variant of this would be to repeat the test with quadratic
1710 elements. All you need to do is to set the polynomial degree of the finite
1711 element to two in the constructor
1713
1714 <li>Convergence of the mean: A different way to see that the solution
1715 actually converges (to something &mdash; we can't tell whether it's really
1716 the correct value!) is to compute the mean of the solution. To this end, add
1718 @code
1719 std::cout << "Mean value: "
1720 << VectorTools::compute_mean_value (dof_handler,
1721 QGauss<2>(fe.degree + 1),
1722 solution,
1723 0)
1724 << std::endl;
1725 @endcode
1727 parameters mean, while the first and third should be obvious. Doing the same
1728 study again where we change the number of global refinement steps, we get
1730 <table align="center" class="doxtable">
1731 <tr> <th># of refinements</th> <th>@f$\int_\Omega u_h(x)\; dx@f$</th> </tr>
1732 <tr> <td>0</td> <td>0.09375000</td> </tr>
1733 <tr> <td>1</td> <td>0.12790179</td> </tr>
1734 <tr> <td>2</td> <td>0.13733440</td> </tr>
1735 <tr> <td>3</td> <td>0.13976069</td> </tr>
1736 <tr> <td>4</td> <td>0.14037251</td> </tr>
1737 <tr> <td>5</td> <td>0.14052586</td> </tr>
1738 <tr> <td>6</td> <td>0.14056422</td> </tr>
1739 <tr> <td>7</td> <td>0.14057382</td> </tr>
1740 <tr> <td>8</td> <td>0.14057622</td> </tr>
1741 </table>
1743 factor of four, indicating convergence as @f${\cal O}(h^2)@f$.
1744</ul>
1745
1746
1747
1748<a name="step_3-UsingHDF5tooutputthesolutionandadditionaldata"></a><h3>Using %HDF5 to output the solution and additional data</h3>
1749
1750
1752languages (e.g. R or Python). It is not difficult to get deal.II to
1754postprocess some of the data generated by this program. Here are some
1755ideas on what is possible.
1756
1757
1758<a name="step_3-Changingtheoutputtoh5"></a><h4> Changing the output to .h5</h4>
1759
1760
1761To fully make use of the automation we first need to introduce a private variable for the number of
1762global refinement steps <code>unsigned int n_refinement_steps </code>, which will be used for the output filename.
1763In <code>make_grid()</code> we then replace <code>triangulation.refine_global(5);</code> with
1764@code
1765n_refinement_steps = 5;
1766triangulation.refine_global(n_refinement_steps);
1767@endcode
1769namespace (for interfacing to general-purpose data files)
1772Although the HDF5 deal.II binding supports both serial and MPI, the %HDF5 DataOut binding
1773only supports parallel output.
1774For this reason we need to initialize an MPI
1776@code
1777int main(int argc, char* argv[])
1778{
1780 ...
1781}
1782@endcode
1785@code
1786const std::string filename_h5 = "solution_" + std::to_string(n_refinement_steps) + ".h5";
1787DataOutBase::DataOutFilterFlags flags(true, true);
1789data_out.write_filtered_data(data_filter);
1790data_out.write_hdf5_parallel(data_filter, filename_h5, MPI_COMM_WORLD);
1791@endcode
1796
1797
1798<a name="step_3-Addingthepointvalueandthemeanseeextensionaboveintotheh5file"></a><h4> Adding the point value and the mean (see extension above) into the .h5 file</h4>
1799
1800
1803of our experiment in a single result file, which can then be read and
1807
1809@code
1810#include <deal.II/base/hdf5.h>
1811@endcode
1814solution at a particular point, as well as the mean value of the
1815solution, to our %HDF5 file :
1816@code
1819point_value[0] = VectorTools::point_value(dof_handler, solution,
1820 Point<2>(1./3, 1./3));
1821data_file.write_dataset("point_value", point_value);
1822Vector<double> mean_value(1);
1823mean_value[0] = VectorTools::compute_mean_value(dof_handler,
1824 QGauss<2>(fe.degree + 1),
1825 solution, 0);
1826data_file.write_dataset("mean_value",mean_value);
1827@endcode
1828
1829
1830
1831<a name="step_3-UsingRandggplot2togenerateplots"></a><h3> Using R and ggplot2 to generate plots</h3>
1832
1834
1837how this can, in particular, be done with the
1838<a href="https://en.wikipedia.org/wiki/R_(programming_language)">R
1842<a href="https://datacarpentry.org/R-ecology-lesson/index.html">here</a>.
1843Furthermore, since most search engines struggle with searches of the form "R + topic",
1845href="http://rseek.org">RSeek </a> instead.
1846
1847The most prominent difference between R and other languages is that
1850To open the `.h5` file in R you have to install the <a href="https://bioconductor.org/packages/release/bioc/html/rhdf5.html">rhdf5</a> package, which is a part of the Bioconductor package.
1851
1853@code{.r}
1854library(rhdf5) # library for handling HDF5 files
1855library(ggplot2) # main plotting library
1856library(grDevices) # needed for output to PDF
1857library(viridis) # contains good colormaps for sequential data
1858
1859refinement <- 5
1860h5f <- H5Fopen(paste("solution_",refinement,".h5",sep=""))
1861print(h5f)
1862@endcode
1863This gives the following output
1864@code{.unparsed}
1865HDF5 FILE
1866 name /
1868
1869 name otype dclass dim
18700 cells H5I_DATASET INTEGER x 1024
18711 mean_value H5I_DATASET FLOAT 1
18722 nodes H5I_DATASET FLOAT x 1089
18744 solution H5I_DATASET FLOAT x 1089
1875@endcode
1877<code>dim(h5f\$cells)</code> gives us the dimensions of the matrix
1878that is used to store our cells.
1879We can see the following three matrices, as well as the two
1880additional data points we added.
1881<ul>
1882<li> <code>cells</code>: a 4x1024 matrix that stores the (C++) vertex indices for each cell
1883<li> <code>nodes</code>: a 2x1089 matrix storing the position values (x,y) for our cell vertices
1884<li> <code>solution</code>: a 1x1089 matrix storing the values of our solution at each vertex
1885</ul>
1889
1890<code>nodes</code> and <code>cells</code> contain all the information we need to plot our grid.
1892@code{.r}
1893# Counting in R starts at 1 instead of 0, so we need to increment all
1894# vertex indices by one:
1895cell_ids <- h5f@f$cells+1
1896
1897# Store the x and y positions of each vertex in one big vector in a
1898# cell by cell fashion (every 4 entries belong to one cell):
1899cells_x <- h5f@f$nodes[1,][cell_ids]
1900cells_y <- h5f@f$nodes[2,][cell_ids]
1901
1902# Construct a vector that stores the matching cell by cell grouping
1903# (1,1,1,1,2,2,2,2,...):
1904groups <- rep(1:ncol(cell_ids),each=4)
1905
1907meshdata <- data.frame(x = cells_x, y = cells_y, id = groups)
1908@endcode
1909
1911@code{.r}
1912pdf (paste("grid_",refinement,".pdf",sep=""),width = 5,height = 5) # Open new PDF file
1914 # object, at first only data
1915
1916plt <- plt + geom_polygon(fill="white",colour="black") # Actual plotting of the grid as polygons
1917plt <- plt + ggtitle(paste("grid at refinement level #",refinement))
1918
1919print(plt) # Show the current state of the plot/add it to the pdf
1920dev.off() # Close PDF file
1921@endcode
1922
1924you get the idea):
1925<table width="60%" align="center">
1926 <tr>
1927 <td align="center">
1928 <img src="https://www.dealii.org/images/steps/developer/step-3.extensions.grid_5.png" alt="Grid after 5 refinement steps of step-3">
1929 </td>
1930 </tr>
1931</table>
1932
1936This function needs a structured grid, i.e. uniform in x and y directions.
1939@code{.r}
1940pdf (paste("pseudocolor_",refinement,".pdf",sep=""),width = 5,height = 4.2) # Open new PDF file
1941colordata <- data.frame(x = h5f@f$nodes[1,],y = h5f@f$nodes[2,] , solution = h5f@f$solution[1,])
1942plt <- ggplot(colordata,aes(x=x,y=y,fill=solution))
1945plt <- plt + ggtitle(paste("solution at refinement level #",refinement))
1946
1947print(plt)
1948dev.off()
1949H5Fclose(h5f) # Close the HDF5 file
1950@endcode
1952<table width="60%" align="center">
1953 <tr>
1954 <td align="center">
1955 <img src="https://www.dealii.org/images/steps/developer/step-3.extensions.pseudocolor_5.png" alt="Solution after 5 refinement steps of step-3">
1956 </td>
1957 </tr>
1958</table>
1959
1961starting from 1.
1963@code{.r}
1965
1966# First we initiate all vectors with the results of the first level
1967h5f <- H5Fopen("solution_1.h5")
1968dofs <- dim(h5f@f$solution)[2]
1972
1973for (reflevel in 2:n_ref)
1974{
1975 h5f <- H5Fopen(paste("solution_",reflevel,".h5",sep=""))
1976 dofs <- c(dofs,dim(h5f\$solution)[2])
1977 mean <- c(mean,h5f\$mean_value)
1978 point <- c(point,h5f\$point_value)
1979 H5Fclose(h5f)
1980}
1981@endcode
1984@code{.r}
1985# Calculate the error w.r.t. our maximum refinement step
1986mean_error <- abs(mean[1:n_ref-1]-mean[n_ref])
1988
1990dofs <- dofs[1:n_ref-1]
1991convdata <- data.frame(dofs = dofs, mean_value= mean_error, point_value = point_error)
1992@endcode
1993Now we have all the data available to generate our plots.
1996@code
1997pdf (paste("convergence.pdf",sep=""),width = 5,height = 4.2)
1998plt <- ggplot(convdata,mapping=aes(x = dofs, y = mean_value))
1999plt <- plt+geom_line()
2000plt <- plt+labs(x="#DoFs",y = "mean value error")
2002print(plt)
2003
2004plt <- ggplot(convdata,mapping=aes(x = dofs, y = point_value))
2005plt <- plt+geom_line()
2006plt <- plt+labs(x="#DoFs",y = "point value error")
2008print(plt)
2009
2010dev.off()
2011@endcode
2014to zero:
2015<table style="width:50%" align="center">
2016 <tr>
2017 <td><img src="https://www.dealii.org/images/steps/developer/step-3.extensions.convergence_mean.png" alt=""></td>
2018 <td><img src="https://www.dealii.org/images/steps/developer/step-3.extensions.convergence_point.png" alt=""></td>
2019 </tr>
2020</table>
2021
2022<a name="step_3-UsingPythontogenerateplots"></a><h3> Using Python to generate plots</h3>
2023
2024
2027@code{.py}
2028import numpy as np # to work with multidimensional arrays
2029import h5py # to work with %HDF5 files
2030
2031import pandas as pd # for data frames
2033from matplotlib.patches import Polygon
2034
2037
2038@endcode
2039The %HDF5 solution file corresponding to `refinement = 5` can be opened as
2040@code{.py}
2041refinement = 5
2042filename = "solution_%d.h5" % refinement
2043file = h5py.File(filename, "r")
2044@endcode
2045The following prints out the tables in the %HDF5 file
2046@code{.py}
2047for key, value in file.items():
2048 print(key, " : ", value)
2049@endcode
2050which prints out
2051@code{.unparsed}
2052cells : <HDF5 dataset "cells": shape (1024, 4), type "<u4">
2053mean_value : <HDF5 dataset "mean_value": shape (1,), type "<f8">
2054nodes : <HDF5 dataset "nodes": shape (1089, 2), type "<f8">
2055point_value : <HDF5 dataset "point_value": shape (1,), type "<f8">
2056solution : <HDF5 dataset "solution": shape (1089, 1), type "<f8">
2057@endcode
2058There are @f$(32+1)\times(32+1) = 1089@f$ nodes.
2059The coordinates of these nodes @f$(x,y)@f$ are stored in the table `nodes` in the %HDF5 file.
2060There are a total of @f$32\times 32 = 1024@f$ cells. The nodes which make up each cell are
2061marked in the table `cells` in the %HDF5 file.
2062
2064@code{.py}
2065nodes = np.array(file["/nodes"])
2066cells = np.array(file["/cells"])
2067solution = np.array(file["/solution"])
2068
2069x, y = nodes.T
2070@endcode
2071
2073@code{.py}
2074cell_x = x[cells.flatten()]
2075cell_y = y[cells.flatten()]
2076@endcode
2077The following tags the cell ids. Each four entries correspond to one cell.
2079@code{.py}
2080n_cells = cells.shape[0]
2081cell_ids = np.repeat(np.arange(n_cells), 4)
2082meshdata = pd.DataFrame({"x": cell_x, "y": cell_y, "ids": cell_ids})
2083@endcode
2085@code{.unparsed}
2086print(meshdata)
2087
2088 x y ids
20890 0.00000 0.00000 0
20901 0.03125 0.00000 0
20912 0.03125 0.03125 0
20923 0.00000 0.03125 0
20934 0.03125 0.00000 1
2094... ... ... ...
20954091 0.93750 1.00000 1022
20964092 0.96875 0.96875 1023
20974093 1.00000 0.96875 1023
20984094 1.00000 1.00000 1023
20994095 0.96875 1.00000 1023
2100
21014096 rows × 3 columns
2102@endcode
2103
2104To plot the mesh, we loop over all cells and connect the first and last node to complete the polygon
2105@code{.py}
2106fig, ax = plt.subplots()
2107ax.set_aspect("equal", "box")
2108ax.set_title("grid at refinement level #%d" % refinement)
2109
2110for cell_id, cell in meshdata.groupby(["ids"]):
2111 cell = pd.concat([cell, cell.head(1)])
2112 plt.plot(cell["x"], cell["y"], c="k")
2113@endcode
2115@code{.py}
2116fig, ax = plt.subplots()
2117ax.set_aspect("equal", "box")
2118ax.set_title("grid at refinement level #%d" % refinement)
2119for cell_id, cell in meshdata.groupby(["ids"]):
2120 patch = Polygon(cell[["x", "y"]], facecolor="w", edgecolor="k")
2122@endcode
2123
2124To plot the solution, we first create a finer grid and then interpolate the solution values
2125into the grid and then plot it.
2126@code{.py}
2127nx = int(np.sqrt(n_cells)) + 1
2128nx *= 10
2129xg = np.linspace(x.min(), x.max(), nx)
2130yg = np.linspace(y.min(), y.max(), nx)
2131
2132xgrid, ygrid = np.meshgrid(xg, yg)
2133solution_grid = griddata((x, y), solution.flatten(), (xgrid, ygrid), method="linear")
2134
2135fig = plt.figure()
2136ax = fig.add_subplot(1, 1, 1)
2137ax.set_title("solution at refinement level #%d" % refinement)
2138c = ax.pcolor(xgrid, ygrid, solution_grid, cmap=cm.viridis)
2139fig.colorbar(c, ax=ax)
2140
2141plt.show()
2142@endcode
2143
2144To check the convergence of `mean_value` and `point_value`
2145we loop over data of all refinements and store into vectors <code>mean_values</code> and <code>mean_values</code>
2146@code{.py}
2147mean_values = np.zeros((8,))
2148point_values = np.zeros((8,))
2149dofs = np.zeros((8,))
2150
2151for refinement in range(1, 9):
2152 filename = "solution_%d.h5" % refinement
2153 file = h5py.File(filename, "r")
2154 mean_values[refinement - 1] = np.array(file["/mean_value"])[0]
2155 point_values[refinement - 1] = np.array(file["/point_value"])[0]
2156 dofs[refinement - 1] = np.array(file["/solution"]).shape[0]
2157@endcode
2158
2159Following is the plot of <code>mean_values</code> on `log-log` scale
2160@code{.py}
2161mean_error = np.abs(mean_values[1:] - mean_values[:-1])
2162plt.loglog(dofs[:-1], mean_error)
2163plt.grid()
2164plt.xlabel("#DoFs")
2165plt.ylabel("mean value error")
2166plt.show()
2167@endcode
2168
2169Following is the plot of <code>point_values</code> on `log-log` scale
2170@code{.py}
2171point_error = np.abs(point_values[1:] - point_values[:-1])
2172plt.loglog(dofs[:-1], point_error)
2173plt.grid()
2174plt.xlabel("#DoFs")
2175plt.ylabel("point value error")
2176plt.show()
2177@endcode
2178
2179A Python package which mimics the `R` ggplot2 (which is based on specifying the grammar of the graphics) is
2180<a href="https://plotnine.org/">plotnine</a>.
2181@code{.py}
2182We need to import the following from the <code>plotnine</code> package
2183from plotnine import (
2184 ggplot,
2185 aes,
2186 geom_raster,
2187 geom_polygon,
2188 geom_line,
2189 labs,
2190 scale_x_log10,
2191 scale_y_log10,
2192 ggtitle,
2193)
2194@endcode
2195Then plot the mesh <code>meshdata</code> dataframe
2196@code{.py}
2197plot = (
2198 ggplot(meshdata, aes(x="x", y="y", group="ids"))
2199 + geom_polygon(fill="white", colour="black")
2200 + ggtitle("grid at refinement level #%d" % refinement)
2201)
2202
2203print(plot)
2204@endcode
2205Collect the solution into a dataframe
2206@code{.py}
2207colordata = pd.DataFrame({"x": x, "y": y, "solution": solution.flatten()})
2208@endcode
2209Plot of the solution
2210@code{.py}
2211plot = (
2212 ggplot(colordata, aes(x="x", y="y", fill="solution"))
2213 + geom_raster(interpolate=True)
2214 + ggtitle("solution at refinement level #%d" % refinement)
2215)
2216
2217print(plot)
2218@endcode
2219
2220Collect the convergence data into a dataframe
2221@code{.py}
2222convdata = pd.DataFrame(
2223 {"dofs": dofs[:-1], "mean_value": mean_error, "point_value": point_error}
2224)
2225
2226@endcode
2227Following is the plot of <code>mean_values</code> on `log-log` scale
2228@code{.py}
2229plot = (
2230 ggplot(convdata, mapping=aes(x="dofs", y="mean_value"))
2231 + geom_line()
2232 + labs(x="#DoFs", y="mean value error")
2233 + scale_x_log10()
2234 + scale_y_log10()
2235)
2236
2237plot.save("mean_error.pdf", dpi=60)
2238print(plot)
2239@endcode
2240
2241Following is the plot of <code>point_values</code> on `log-log` scale
2242@code{.py}
2243plot = (
2244 ggplot(convdata, mapping=aes(x="dofs", y="point_value"))
2245 + geom_line()
2246 + labs(x="#DoFs", y="point value error")
2247 + scale_x_log10()
2248 + scale_y_log10()
2249)
2250
2251plot.save("point_error.pdf", dpi=60)
2252print(plot)
2253@endcode
2254 *
2255 *
2256<a name="step_3-PlainProg"></a>
2257<h1> The plain program</h1>
2258@include "step-3.cc"
2259*/
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition fe_q.h:554
void write_dataset(const std::string &name, const Container &data) const
Definition hdf5.h:2244
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
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
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
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
std::vector< index_type > data
Definition mpi.cc:740
const Event initial
Definition event.cc:68
void interpolate(const DoFHandler< dim, spacedim > &dof1, const InVector &u1, const DoFHandler< dim, spacedim > &dof2, OutVector &u2)
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition hdf5.h:345
constexpr char O
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr char T
constexpr types::blas_int zero
constexpr char A
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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)
VectorType::value_type * end(VectorType &V)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
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={})
Number compute_mean_value(const hp::MappingCollection< dim, spacedim > &mapping_collection, const DoFHandler< dim, spacedim > &dof, const hp::QCollection< dim > &q_collection, const ReadVector< Number > &v, const unsigned int component)
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14889
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)
STL namespace.
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation