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