661 *
const unsigned int component = 0)
const override
664 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
672 *
class InitialValuesV :
public Function<dim>
676 *
const unsigned int component = 0) const override
688 * Secondly, we have the right hand side forcing term. Boring as we are, we
689 * choose zero here as well:
693 *
class RightHandSide :
public Function<dim>
697 *
const unsigned int component = 0) const override
709 * Finally, we have boundary
values for @f$u@f$ and @f$v@f$. They are as described
710 * in the introduction, one being the time derivative of the other:
714 *
class BoundaryValuesU :
public Function<dim>
718 *
const unsigned int component = 0) const override
723 *
if ((this->
get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
734 *
class BoundaryValuesV :
public Function<dim>
738 *
const unsigned int component = 0) const override
743 *
if ((this->
get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
756 * <a name=
"step_23-ImplementationofthecodeWaveEquationcodeclass"></a>
757 * <h3>Implementation of the <code>WaveEquation</code>
class</h3>
761 * The implementation of the actual logic is actually fairly short, since we
762 * relegate things like assembling the matrices and right hand side vectors
763 * to the library. The rest boils down to not much more than 130 lines of
764 * actual code, a significant fraction of which is boilerplate code that can
765 * be taken from previous example programs (
e.g. the functions that solve
766 * linear systems, or that generate output).
770 * Let
's start with the constructor (for an explanation of the choice of
771 * time step, see the section on Courant, Friedrichs, and Lewy in the
776 * WaveEquation<dim>::WaveEquation()
778 * , dof_handler(triangulation)
779 * , time_step(1. / 64)
781 * , timestep_number(1)
789 * <a name="step_23-WaveEquationsetup_system"></a>
790 * <h4>WaveEquation::setup_system</h4>
794 * The next function is the one that sets up the mesh, DoFHandler, and
795 * matrices and vectors at the beginning of the program, i.e. before the
796 * first time step. The first few lines are pretty much standard if you've
797 * read through the tutorial programs at least up to @ref step_6
"step-6":
801 *
void WaveEquation<dim>::setup_system()
806 * std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
809 * dof_handler.distribute_dofs(fe);
811 * std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
817 * sparsity_pattern.copy_from(dsp);
821 * Then comes a block where we have to initialize the 3 matrices we need
822 * in the course of the program: the mass
matrix, the Laplace
matrix, and
823 * the
matrix @f$M+k^2\theta^2A@f$ used when solving
for @f$U^n@f$ in each time
828 * When setting up these matrices, note that they all make use of the same
829 * sparsity pattern
object. Finally, the reason why matrices and sparsity
830 * patterns are separate objects in deal.II (unlike in many other finite
831 * element or linear algebra classes) becomes clear: in a significant
832 * fraction of applications, one has to hold several matrices that happen
833 * to have the same sparsity pattern, and there is no reason
for them not
834 * to share
this information, rather than re-building and wasting memory
835 * on it several times.
839 * After initializing all of these matrices, we call library
functions
840 * that build the Laplace and mass matrices. All they need is a
DoFHandler
841 *
object and a quadrature formula
object that is to be used
for numerical
842 * integration. Note that in many respects these
functions are better than
843 * what we would usually
do in application programs,
for example because
844 * they automatically parallelize building the matrices
if multiple
845 * processors are available in a machine:
for more information see the
847 * @ref threads
"Parallel computing with multiple processors"
848 * topic. The matrices
for solving linear systems will be filled in the
849 *
run() method because we need to re-apply boundary conditions every time
853 * mass_matrix.reinit(sparsity_pattern);
854 * laplace_matrix.reinit(sparsity_pattern);
855 * matrix_u.reinit(sparsity_pattern);
856 * matrix_v.reinit(sparsity_pattern);
859 *
QGauss<dim>(fe.degree + 1),
862 *
QGauss<dim>(fe.degree + 1),
867 * The rest of the function is spent on setting vector sizes to the
868 * correct value. The final line closes the hanging node constraints
869 *
object. Since we work on a uniformly refined mesh, no constraints exist
870 * or have been computed (i.e. there was no need to call
871 *
DoFTools::make_hanging_node_constraints as in other programs), but we
872 * need a constraints
object in one place further down below anyway.
875 * solution_u.reinit(dof_handler.n_dofs());
876 * solution_v.reinit(dof_handler.n_dofs());
877 * old_solution_u.reinit(dof_handler.n_dofs());
878 * old_solution_v.reinit(dof_handler.n_dofs());
879 * system_rhs.reinit(dof_handler.n_dofs());
881 * constraints.close();
889 * <a name="step_23-WaveEquationsolve_uandWaveEquationsolve_v"></a>
890 * <h4>WaveEquation::solve_u and WaveEquation::solve_v</h4>
894 * The next two functions deal with solving the linear systems associated
895 * with the equations for @f$U^n@f$ and @f$V^n@f$. Both are not particularly
896 * interesting as they pretty much follow the scheme used in all the
897 * previous tutorial programs.
901 * One can make little experiments with preconditioners for the two matrices
902 * we have to
invert. As it turns out, however, for the matrices at hand
903 * here, using Jacobi or SSOR preconditioners reduces the number of
904 * iterations necessary to solve the linear system slightly, but due to the
905 * cost of applying the preconditioner it is no win in terms of run-time. It
906 * is not much of a loss either, but let's keep it simple and just do
911 *
void WaveEquation<dim>::solve_u()
913 *
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
918 * std::cout <<
" u-equation: " << solver_control.last_step()
919 * <<
" CG iterations." << std::endl;
925 *
void WaveEquation<dim>::solve_v()
927 *
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
932 * std::cout <<
" v-equation: " << solver_control.last_step()
933 * <<
" CG iterations." << std::endl;
941 * <a name=
"step_23-WaveEquationoutput_results"></a>
942 * <h4>WaveEquation::output_results</h4>
946 * Likewise, the following function is pretty much what we
've done
947 * before. The only thing worth mentioning is how here we generate a string
948 * representation of the time step number padded with leading zeros to 3
949 * character length using the Utilities::int_to_string function's
second
954 *
void WaveEquation<dim>::output_results() const
959 * data_out.add_data_vector(solution_u,
"U");
960 * data_out.add_data_vector(solution_v,
"V");
962 * data_out.build_patches();
964 *
const std::string filename =
968 * Like @ref step_15
"step-15", since we write output at every time step (and the system
969 * we have to solve is relatively easy), we instruct
DataOut to use the
970 * zlib compression algorithm that is optimized
for speed instead of disk
971 * usage since otherwise plotting the output becomes a bottleneck:
976 * data_out.set_flags(vtk_flags);
977 * std::ofstream output(filename);
978 * data_out.write_vtu(output);
986 * <a name=
"step_23-WaveEquationrun"></a>
987 * <h4>WaveEquation::run</h4>
991 * The following is really the only interesting function of the program. It
992 * contains the
loop over all time steps, but before we get to that we have
993 * to set up the grid,
DoFHandler, and matrices. In addition, we have to
996 * continuous function and computes the @f$L^2@f$ projection of
this function
997 * onto the finite element space described by the
DoFHandler object. Can
't
998 * be any simpler than that:
1001 * template <int dim>
1002 * void WaveEquation<dim>::run()
1006 * VectorTools::project(dof_handler,
1008 * QGauss<dim>(fe.degree + 1),
1009 * InitialValuesU<dim>(),
1011 * VectorTools::project(dof_handler,
1013 * QGauss<dim>(fe.degree + 1),
1014 * InitialValuesV<dim>(),
1019 * The next thing is to loop over all the time steps until we reach the
1020 * end time (@f$T=5@f$ in this case). In each time step, we first have to
1021 * solve for @f$U^n@f$, using the equation @f$(M^n + k^2\theta^2 A^n)U^n =@f$
1022 * @f$(M^{n,n-1} - k^2\theta(1-\theta) A^{n,n-1})U^{n-1} + kM^{n,n-1}V^{n-1}
1023 * +@f$ @f$k\theta \left[k \theta F^n + k(1-\theta) F^{n-1} \right]@f$. Note
1024 * that we use the same mesh for all time steps, so that @f$M^n=M^{n,n-1}=M@f$
1025 * and @f$A^n=A^{n,n-1}=A@f$. What we therefore have to do first is to add up
1026 * @f$MU^{n-1} - k^2\theta(1-\theta) AU^{n-1} + kMV^{n-1}@f$ and the forcing
1027 * terms, and put the result into the <code>system_rhs</code> vector. (For
1028 * these additions, we need a temporary vector that we declare before the
1029 * loop to avoid repeated memory allocations in each time step.)
1033 * The one thing to realize here is how we communicate the time variable
1034 * to the object describing the right hand side: each object derived from
1035 * the Function class has a time field that can be set using the
1036 * Function::set_time and read by Function::get_time. In essence, using
1037 * this mechanism, all functions of space and time are therefore
1038 * considered functions of space evaluated at a particular time. This
1039 * matches well what we typically need in finite element programs, where
1040 * we almost always work on a single time step at a time, and where it
1041 * never happens that, for example, one would like to evaluate a
1042 * space-time function for all times at any given spatial location.
1045 * Vector<double> tmp(solution_u.size());
1046 * Vector<double> forcing_terms(solution_u.size());
1048 * for (; time <= 5; time += time_step, ++timestep_number)
1050 * std::cout << "Time step " << timestep_number << " at t=" << time
1053 * mass_matrix.vmult(system_rhs, old_solution_u);
1055 * mass_matrix.vmult(tmp, old_solution_v);
1056 * system_rhs.add(time_step, tmp);
1058 * laplace_matrix.vmult(tmp, old_solution_u);
1059 * system_rhs.add(-theta * (1 - theta) * time_step * time_step, tmp);
1061 * RightHandSide<dim> rhs_function;
1062 * rhs_function.set_time(time);
1063 * VectorTools::create_right_hand_side(dof_handler,
1064 * QGauss<dim>(fe.degree + 1),
1067 * forcing_terms = tmp;
1068 * forcing_terms *= theta * time_step;
1070 * rhs_function.set_time(time - time_step);
1071 * VectorTools::create_right_hand_side(dof_handler,
1072 * QGauss<dim>(fe.degree + 1),
1076 * forcing_terms.add((1 - theta) * time_step, tmp);
1078 * system_rhs.add(theta * time_step, forcing_terms);
1082 * After so constructing the right hand side vector of the first
1083 * equation, all we have to do is apply the correct boundary
1084 * values. As for the right hand side, this is a space-time function
1085 * evaluated at a particular time, which we interpolate at boundary
1086 * nodes and then use the result to apply boundary values as we
1087 * usually do. The result is then handed off to the solve_u()
1092 * BoundaryValuesU<dim> boundary_values_u_function;
1093 * boundary_values_u_function.set_time(time);
1095 * std::map<types::global_dof_index, double> boundary_values;
1096 * VectorTools::interpolate_boundary_values(dof_handler,
1098 * boundary_values_u_function,
1103 * The matrix for solve_u() is the same in every time steps, so one
1104 * could think that it is enough to do this only once at the
1105 * beginning of the simulation. However, since we need to apply
1106 * boundary values to the linear system (which eliminate some matrix
1107 * rows and columns and give contributions to the right hand side),
1108 * we have to refill the matrix in every time steps before we
1109 * actually apply boundary data. The actual content is very simple:
1110 * it is the sum of the mass matrix and a weighted Laplace matrix:
1113 * matrix_u.copy_from(mass_matrix);
1114 * matrix_u.add(theta * theta * time_step * time_step, laplace_matrix);
1115 * MatrixTools::apply_boundary_values(boundary_values,
1125 * The second step, i.e. solving for @f$V^n@f$, works similarly, except
1126 * that this time the matrix on the left is the mass matrix (which we
1127 * copy again in order to be able to apply boundary conditions, and
1128 * the right hand side is @f$MV^{n-1} - k\left[ \theta A U^n +
1129 * (1-\theta) AU^{n-1}\right]@f$ plus forcing terms. Boundary values
1130 * are applied in the same way as before, except that now we have to
1131 * use the BoundaryValuesV class:
1134 * laplace_matrix.vmult(system_rhs, solution_u);
1135 * system_rhs *= -theta * time_step;
1137 * mass_matrix.vmult(tmp, old_solution_v);
1138 * system_rhs += tmp;
1140 * laplace_matrix.vmult(tmp, old_solution_u);
1141 * system_rhs.add(-time_step * (1 - theta), tmp);
1143 * system_rhs += forcing_terms;
1146 * BoundaryValuesV<dim> boundary_values_v_function;
1147 * boundary_values_v_function.set_time(time);
1149 * std::map<types::global_dof_index, double> boundary_values;
1150 * VectorTools::interpolate_boundary_values(dof_handler,
1152 * boundary_values_v_function,
1154 * matrix_v.copy_from(mass_matrix);
1155 * MatrixTools::apply_boundary_values(boundary_values,
1164 * Finally, after both solution components have been computed, we
1165 * output the result, compute the energy in the solution, and go on to
1166 * the next time step after shifting the present solution into the
1167 * vectors that hold the solution at the previous time step. Note the
1168 * function SparseMatrix::matrix_norm_square that can compute
1169 * @f$\left<V^n,MV^n\right>@f$ and @f$\left<U^n,AU^n\right>@f$ in one step,
1170 * saving us the expense of a temporary vector and several lines of
1176 * std::cout << " Total energy: "
1177 * << (mass_matrix.matrix_norm_square(solution_v) +
1178 * laplace_matrix.matrix_norm_square(solution_u)) /
1182 * old_solution_u = solution_u;
1183 * old_solution_v = solution_v;
1186 * } // namespace Step23
1192 * <a name="step_23-Thecodemaincodefunction"></a>
1193 * <h3>The <code>main</code> function</h3>
1197 * What remains is the main function of the program. There is nothing here
1198 * that hasn't been shown in several of the previous programs:
1205 *
using namespace Step23;
1207 * WaveEquation<2> wave_equation_solver;
1208 * wave_equation_solver.run();
1210 *
catch (std::exception &exc)
1212 * std::cerr << std::endl
1214 * <<
"----------------------------------------------------"
1216 * std::cerr <<
"Exception on processing: " << std::endl
1217 * << exc.what() << std::endl
1218 * <<
"Aborting!" << std::endl
1219 * <<
"----------------------------------------------------"
1226 * std::cerr << std::endl
1228 * <<
"----------------------------------------------------"
1230 * std::cerr <<
"Unknown exception!" << std::endl
1231 * <<
"Aborting!" << std::endl
1232 * <<
"----------------------------------------------------"
1240<a name=
"step_23-Results"></a><h1>Results</h1>
1243When the program is
run, it produces the following output:
1245Number of active cells: 16384
1246Number of degrees of freedom: 16641
1248Time step 1 at t=0.015625
1249 u-equation: 8 CG iterations.
1250 v-equation: 22 CG iterations.
1251 Total energy: 1.17887
1252Time step 2 at t=0.03125
1253 u-equation: 8 CG iterations.
1254 v-equation: 20 CG iterations.
1255 Total energy: 2.9655
1256Time step 3 at t=0.046875
1257 u-equation: 8 CG iterations.
1258 v-equation: 21 CG iterations.
1259 Total energy: 4.33761
1260Time step 4 at t=0.0625
1261 u-equation: 7 CG iterations.
1262 v-equation: 21 CG iterations.
1263 Total energy: 5.35499
1264Time step 5 at t=0.078125
1265 u-equation: 7 CG iterations.
1266 v-equation: 21 CG iterations.
1267 Total energy: 6.18652
1268Time step 6 at t=0.09375
1269 u-equation: 7 CG iterations.
1270 v-equation: 20 CG iterations.
1271 Total energy: 6.6799
1275Time step 31 at t=0.484375
1276 u-equation: 7 CG iterations.
1277 v-equation: 20 CG iterations.
1278 Total energy: 21.9068
1279Time step 32 at t=0.5
1280 u-equation: 7 CG iterations.
1281 v-equation: 20 CG iterations.
1282 Total energy: 23.3394
1283Time step 33 at t=0.515625
1284 u-equation: 7 CG iterations.
1285 v-equation: 20 CG iterations.
1286 Total energy: 23.1019
1290Time step 319 at t=4.98438
1291 u-equation: 7 CG iterations.
1292 v-equation: 20 CG iterations.
1293 Total energy: 23.1019
1295 u-equation: 7 CG iterations.
1296 v-equation: 20 CG iterations.
1297 Total energy: 23.1019
1300What we see immediately is that the energy is a
constant at least after
1301@f$t=\frac 12@f$ (until which the boundary source term @f$g@f$ is nonzero, injecting
1302energy into the system).
1304In addition to the screen output, the program writes the solution of each time
1305step to an output file. If we process them adequately and paste them into a
1306movie, we get the following:
1308<img src=
"https://www.dealii.org/images/steps/developer/step-23.movie.gif" alt=
"Animation of the solution of step 23.">
1310The movie shows the generated wave nice traveling through the domain and back,
1311being reflected at the clamped boundary. Some numerical noise is trailing the
1312wave, an artifact of a too-large mesh size that can be reduced by reducing the
1313mesh width and the time step.
1316<a name=
"step-23-extensions"></a>
1317<a name=
"step_23-Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
1320If you want to explore a bit,
try out some of the following things:
1322 <li>Varying @f$\theta@f$. This gives different time stepping schemes, some of
1323 which are stable
while others are not. Take a look at how the energy
1326 <li>Different
initial and boundary conditions, right hand sides.
1328 <li>More complicated domains or more refined meshes. Remember that the time
1329 step needs to be bounded by the mesh width, so changing the mesh should
1330 always involve also changing the time step. We will come back to
this issue
1331 in @ref step_24
"step-24".
1333 <li>Variable coefficients: In real media, the wave speed is often
1334 variable. In particular, the
"real" wave equation in realistic media would
1337 \rho(x) \frac{\partial^2 u}{\partial t^2}
1342 where @f$\rho(x)@f$ is the density of the material, and @f$a(x)@f$ is related to the
1343 stiffness coefficient. The wave speed is then @f$c=\
sqrt{a/\rho}@f$.
1345 To make such a change, we would have to compute the mass and Laplace
1346 matrices with a variable coefficient. Fortunately,
this isn
't too hard: the
1347 functions MatrixCreator::create_laplace_matrix and
1348 MatrixCreator::create_mass_matrix have additional default parameters that can
1349 be used to pass non-constant coefficient functions to them. The required
1350 changes are therefore relatively small. On the other hand, care must be
1351 taken again to make sure the time step is within the allowed range.
1353 <li>In the in-code comments, we discussed the fact that the matrices for
1354 solving for @f$U^n@f$ and @f$V^n@f$ need to be reset in every time because of
1355 boundary conditions, even though the actual content does not change. It is
1356 possible to avoid copying by not eliminating columns in the linear systems,
1357 which is implemented by appending a @p false argument to the call:
1359 MatrixTools::apply_boundary_values(boundary_values,
1366 <li>deal.II being a library that supports adaptive meshes it would of course be
1367 nice if this program supported change the mesh every few time steps. Given the
1368 structure of the solution — a wave that travels through the domain —
1369 it would seem appropriate if we only refined the mesh where the wave currently is,
1370 and not simply everywhere. It is intuitively clear that we should be able to
1371 save a significant amount of cells this way. (Though upon further thought one
1372 realizes that this is really only the case in the initial stages of the simulation.
1373 After some time, for wave phenomena, the domain is filled with reflections of
1374 the initial wave going in every direction and filling every corner of the domain.
1375 At this point, there is in general little one can gain using local mesh
1378 To make adaptively changing meshes possible, there are basically two routes.
1379 The "correct" way would be to go back to the weak form we get using Rothe's
1380 method. For example, the
first of the two equations to be solved in each time
1381 step looked like
this:
1383 (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
1384 (u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla u^{n-1},\nabla \varphi)
1389 \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
1392 Now, note that we solve
for @f$u^n@f$ on mesh @f${\mathbb T}^n@f$, and
1393 consequently the test
functions @f$\varphi@f$ have to be from the space
1394 @f$V_h^n@f$ as well. As discussed in the introduction, terms like
1395 @f$(u^{n-1},\varphi)@f$ then require us to integrate the solution of the
1396 previous step (which may have been computed on a different mesh
1397 @f${\mathbb T}^{n-1}@f$) against the test functions of the current mesh,
1398 leading to a matrix @f$M^{n,n-1}@f$. This process of integrating shape
1399 functions from different meshes is, at best, awkward. It can be done
1400 but because it is difficult to ensure that @f${\mathbb T}^{n-1}@f$ and
1401 @f${\mathbb T}^{n}@f$ differ by at most one
level of refinement, one
1402 has to recursively match cells from both meshes. It is feasible to
1403 do this, but it leads to lengthy and not entirely obvious code.
1405 The
second approach is the following: whenever we change the mesh,
1406 we simply
interpolate the solution from the last time step on the old
1408 instead of the equation above, we would solve
1410 (u^n,\varphi) + k^2\theta^2(\nabla u^n,\nabla \varphi) &=&
1411 (I^n u^{n-1},\varphi) - k^2\theta(1-\theta)(\nabla I^n u^{n-1},\nabla \varphi)
1413 k(I^n v^{n-1},\varphi)
1416 \theta (f^n,\varphi) + (1-\theta) (f^{n-1},\varphi)
1419 where @f$I^n@f$ interpolates a given function onto mesh @f${\mathbb T}^n@f$.
1420 This is a much simpler approach because, in each time step, we no
1421 longer have to worry whether @f$u^{n-1},v^{n-1}@f$ were computed on the
1422 same mesh as we are
using now or on a different mesh. Consequently,
1423 the only changes to the code necessary are the addition of a function
1424 that computes the error, marks cells
for refinement, sets up a
1426 rebuilds matrices and right hand side vectors on the
new mesh. Neither
1427 the
functions building the matrices and right hand sides, nor the
1428 solvers need to be changed.
1430 While
this second approach is, strictly speaking,
1431 not quite correct in the Rothe framework (it introduces an addition source
1432 of error, namely the interpolation), it is nevertheless what
1433 almost everyone solving time dependent equations does. We will use
this
1434 method in @ref step_31
"step-31",
for example.
1438<a name=
"step_23-PlainProg"></a>
1439<h1> The plain program</h1>
1440@include
"step-23.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
#define Assert(cond, exc)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
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())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
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)
static constexpr double PI
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DataOutBase::CompressionLevel compression_level
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)