Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
step-8.h
Go to the documentation of this file.
1 false);
771 *   sparsity_pattern.copy_from(dsp);
772 *  
773 *   system_matrix.reinit(sparsity_pattern);
774 *   }
775 *  
776 *  
777 * @endcode
778 *
779 *
780 * <a name="step_8-ElasticProblemassemble_system"></a>
781 * <h4>ElasticProblem::assemble_system</h4>
782 *
783
784 *
785 * The big changes in this program are in the creation of matrix and right
786 * hand side, since they are problem-dependent. We will go through that
787 * process step-by-step, since it is a bit more complicated than in previous
788 * examples.
789 *
790
791 *
792 * The first parts of this function are the same as before, however: setting
793 * up a suitable quadrature formula, initializing an FEValues object for the
794 * (vector-valued) finite element we use as well as the quadrature object,
795 * and declaring a number of auxiliary arrays. In addition, we declare the
796 * ever same two abbreviations: <code>n_q_points</code> and
797 * <code>dofs_per_cell</code>. The number of degrees of freedom per cell we
798 * now obviously ask from the composed finite element rather than from the
799 * underlying scalar Q1 element. Here, it is <code>dim</code> times the
800 * number of degrees of freedom per cell of the Q1 element, though this is
801 * not explicit knowledge we need to care about:
802 *
803 * @code
804 *   template <int dim>
805 *   void ElasticProblem<dim>::assemble_system()
806 *   {
807 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
808 *  
809 *   FEValues<dim> fe_values(fe,
810 *   quadrature_formula,
813 *  
814 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
815 *   const unsigned int n_q_points = quadrature_formula.size();
816 *  
817 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
818 *   Vector<double> cell_rhs(dofs_per_cell);
819 *  
820 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
821 *  
822 * @endcode
823 *
824 * As was shown in previous examples as well, we need a place where to
825 * store the values of the coefficients at all the quadrature points on a
826 * cell. In the present situation, we have two coefficients, lambda and
827 * mu.
828 *
829 * @code
830 *   std::vector<double> lambda_values(n_q_points);
831 *   std::vector<double> mu_values(n_q_points);
832 *  
833 * @endcode
834 *
835 * Well, we could as well have omitted the above two arrays since we will
836 * use constant coefficients for both lambda and mu, which can be declared
837 * like this. They both represent functions always returning the constant
838 * value 1.0. Although we could omit the respective factors in the
839 * assemblage of the matrix, we use them here for purpose of
840 * demonstration.
841 *
842 * @code
844 *  
845 * @endcode
846 *
847 * Like the two constant functions above, we will call the function
848 * right_hand_side just once per cell to make things simpler.
849 *
850 * @code
851 *   std::vector<Tensor<1, dim>> rhs_values(n_q_points);
852 *  
853 * @endcode
854 *
855 * Now we can begin with the loop over all cells:
856 *
857 * @code
858 *   for (const auto &cell : dof_handler.active_cell_iterators())
859 *   {
860 *   fe_values.reinit(cell);
861 *  
862 *   cell_matrix = 0;
863 *   cell_rhs = 0;
864 *  
865 * @endcode
866 *
867 * Next we get the values of the coefficients at the quadrature
868 * points. Likewise for the right hand side:
869 *
870 * @code
871 *   lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
872 *   mu.value_list(fe_values.get_quadrature_points(), mu_values);
873 *   right_hand_side(fe_values.get_quadrature_points(), rhs_values);
874 *  
875 * @endcode
876 *
877 * Then assemble the entries of the local @ref GlossStiffnessMatrix "stiffness matrix" and right
878 * hand side vector. This follows almost one-to-one the pattern
879 * described in the introduction of this example. One of the few
880 * comments in place is that we can compute the number
881 * <code>comp(i)</code>, i.e. the index of the only nonzero vector
882 * component of shape function <code>i</code> using the
883 * <code>fe.system_to_component_index(i).first</code> function call
884 * below.
885 *
886
887 *
888 * (By accessing the <code>first</code> variable of the return value
889 * of the <code>system_to_component_index</code> function, you might
890 * already have guessed that there is more in it. In fact, the
891 * function returns a <code>std::pair@<unsigned int, unsigned
892 * int@></code>, of which the first element is <code>comp(i)</code>
893 * and the second is the value <code>base(i)</code> also noted in the
894 * introduction, i.e. the index of this shape function within all the
895 * shape functions that are nonzero in this component,
896 * i.e. <code>base(i)</code> in the diction of the introduction. This
897 * is not a number that we are usually interested in, however.)
898 *
899
900 *
901 * With this knowledge, we can assemble the local matrix
902 * contributions:
903 *
904 * @code
905 *   for (const unsigned int i : fe_values.dof_indices())
906 *   {
907 *   const unsigned int component_i =
908 *   fe.system_to_component_index(i).first;
909 *  
910 *   for (const unsigned int j : fe_values.dof_indices())
911 *   {
912 *   const unsigned int component_j =
913 *   fe.system_to_component_index(j).first;
914 *  
915 *   for (const unsigned int q_point :
916 *   fe_values.quadrature_point_indices())
917 *   {
918 *   cell_matrix(i, j) +=
919 * @endcode
920 *
921 * The first term is @f$(\lambda \partial_i u_i, \partial_j
922 * v_j) + (\mu \partial_i u_j, \partial_j v_i)@f$. Note
923 * that <code>shape_grad(i,q_point)</code> returns the
924 * gradient of the only nonzero component of the i-th
925 * shape function at quadrature point q_point. The
926 * component <code>comp(i)</code> of the gradient, which
927 * is the derivative of this only nonzero vector
928 * component of the i-th shape function with respect to
929 * the comp(i)th coordinate is accessed by the appended
930 * brackets.
931 *
932 * @code
933 *   (
934 *   (fe_values.shape_grad(i, q_point)[component_i] *
935 *   fe_values.shape_grad(j, q_point)[component_j] *
936 *   lambda_values[q_point])
937 *   +
938 *   (fe_values.shape_grad(i, q_point)[component_j] *
939 *   fe_values.shape_grad(j, q_point)[component_i] *
940 *   mu_values[q_point])
941 *   +
942 * @endcode
943 *
944 * The second term is @f$(\mu \nabla u_i, \nabla
945 * v_j)@f$. We need not access a specific component of
946 * the gradient, since we only have to compute the
947 * scalar product of the two gradients, of which an
948 * overloaded version of <tt>operator*</tt> takes
949 * care, as in previous examples.
950 *
951
952 *
953 * Note that by using the <tt>?:</tt> operator, we only
954 * do this if <tt>component_i</tt> equals
955 * <tt>component_j</tt>, otherwise a zero is added
956 * (which will be optimized away by the compiler).
957 *
958 * @code
959 *   ((component_i == component_j) ?
960 *   (fe_values.shape_grad(i, q_point) *
961 *   fe_values.shape_grad(j, q_point) *
962 *   mu_values[q_point]) :
963 *   0)
964 *   ) *
965 *   fe_values.JxW(q_point);
966 *   }
967 *   }
968 *   }
969 *  
970 * @endcode
971 *
972 * Assembling the right hand side is also just as discussed in the
973 * introduction:
974 *
975 * @code
976 *   for (const unsigned int i : fe_values.dof_indices())
977 *   {
978 *   const unsigned int component_i =
979 *   fe.system_to_component_index(i).first;
980 *  
981 *   for (const unsigned int q_point :
982 *   fe_values.quadrature_point_indices())
983 *   cell_rhs(i) += fe_values.shape_value(i, q_point) *
984 *   rhs_values[q_point][component_i] *
985 *   fe_values.JxW(q_point);
986 *   }
987 *  
988 * @endcode
989 *
990 * The transfer from local degrees of freedom into the global matrix
991 * and right hand side vector does not depend on the equation under
992 * consideration, and is thus the same as in all previous
993 * examples.
994 *
995 * @code
996 *   cell->get_dof_indices(local_dof_indices);
997 *   constraints.distribute_local_to_global(
998 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
999 *   }
1000 *   }
1001 *  
1002 *  
1003 *  
1004 * @endcode
1005 *
1006 *
1007 * <a name="step_8-ElasticProblemsolve"></a>
1008 * <h4>ElasticProblem::solve</h4>
1009 *
1010
1011 *
1012 * The solver does not care about where the system of equations comes from, as
1013 * long as it is positive definite and symmetric (which are the
1014 * requirements for the use of the CG solver), which the system indeed
1015 * is. Therefore, we need not change anything.
1016 *
1017 * @code
1018 *   template <int dim>
1019 *   void ElasticProblem<dim>::solve()
1020 *   {
1021 *   SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
1022 *   SolverCG<Vector<double>> cg(solver_control);
1023 *  
1024 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1025 *   preconditioner.initialize(system_matrix, 1.2);
1026 *  
1027 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
1028 *  
1029 *   constraints.distribute(solution);
1030 *   }
1031 *  
1032 *  
1033 * @endcode
1034 *
1035 *
1036 * <a name="step_8-ElasticProblemrefine_grid"></a>
1037 * <h4>ElasticProblem::refine_grid</h4>
1038 *
1039
1040 *
1041 * The function that does the refinement of the grid is the same as in the
1042 * @ref step_6 "step-6" example. The quadrature formula is adapted to the linear elements
1043 * again. Note that the error estimator by default adds up the estimated
1044 * obtained from all components of the finite element solution, i.e., it
1045 * uses the displacement in all directions with the same weight. If we would
1046 * like the grid to be adapted to the x-displacement only, we could pass the
1047 * function an additional parameter which tells it to do so and do not
1048 * consider the displacements in all other directions for the error
1049 * indicators. However, for the current problem, it seems appropriate to
1050 * consider all displacement components with equal weight.
1051 *
1052 * @code
1053 *   template <int dim>
1054 *   void ElasticProblem<dim>::refine_grid()
1055 *   {
1056 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1057 *  
1058 *   KellyErrorEstimator<dim>::estimate(dof_handler,
1059 *   QGauss<dim - 1>(fe.degree + 1),
1060 *   {},
1061 *   solution,
1062 *   estimated_error_per_cell);
1063 *  
1065 *   estimated_error_per_cell,
1066 *   0.3,
1067 *   0.03);
1068 *  
1069 *   triangulation.execute_coarsening_and_refinement();
1070 *   }
1071 *  
1072 *  
1073 * @endcode
1074 *
1075 *
1076 * <a name="step_8-ElasticProblemoutput_results"></a>
1077 * <h4>ElasticProblem::output_results</h4>
1078 *
1079
1080 *
1081 * The output happens mostly as has been shown in previous examples
1082 * already. The only difference is that the solution function is vector
1083 * valued. The DataOut class takes care of this automatically, but we have
1084 * to give each component of the solution vector a different name.
1085 *
1086
1087 *
1088 * To do this, the DataOut::add_vector() function wants a vector of
1089 * strings. Since the number of components is the same as the number
1090 * of dimensions we are working in, we use the <code>switch</code>
1091 * statement below.
1092 *
1093
1094 *
1095 * We note that some graphics programs have restriction on what
1096 * characters are allowed in the names of variables. deal.II therefore
1097 * supports only the minimal subset of these characters that is supported
1098 * by all programs. Basically, these are letters, numbers, underscores,
1099 * and some other characters, but in particular no whitespace and
1100 * minus/hyphen. The library will throw an exception otherwise, at least
1101 * if in debug mode.
1102 *
1103
1104 *
1105 * After listing the 1d, 2d, and 3d case, it is good style to let the
1106 * program die if we run into a case which we did not consider. You have
1107 * previously already seen the use of the `Assert` macro that generates
1108 * aborts the program with an error message if a condition is not satisfied
1109 * (see @ref step_5 "step-5", for example). We could use this in the `default` case
1110 * below, in the form `Assert(false, ExcNotImplemented())` -- in other words,
1111 * the "condition" here is always `false`, and so the assertion always fails
1112 * and always aborts the program whenever it gets to the default statement.
1113 * This is perhaps more difficult to read than necessary, and consequently
1114 * there is a short-cut: `DEAL_II_NOT_IMPLEMENTED()`. It does the same
1115 * as the form above (with the minor difference that it also aborts the
1116 * program in release mode). It is written in all-caps because that makes
1117 * it stand out visually (and also because it is not actually a function,
1118 * but a macro).
1119 *
1120 * @code
1121 *   template <int dim>
1122 *   void ElasticProblem<dim>::output_results(const unsigned int cycle) const
1123 *   {
1124 *   DataOut<dim> data_out;
1125 *   data_out.attach_dof_handler(dof_handler);
1126 *  
1127 *   std::vector<std::string> solution_names;
1128 *   switch (dim)
1129 *   {
1130 *   case 1:
1131 *   solution_names.emplace_back("displacement");
1132 *   break;
1133 *   case 2:
1134 *   solution_names.emplace_back("x_displacement");
1135 *   solution_names.emplace_back("y_displacement");
1136 *   break;
1137 *   case 3:
1138 *   solution_names.emplace_back("x_displacement");
1139 *   solution_names.emplace_back("y_displacement");
1140 *   solution_names.emplace_back("z_displacement");
1141 *   break;
1142 *   default:
1144 *   }
1145 *  
1146 * @endcode
1147 *
1148 * After setting up the names for the different components of the
1149 * solution vector, we can add the solution vector to the list of
1150 * data vectors scheduled for output. Note that the following
1151 * function takes a vector of strings as second argument, whereas
1152 * the one which we have used in all previous examples took a
1153 * single string there (which was the right choice because
1154 * we had only a single solution variable in all previous examples).
1155 *
1156 * @code
1157 *   data_out.add_data_vector(solution, solution_names);
1158 *   data_out.build_patches();
1159 *  
1160 *   std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
1161 *   data_out.write_vtk(output);
1162 *   }
1163 *  
1164 *  
1165 *  
1166 * @endcode
1167 *
1168 *
1169 * <a name="step_8-ElasticProblemrun"></a>
1170 * <h4>ElasticProblem::run</h4>
1171 *
1172
1173 *
1174 * The <code>run</code> function does the same things as in @ref step_6 "step-6", for
1175 * example. This time, we use the square [-1,1]^d as domain, and we refine
1176 * it globally four times before starting the first iteration.
1177 *
1178
1179 *
1180 * The reason for refining is a bit accidental: we use the QGauss
1181 * quadrature formula with two points in each direction for integration of the
1182 * right hand side; that means that there are four quadrature points on each
1183 * cell (in 2d). If we only refine the initial grid once globally, then there
1184 * will be only four quadrature points in each direction on the
1185 * domain. However, the right hand side function was chosen to be rather
1186 * localized and in that case, by pure chance, it happens that all quadrature
1187 * points lie at points where the right hand side function is zero (in
1188 * mathematical terms, the quadrature points happen to be at points outside
1189 * the <i>support</i> of the right hand side function). The right hand side
1190 * vector computed with quadrature will then contain only zeroes (even though
1191 * it would of course be nonzero if we had computed the right hand side vector
1192 * exactly using the integral) and the solution of the system of
1193 * equations is the zero vector, i.e., a finite element function that is zero
1194 * everywhere. In a sense, we
1195 * should not be surprised that this is happening since we have chosen
1196 * an initial grid that is totally unsuitable for the problem at hand.
1197 *
1198
1199 *
1200 * The unfortunate thing is that if the discrete solution is constant, then
1201 * the error indicators computed by the KellyErrorEstimator class are zero
1202 * for each cell as well, and the call to
1203 * Triangulation::refine_and_coarsen_fixed_number() will not flag any cells
1204 * for refinement (why should it if the indicated error is zero for each
1205 * cell?). The grid in the next iteration will therefore consist of four
1206 * cells only as well, and the same problem occurs again.
1207 *
1208
1209 *
1210 * The conclusion needs to be: while of course we will not choose the
1211 * initial grid to be well-suited for the accurate solution of the problem,
1212 * we must at least choose it such that it has the chance to capture the
1213 * important features of the solution. In this case, it needs to be able to
1214 * see the right hand side. Thus, we refine globally four times. (Any larger
1215 * number of global refinement steps would of course also work.)
1216 *
1217 * @code
1218 *   template <int dim>
1219 *   void ElasticProblem<dim>::run()
1220 *   {
1221 *   for (unsigned int cycle = 0; cycle < 8; ++cycle)
1222 *   {
1223 *   std::cout << "Cycle " << cycle << ':' << std::endl;
1224 *  
1225 *   if (cycle == 0)
1226 *   {
1228 *   triangulation.refine_global(4);
1229 *   }
1230 *   else
1231 *   refine_grid();
1232 *  
1233 *   std::cout << " Number of active cells: "
1234 *   << triangulation.n_active_cells() << std::endl;
1235 *  
1236 *   setup_system();
1237 *  
1238 *   std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
1239 *   << std::endl;
1240 *  
1241 *   assemble_system();
1242 *   solve();
1243 *   output_results(cycle);
1244 *   }
1245 *   }
1246 *   } // namespace Step8
1247 *  
1248 * @endcode
1249 *
1250 *
1251 * <a name="step_8-Thecodemaincodefunction"></a>
1252 * <h3>The <code>main</code> function</h3>
1253 *
1254
1255 *
1256 * After closing the <code>Step8</code> namespace in the last line above, the
1257 * following is the main function of the program and is again exactly like in
1258 * @ref step_6 "step-6" (apart from the changed class names, of course).
1259 *
1260 * @code
1261 *   int main()
1262 *   {
1263 *   try
1264 *   {
1265 *   Step8::ElasticProblem<2> elastic_problem_2d;
1266 *   elastic_problem_2d.run();
1267 *   }
1268 *   catch (std::exception &exc)
1269 *   {
1270 *   std::cerr << std::endl
1271 *   << std::endl
1272 *   << "----------------------------------------------------"
1273 *   << std::endl;
1274 *   std::cerr << "Exception on processing: " << std::endl
1275 *   << exc.what() << std::endl
1276 *   << "Aborting!" << std::endl
1277 *   << "----------------------------------------------------"
1278 *   << std::endl;
1279 *  
1280 *   return 1;
1281 *   }
1282 *   catch (...)
1283 *   {
1284 *   std::cerr << std::endl
1285 *   << std::endl
1286 *   << "----------------------------------------------------"
1287 *   << std::endl;
1288 *   std::cerr << "Unknown exception!" << std::endl
1289 *   << "Aborting!" << std::endl
1290 *   << "----------------------------------------------------"
1291 *   << std::endl;
1292 *   return 1;
1293 *   }
1294 *  
1295 *   return 0;
1296 *   }
1297 * @endcode
1298<a name="step_8-Results"></a><h1>Results</h1>
1299
1300
1301
1302There is not much to be said about the results of this program, other than
1303that they look nice. All images were made using VisIt from the
1304output files that the program wrote to disk. The first two pictures show
1305the @f$x@f$- and @f$y@f$-displacements as scalar components:
1306
1307<table width="100%">
1308<tr>
1309<td>
1310<img src="https://www.dealii.org/images/steps/developer/step-8.x.png" alt="">
1311</td>
1312<td>
1313<img src="https://www.dealii.org/images/steps/developer/step-8.y.png" alt="">
1314</td>
1315</tr>
1316</table>
1317
1318
1319You can clearly see the sources of @f$x@f$-displacement around @f$x=0.5@f$ and
1320@f$x=-0.5@f$, and of @f$y@f$-displacement at the origin.
1321
1322
1323<a name="step-8-extensions"></a>
1324<a name="step_8-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1325
1326
1327What one frequently would like to do is to show the displacement as a vector
1328field, i.e., vectors that for each point illustrate the direction and magnitude
1329of displacement. Unfortunately, that's a bit more involved. To understand why
1330this is so, remember that we have defined our finite element as just a
1331collection of `dim` individual components when we said `fe(FE_Q<dim>(1) ^ dim)`
1332in the constructor. From the perspective of deal.II, these components could
1333just be a pressure and a concentration (i.e., two scalar and unrelated
1334quantities). What we *want* to express somehow is that the components of the finite
1335element actually are the parts of a *vector-valued* quantity, namely the
1336displacement. Absent having specific knowledge about this, the
1337DataOut class assumes that all individual variables we output are separate
1338scalars, and VisIt and Paraview then faithfully assume that this is indeed what it is. In
1339other words, once we have written the data as scalars, there is nothing in
1340these programs that allows us to paste these two scalar fields back together as a
1341vector field. (This is not entirely true: It is possible to *define* a
1342vector-valued field in both of these programs whose components are scalar fields
1343that already exist in an output field; but it is not easy to find how
1344exactly one would do that.)
1345
1346Where we would have to attack this problem is at its root,
1347namely in <code>ElasticProblem::output_results()</code>. The @ref step_22 "step-22" program
1348discussed in some detail how to do that for a more general situation. That said,
1349the following small variation of the `output_results()`, in which we essentially
1350only add the `data_component_interpretation` variable and use the same
1351string `"displacement"` `dim` times, will do the trick:
1352@code
1353 template <int dim>
1354 void ElasticProblem<dim>::output_results(const unsigned int cycle) const
1355 {
1356 DataOut<dim> data_out;
1357 data_out.attach_dof_handler(dof_handler);
1358
1359 const std::vector<std::string> solution_names(dim, "displacement");
1360
1361 const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1362 data_component_interpretation(
1363 dim, DataComponentInterpretation::component_is_part_of_vector);
1364
1365 data_out.add_data_vector(solution,
1366 solution_names,
1367 DataOut<dim>::type_dof_data,
1368 data_component_interpretation);
1369 data_out.build_patches();
1370
1371 std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
1372 data_out.write_vtk(output);
1373 }
1374@endcode
1375As mentioned, @ref step_22 "step-22" discusses what exactly this does. In any case,
1376the displacement vector field that is then output then looks like this
1377(VisIt and Paraview randomly select a few
1378hundred vertices from which to draw the vectors; drawing them from each
1379individual vertex would make the picture unreadable):
1380
1381<img src="https://www.dealii.org/images/steps/developer/step-8.vectors.png" alt="">
1382
1383We note that one may have intuitively expected the
1384solution to be symmetric about the @f$x@f$- and @f$y@f$-axes since the @f$x@f$- and
1385@f$y@f$-forces are symmetric with respect to these axes. However, the force
1386considered as a vector is not symmetric and consequently neither is
1387the solution.
1388 *
1389 *
1390<a name="step_8-PlainProg"></a>
1391<h1> The plain program</h1>
1392@include "step-8.cc"
1393*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
#define Assert(cond, exc)
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.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Definition mpi.cc:746
const Event initial
Definition event.cc:68
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)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
constexpr types::blas_int zero
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 > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
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)
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)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation