416 *
vector_value_list(
const std::vector<
Point<dim>> &points,
421 *
for (
unsigned int p = 0; p < points.size(); ++p)
429 * <a name=
"step_29-ThecodeParameterReadercodeclass"></a>
448 *
void declare_parameters();
465 * <a name=
"step_29-codeParameterReaderdeclare_parameterscode"></a>
466 * <
h4><
code>ParameterReader::declare_parameters</
code></
h4>
478 *
void ParameterReader::declare_parameters()
491 *
prm.enter_subsection(
"Mesh & geometry parameters");
493 *
prm.declare_entry(
"Number of refinements",
496 *
"Number of global mesh refinement steps "
497 *
"applied to initial coarse grid");
499 *
prm.declare_entry(
"Focal distance",
502 *
"Distance of the focal point of the lens "
505 *
prm.leave_subsection();
516 *
prm.enter_subsection(
"Physical constants");
522 *
prm.leave_subsection();
532 *
prm.enter_subsection(
"Output parameters");
534 *
prm.declare_entry(
"Output filename",
537 *
"Name of the output file (without extension)");
556 *
template-parameter-
dependent class.)
To find out what parameters
567 *
prm.leave_subsection();
573 * <a name=
"step_29-codeParameterReaderread_parameterscode"></a>
587 *
declare_parameters();
597 * <a name=
"step_29-ThecodeComputeIntensitycodeclass"></a>
668 * compute @
f$|
u|@
f$,
so we're good with the update_values flag.
672 * ComputeIntensity<dim>::ComputeIntensity()
673 * : DataPostprocessorScalar<dim>("Intensity", update_values)
679 * The actual postprocessing happens in the following function. Its input is
680 * an object that stores values of the function (which is here vector-valued)
681 * representing the data vector given to DataOut::add_data_vector, evaluated
682 * at all evaluation points where we generate output, and some tensor objects
718 *
const std::complex<double>
u(
inputs.solution_values[p](0),
719 *
inputs.solution_values[p](1));
729 * <a name=
"step_29-ThecodeUltrasoundProblemcodeclass"></a>
785 *
, fe(
FE_Q<dim>(1) ^ 2)
791 * <a name=
"step_29-codeUltrasoundProblemmake_gridcode"></a>
810 *
std::
cout <<
"Generating grid... ";
820 *
prm.enter_subsection(
"Mesh & geometry parameters");
823 *
const unsigned int n_refinements = prm.get_integer(
"Number of refinements");
825 *
prm.leave_subsection();
863 *
if (face->at_boundary() &&
866 *
face->set_boundary_id(1);
867 *
face->set_manifold_id(1);
895 *
std::cout <<
"done (" << timer.cpu_time() <<
"s)" << std::endl;
897 *
std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
905 * <a name=
"step_29-codeUltrasoundProblemsetup_systemcode"></a>
919 *
std::cout <<
"Setting up system... ";
922 *
dof_handler.distribute_dofs(fe);
926 *
sparsity_pattern.copy_from(
dsp);
928 *
system_matrix.reinit(sparsity_pattern);
930 *
solution.reinit(dof_handler.n_dofs());
933 *
std::cout <<
"done (" << timer.cpu_time() <<
"s)" << std::endl;
935 *
std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
943 * <a name=
"step_29-codeUltrasoundProblemassemble_systemcode"></a>
949 * right
hand side vector:
955 *
std::cout <<
"Assembling system matrix... ";
968 *
prm.enter_subsection(
"Physical constants");
970 *
const double omega = prm.get_double(
"omega"), c = prm.get_double(
"c");
972 *
prm.leave_subsection();
987 *
dofs_per_cell = fe.n_dofs_per_cell();
993 * need the values and gradients of the shape functions, and of course the
994 * quadrature weights. For the terms involving the boundary integrals,
995 * only shape function values and the quadrature weights are necessary.
998 * FEValues<dim> fe_values(fe,
999 * quadrature_formula,
1000 * update_values | update_gradients |
1001 * update_JxW_values);
1003 * FEFaceValues<dim> fe_face_values(fe,
1004 * face_quadrature_formula,
1005 * update_values | update_JxW_values);
1009 * As usual, the system matrix is assembled cell by cell, and we need a
1010 * matrix for storing the local cell contributions as well as an index
1011 * vector to transfer the cell contributions to the appropriate location
1012 * in the global system matrix after.
1015 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1016 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1018 * for (const auto &cell : dof_handler.active_cell_iterators())
1022 * On each cell, we first need to reset the local contribution matrix
1023 * and request the FEValues object to compute the shape functions for
1028 * fe_values.reinit(cell);
1030 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1032 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1036 * At this point, it is important to keep in mind that we are
1037 * dealing with a finite element system with two
1038 * components. Due to the way we constructed this FESystem,
1039 * namely as the Cartesian product of two scalar finite
1040 * element fields, each shape function has only a single
1041 * nonzero component (they are, in deal.II lingo, @ref
1042 * GlossPrimitive "primitive"). Hence, each shape function
1043 * can be viewed as one of the @f$\phi@f$'s
or @
f$\psi@
f$'s from the
1044 * introduction, and similarly the corresponding degrees of
1045 * freedom can be attributed to either @f$\alpha@f$ or @f$\beta@f$.
1046 * As we iterate through all the degrees of freedom on the
1047 * current cell however, they do not come in any particular
1048 * order, and so we cannot decide right away whether the DoFs
1049 * with index @f$i@f$ and @f$j@f$ belong to the real or imaginary part
1050 * of our solution. On the other hand, if you look at the
1051 * form of the system matrix in the introduction, this
1052 * distinction is crucial since it will determine to which
1053 * block in the system matrix the contribution of the current
1054 * pair of DoFs will go and hence which quantity we need to
1055 * compute from the given two shape functions. Fortunately,
1056 * the FESystem object can provide us with this information,
1057 * namely it has a function
1058 * FESystem::system_to_component_index(), that for each local
1059 * DoF index returns a pair of integers of which the first
1060 * indicates to which component of the system the DoF
1061 * belongs. The second integer of the pair indicates which
1062 * index the DoF has in the scalar base finite element field,
1063 * but this information is not relevant here. If you want to
1064 * know more about this function and the underlying scheme
1065 * behind primitive vector valued elements, take a look at
1066 * @ref step_8 "step-8" or the @ref vector_valued topic, where these topics
1067 * are explained in depth.
1070 * if (fe.system_to_component_index(i).first ==
1071 * fe.system_to_component_index(j).first)
1075 * If both DoFs @f$i@f$ and @f$j@f$ belong to same component,
1076 * i.e. their shape functions are both @f$\phi@f$'s
or both
1077 * @
f$\psi@
f$'s, the contribution will end up in one of the
1078 * diagonal blocks in our system matrix, and since the
1079 * corresponding entries are computed by the same formula,
1080 * we do not bother if they actually are @f$\phi@f$ or @f$\psi@f$
1081 * shape functions. We can simply compute the entry by
1082 * iterating over all quadrature points and adding up
1083 * their contributions, where values and gradients of the
1084 * shape functions are supplied by our FEValues object.
1090 * for (unsigned int q_point = 0; q_point < n_q_points;
1092 * cell_matrix(i, j) +=
1093 * (((fe_values.shape_value(i, q_point) *
1094 * fe_values.shape_value(j, q_point)) *
1095 * (-omega * omega) +
1096 * (fe_values.shape_grad(i, q_point) *
1097 * fe_values.shape_grad(j, q_point)) *
1099 * fe_values.JxW(q_point));
1103 * You might think that we would have to specify which
1104 * component of the shape function we'd like to evaluate
1127 *
for (
const auto face_no : cell->face_indices())
1128 *
if (cell->face(
face_no)->at_boundary() &&
1149 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1150 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1151 *
if ((fe.system_to_component_index(i).first !=
1152 *
fe.system_to_component_index(
j).first) &&
1153 *
fe.has_support_on_face(i,
face_no) &&
1154 *
fe.has_support_on_face(
j,
face_no))
1159 * it we would simply add up terms to the local cell
1160 * matrix that happen to be zero because at least one of
1161 * the shape functions happens to be zero. However, we can
1162 * save that work by adding the checks above.
1166 * In either case, these DoFs will contribute to the
1167 * boundary integrals in the off-diagonal blocks of the
1168 * system matrix. To compute the integral, we loop over
1169 * all the quadrature points on the face and sum up the
1170 * contribution weighted with the quadrature weights that
1171 * the face quadrature rule provides. In contrast to the
1172 * entries on the diagonal blocks, here it does matter
1173 * which one of the shape functions is a @f$\psi@f$ and which
1174 * one is a @f$\phi@f$, since that will determine the sign of
1175 * the entry. We account for this by a simple conditional
1176 * statement that determines the correct sign. Since we
1177 * already checked that DoF @f$i@f$ and @f$j@f$ belong to
1178 * different components, it suffices here to test for one
1179 * of them to which component it belongs.
1182 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1184 * cell_matrix(i, j) +=
1185 * ((fe.system_to_component_index(i).first == 0) ? -1 :
1187 * fe_face_values.shape_value(i, q_point) *
1188 * fe_face_values.shape_value(j, q_point) * c * omega *
1189 * fe_face_values.JxW(q_point);
1194 * Now we are done with this cell and have to transfer its
1195 * contributions from the local to the global system matrix. To this
1196 * end, we first get a list of the global indices of the this cells
1200 * cell->get_dof_indices(local_dof_indices);
1205 * ...and then add the entries to the system matrix one by one:
1208 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1209 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1210 * system_matrix.add(local_dof_indices[i],
1211 * local_dof_indices[j],
1212 * cell_matrix(i, j));
1218 * The only thing left are the Dirichlet boundary values on @f$\Gamma_1@f$,
1219 * which is characterized by the boundary indicator 1. The Dirichlet
1220 * values are provided by the <code>DirichletBoundaryValues</code> class
1224 * std::map<types::global_dof_index, double> boundary_values;
1225 * VectorTools::interpolate_boundary_values(dof_handler,
1227 * DirichletBoundaryValues<dim>(),
1230 * MatrixTools::apply_boundary_values(boundary_values,
1236 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1244 * <a name="step_29-codeUltrasoundProblemsolvecode"></a>
1245 * <h4><code>UltrasoundProblem::solve</code></h4>
1249 * As already mentioned in the introduction, the system matrix is neither
1250 * symmetric nor definite, and so it is not quite obvious how to come up
1251 * with an iterative solver and a preconditioner that do a good job on this
1252 * matrix. (For more on this topic, see also the
1253 * <a href="#extensions">Possibilities for extensions</a> section below.)
1254 * We chose instead to go a different way and solve the linear
1255 * system with the sparse LU decomposition provided by UMFPACK. This is
1256 * often a good first choice for 2d problems and works reasonably well even
1257 * for moderately large numbers of DoFs. The deal.II interface to UMFPACK
1258 * is implemented in the SparseDirectUMFPACK class, which is very easy to
1259 * use and allows us to solve our linear system with just 3 lines of code.
1263 * Note again that for compiling this example program, you need to have the
1264 * deal.II library built with UMFPACK support.
1267 * template <int dim>
1268 * void UltrasoundProblem<dim>::solve()
1270 * std::cout << "Solving linear system... ";
1275 * The code to solve the linear system is short: First, we allocate an
1276 * object of the right type. The following call to
1277 * SparseDirectUMFPACK::solve() takes as argument the matrix to decompose,
1278 * and a vector that upon input equals the right hand side of the linear
1279 * system to be solved, and upon output contains the solution of the linear
1280 * system. To satisfy this input/output requirement, we first assign the
1281 * right hand side vector to the `solution` variable.
1284 * SparseDirectUMFPACK A_direct;
1286 * solution = system_rhs;
1287 * A_direct.solve(system_matrix, solution);
1290 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1298 * <a name="step_29-codeUltrasoundProblemoutput_resultscode"></a>
1299 * <h4><code>UltrasoundProblem::output_results</code></h4>
1303 * Here we output our solution @f$v@f$ and @f$w@f$ as well as the derived quantity
1304 * @f$|u|@f$ in the format specified in the parameter file. Most of the work for
1305 * deriving @f$|u|@f$ from @f$v@f$ and @f$w@f$ was already done in the implementation of
1306 * the <code>ComputeIntensity</code> class, so that the output routine is
1307 * rather straightforward and very similar to what is done in the previous
1311 * template <int dim>
1312 * void UltrasoundProblem<dim>::output_results() const
1314 * std::cout << "Generating output... ";
1319 * Define objects of our <code>ComputeIntensity</code> class and a DataOut
1323 * ComputeIntensity<dim> intensities;
1324 * DataOut<dim> data_out;
1326 * data_out.attach_dof_handler(dof_handler);
1330 * Next we query the output-related parameters from the ParameterHandler.
1331 * The DataOut::parse_parameters call acts as a counterpart to the
1332 * DataOutInterface<1>::declare_parameters call in
1333 * <code>ParameterReader::declare_parameters</code>. It collects all the
1334 * output format related parameters from the ParameterHandler and sets the
1335 * corresponding properties of the DataOut object accordingly.
1338 * prm.enter_subsection("Output parameters");
1340 * const std::string output_filename = prm.get("Output filename");
1341 * data_out.parse_parameters(prm);
1343 * prm.leave_subsection();
1347 * Now we put together the filename from the base name provided by the
1348 * ParameterHandler and the suffix which is provided by the DataOut class
1349 * (the default suffix is set to the right type that matches the one set
1350 * in the .prm file through parse_parameters()):
1353 * const std::string filename = output_filename + data_out.default_suffix();
1355 * std::ofstream output(filename);
1359 * The solution vectors @f$v@f$ and @f$w@f$ are added to the DataOut object in the
1363 * std::vector<std::string> solution_names;
1364 * solution_names.emplace_back("Re_u");
1365 * solution_names.emplace_back("Im_u");
1367 * data_out.add_data_vector(solution, solution_names);
1371 * For the intensity, we just call <code>add_data_vector</code> again, but
1372 * this with our <code>ComputeIntensity</code> object as the second
1373 * argument, which effectively adds @f$|u|@f$ to the output data:
1376 * data_out.add_data_vector(solution, intensities);
1380 * The last steps are as before. Note that the actual output format is now
1381 * determined by what is stated in the input file, i.e. one can change the
1382 * output format without having to re-compile this program:
1385 * data_out.build_patches();
1386 * data_out.write(output);
1389 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1397 * <a name="step_29-codeUltrasoundProblemruncode"></a>
1398 * <h4><code>UltrasoundProblem::run</code></h4>
1402 * Here we simply execute our functions one after the other:
1405 * template <int dim>
1406 * void UltrasoundProblem<dim>::run()
1410 * assemble_system();
1414 * } // namespace Step29
1420 * <a name="step_29-Thecodemaincodefunction"></a>
1421 * <h4>The <code>main</code> function</h4>
1425 * Finally the <code>main</code> function of the program. It has the same
1426 * structure as in almost all of the other tutorial programs. The only
1427 * exception is that we define ParameterHandler and
1428 * <code>ParameterReader</code> objects, and let the latter read in the
1429 * parameter values from a textfile called <code>@ref step_29 "step-29".prm</code>. The
1430 * values so read are then handed over to an instance of the UltrasoundProblem
1438 * using namespace dealii;
1439 * using namespace Step29;
1441 * ParameterHandler prm;
1442 * ParameterReader param(prm);
1443 * param.read_parameters("step-29.prm");
1445 * UltrasoundProblem<2> ultrasound_problem(prm);
1446 * ultrasound_problem.run();
1448 * catch (std::exception &exc)
1450 * std::cerr << std::endl
1452 * << "----------------------------------------------------"
1454 * std::cerr << "Exception on processing: " << std::endl
1455 * << exc.what() << std::endl
1456 * << "Aborting!" << std::endl
1457 * << "----------------------------------------------------"
1463 * std::cerr << std::endl
1465 * << "----------------------------------------------------"
1467 * std::cerr << "Unknown exception!" << std::endl
1468 * << "Aborting!" << std::endl
1469 * << "----------------------------------------------------"
1476<a name="step_29-Results"></a><h1>Results</h1>
1479The current program reads its run-time parameters from an input file
1480called <code>step-29.prm</code> that looks like this:
1482subsection Mesh & geometry parameters
1483 # Distance of the focal point of the lens to the x-axis
1484 set Focal distance = 0.3
1486 # Number of global mesh refinement steps applied to initial coarse grid
1487 set Number of refinements = 5
1491subsection Physical constants
1500subsection Output parameters
1501 # Name of the output file (without extension)
1502 set Output file = solution
1504 # A name for the output format to be used
1505 set Output format = vtu
1509As can be seen, we set
1510@f$d=0.3@f$, which amounts to a focus of the transducer lens
1511at @f$x=0.5@f$, @f$y=0.3@f$. The coarse mesh is refined 5 times,
1512resulting in 160x160 cells, and the output is written in vtu
1513format. The parameter reader understands many more parameters
1514pertaining in particular to the generation of output, but we
1515need none of these parameters here and therefore stick with
1516their default values.
1562<table
align=
"center" class=
"doxtable">
1565 <
img src=
"https://www.dealii.org/images/steps/developer/step-29.v.png" alt=
"v = Re(u)">
1568 <
img src=
"https://www.dealii.org/images/steps/developer/step-29.w.png" alt=
"w = Im(u)">
1573 <
img src=
"https://www.dealii.org/images/steps/developer/step-29.intensity.png" alt=
"|u|">
1590<table
align=
"center" class=
"doxtable">
1593 <
img src=
"https://www.dealii.org/images/steps/developer/step-29.surface.png" alt=
"|u|">
1596 <
img src=
"https://www.dealii.org/images/steps/developer/step-29.contours.png" alt=
"|u|">
1609[ 66%] Built target @ref step_29 "step-29"
1610[100%] Run @ref step_29 "step-29" with Release configuration
1611Generating grid... done (0.008767s)
1612 Number of active cells: 25600
1613Setting up system... done (0.007043s)
1614 Number of degrees of freedom: 51842
1615Assembling system matrix... done (0.020942s)
1616Solving linear system... done (0.585853s)
1617Generating output... done (0.092681s)
1618[100%] Built target run
1621[ 66%] Built target @ref step_29 "step-29"
1622[100%] Run @ref step_29 "step-29" with Release configuration
1623Generating grid... done (0.027986s)
1624 Number of active cells: 102400
1625Setting up system... done (0.070022s)
1626 Number of degrees of freedom: 206082
1627Assembling system matrix... done (0.075252s)
1628Solving linear system... done (4.13989s)
1629Generating output... done (0.369556s)
1630[100%] Built target run
1633[ 66%] Built target @ref step_29 "step-29"
1634[100%] Run @ref step_29 "step-29" with Release configuration
1635Generating grid... done (0.186439s)
1636 Number of active cells: 409600
1637Setting up system... done (0.272828s)
1638 Number of degrees of freedom: 821762
1639Assembling system matrix... done (0.316252s)
1640Solving linear system... done (33.1424s)
1641Generating output... done (1.36328s)
1642[100%] Built target run
1645Each time we refine the mesh once more, so the number of cells and degrees
1646of freedom roughly quadruples from each step to the next. As can be seen,
1647generating the grid, setting up degrees of freedom, assembling the
1648linear system, and generating output scale pretty closely to linear,
1649whereas solving the linear system is an operation that requires 8
1650times more time each time the number of degrees of freedom is
1651increased by a factor of 4, i.e. it is @f${\cal O}(N^{3/2})@f$. This can
1652be explained by the fact that (using optimal ordering) the
1653bandwidth of a finite element matrix is @f$B={\cal O}(N^{(dim-1)/dim})@f$,
1654and the effort to solve a banded linear system using LU decomposition
1655is @f${\cal O}(BN)@f$. This also explains why the program does run in 3d
1656as well (after changing the dimension on the
1657<code>UltrasoundProblem</code> object), but scales very badly and
1658takes extraordinary patience before it finishes solving the linear
1659system on a mesh with appreciable resolution, even though all the
1660other parts of the program scale very nicely.
1664<a name="step-29-extensions"></a>
1665<a name="step_29-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1668An obvious possible extension for this program is to run it in 3d
1669— after all, the world around us is three-dimensional, and
1670ultrasound beams propagate in three-dimensional media. You can try
1671this by simply changing the template parameter of the principal class
1672in <code>main()</code> and running it. This won't get
you very far,
1674the parameter file.
You'll simply run out of memory as both the mesh
1675(with its @f$(2^5)^3 \cdot 5^3=2^{15}\cdot 125 \approx 4\cdot 10^6@f$ cells)
1676and in particular the sparse direct solver take too much memory. You
1677can solve with 3 global refinement steps, however, if you have a bit
1678of time: in early 2011, the direct solve takes about half an
1681solution
's waves accurately, and you can see this in plots of the
1682solution. Consequently, this is one of the cases where adaptivity is
1687<a name=
"step_29-PlainProg"></a>
static void declare_parameters(ParameterHandler &prm)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
#define AssertDimension(dim1, dim2)
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)
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
std::vector< index_type > data
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation