146 * #include <deal.II/base/
numbers.h>
147 * #include <deal.II/base/thread_management.h>
148 * #include <deal.II/grid/
tria.h>
149 * #include <deal.II/dofs/dof_handler.h>
150 * #include <deal.II/grid/grid_generator.h>
151 * #include <deal.II/grid/tria_accessor.h>
152 * #include <deal.II/grid/tria_iterator.h>
153 * #include <deal.II/dofs/dof_accessor.h>
154 * #include <deal.II/dofs/dof_renumbering.h>
155 * #include <deal.II/fe/fe_q.h>
156 * #include <deal.II/fe/fe_dgq.h>
157 * #include <deal.II/fe/fe_system.h>
158 * #include <deal.II/dofs/dof_tools.h>
159 * #include <deal.II/fe/fe_values.h>
160 * #include <deal.II/base/quadrature_lib.h>
161 * #include <deal.II/base/function.h>
162 * #include <deal.II/numerics/vector_tools.h>
163 * #include <deal.II/numerics/matrix_tools.h>
164 * #include <deal.II/numerics/derivative_approximation.h>
165 * #include <deal.II/lac/block_vector.h>
166 * #include <deal.II/lac/full_matrix.h>
167 * #include <deal.II/lac/block_sparse_matrix.h>
168 * #include <deal.II/lac/block_sparsity_pattern.h>
169 * #include <deal.II/lac/sparse_direct.h>
170 * #include <deal.II/particles/particle_handler.h>
171 * #include <deal.II/particles/data_out.h>
173 * #include <deal.II/numerics/data_out.h>
175 * #include <iostream>
177 * #include <deal.II/base/logstream.h>
178 * #include <deal.II/grid/grid_refinement.h>
185 * The following is the main
class. It resembles a variation of the @ref step_6
"step-6"
186 * principal
class, with the addition of information-specific stuff. It also
187 * has to deal with solving a vector-valued problem
for (c,lambda,f) as
188 * primal variable, dual variable, and right hand side, as explained
193 *
class InformationDensityMeshRefinement
196 * InformationDensityMeshRefinement ();
200 *
void compute_synthetic_measurements();
201 *
void bounce_measurement_points_to_cell_centers ();
202 *
void setup_system();
203 *
void assemble_system ();
205 *
void compute_information_content ();
206 *
void output_results (
const unsigned int cycle)
const;
207 *
void refine_grid ();
210 *
const double source_radius;
212 * std::vector<Point<dim>> detector_locations;
214 *
const double regularization_parameter;
231 * std::vector<Point<dim>> detector_locations_on_mesh;
232 * std::vector<double> measurement_values;
233 * std::vector<double> noise_level;
240 * InformationDensityMeshRefinement<dim>::InformationDensityMeshRefinement ()
243 * source_radius (0.2),
244 * regularization_parameter (10000),
254 * We have 50 detector points on an outer ring...
257 *
for (
unsigned int i=0; i<50; ++i)
261 * detector_locations.push_back (p);
266 * ...and another 50 detector points on an innner ring:
269 *
for (
unsigned int i=0; i<50; ++i)
273 * detector_locations.push_back (p);
278 * Generate the grid we will work on:
286 * The detector locations are
static, so we can already here
287 * generate a file that contains their locations. We use the
288 * particle framework to
do this,
using detector locations as
289 * particle locations.
295 *
for (
const auto &loc : detector_locations)
298 * new_particle.set_location(loc);
301 * Insert the particle. It is a lie that the particle is in
302 * the
first cell, but
nothing we
do actually cares about the
303 * cell a particle is in.
306 * particle_handler.insert_particle(new_particle,
312 * std::ofstream output(
"detector_locations.vtu");
313 * particle_out.write_vtu(output);
318 * While we
're generating output, also output the source location. Do this
319 * by outputting many (1000) points that indicate the perimeter of the source
323 * Particles::ParticleHandler<dim> particle_handler(triangulation,
324 * StaticMappingQ1<dim>::mapping);
326 * const unsigned int n_points = 1000;
327 * for (unsigned int i=0; i<n_points; ++i)
329 * Point<dim> loc = source_location;
330 * loc[0] += source_radius * std::cos(2*numbers::PI*i/n_points);
331 * loc[1] += source_radius * std::sin(2*numbers::PI*i/n_points);
333 * Particles::Particle<dim> new_particle;
334 * new_particle.set_location(loc);
335 * particle_handler.insert_particle(new_particle,
336 * triangulation.begin_active());
339 * Particles::DataOut<dim> particle_out;
340 * particle_out.build_patches(particle_handler);
341 * std::ofstream output("source_locations.vtu");
342 * particle_out.write_vtu(output);
350 * The following function solves a forward problem on a twice
351 * refined mesh to compute "synthetic data". Refining the mesh
352 * beyond the mesh used for the inverse problem avoids an
357 * void InformationDensityMeshRefinement<dim>::compute_synthetic_measurements ()
359 * std::cout << "Computing synthetic data by solving the forward problem..."
364 * Create a triangulation and DoFHandler that corresponds to a
365 * twice-refined mesh so that we obtain the synthetic data with
366 * higher accuracy than we do on the regular mesh used for all other
370 * Triangulation<dim> forward_triangulation;
371 * forward_triangulation.copy_triangulation (triangulation);
372 * forward_triangulation.refine_global (2);
374 * const FE_Q<dim> forward_fe (fe.base_element(0).degree);
375 * DoFHandler<dim> forward_dof_handler (forward_triangulation);
376 * forward_dof_handler.distribute_dofs (forward_fe);
378 * AffineConstraints<double> constraints;
379 * DoFTools::make_hanging_node_constraints(forward_dof_handler, constraints);
380 * constraints.close();
382 * SparsityPattern sparsity (forward_dof_handler.n_dofs(),
383 * forward_dof_handler.max_couplings_between_dofs());
384 * DoFTools::make_sparsity_pattern (forward_dof_handler, sparsity);
385 * constraints.condense (sparsity);
386 * sparsity.compress ();
388 * SparseMatrix<double> system_matrix (sparsity);
389 * Vector<double> system_rhs (forward_dof_handler.n_dofs());
391 * QGauss<dim> quadrature_formula(3);
392 * FEValues<dim> fe_values (forward_fe, quadrature_formula,
393 * update_values | update_gradients |
394 * update_quadrature_points | update_JxW_values);
396 * const unsigned int dofs_per_cell = forward_fe.dofs_per_cell;
397 * const unsigned int n_q_points = quadrature_formula.size();
401 * First assemble the system matrix and right hand side for the forward
406 * FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
407 * Vector<double> cell_rhs (dofs_per_cell);
408 * std::vector<unsigned int> local_dof_indices (dofs_per_cell);
410 * for (const auto &cell : forward_dof_handler.active_cell_iterators())
412 * fe_values.reinit (cell);
416 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
417 * for (unsigned int i=0; i<dofs_per_cell; ++i)
418 * for (unsigned int j=0; j<dofs_per_cell; ++j)
419 * cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
420 * fe_values.shape_grad(j,q_point)
422 * fe_values.shape_value(i,q_point) *
423 * (velocity * fe_values.shape_grad(j,q_point))
426 * fe_values.JxW(q_point);
427 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
428 * if (fe_values.quadrature_point(q_point).distance (source_location)
430 * for (unsigned int i=0; i<dofs_per_cell; ++i)
433 * fe_values.shape_value (i, q_point) *
434 * fe_values.JxW(q_point);
436 * cell->get_dof_indices (local_dof_indices);
437 * constraints.distribute_local_to_global (cell_matrix,
444 * std::map<unsigned int, double> boundary_values;
445 * VectorTools::interpolate_boundary_values (forward_dof_handler,
447 * Functions::ZeroFunction<dim>(),
449 * Vector<double> tmp (forward_dof_handler.n_dofs());
450 * MatrixTools::apply_boundary_values (boundary_values,
458 * Solve the forward problem and output it into its own VTU file :
461 * SparseDirectUMFPACK A_inverse;
462 * Vector<double> forward_solution (forward_dof_handler.n_dofs());
463 * forward_solution = system_rhs;
464 * A_inverse.solve(system_matrix, forward_solution);
466 * const double max_forward_solution = forward_solution.linfty_norm();
469 * DataOut<dim> data_out;
470 * data_out.attach_dof_handler (forward_dof_handler);
471 * data_out.add_data_vector (forward_solution, "c");
472 * data_out.build_patches (4);
474 * std::ofstream out ("forward-solution.vtu");
475 * data_out.write_vtu (out);
480 * Now evaluate the forward solution at the measurement points:
483 * for (const auto &p : detector_locations)
487 * same 10% noise level for all points
490 * noise_level.push_back (0.1 * max_forward_solution);
492 * const double z_n = VectorTools::point_value(forward_dof_handler, forward_solution, p);
493 * const double eps_n = Utilities::generate_normal_random_number(0, noise_level.back());
495 * measurement_values.push_back (z_n + eps_n);
498 * std::cout << std::endl;
504 * It will make our lives easier if we can always assume that detector
505 * locations are at cell centers, because then we can evaluate the
506 * solution there using a quadrature formula whose sole quadrature
507 * point lies at the center of a cell. That's of course not where the
508 *
"real" detector locations are, but it does not introduce a large
513 *
void InformationDensityMeshRefinement<dim>::bounce_measurement_points_to_cell_centers ()
515 * detector_locations_on_mesh = detector_locations;
516 *
for (
auto &p : detector_locations_on_mesh)
518 * for (const auto &cell :
triangulation.active_cell_iterators())
519 * if (cell->point_inside (p))
530 * The following
functions are all quite standard by what we have
531 * shown in @ref step_4
"step-4", @ref step_6
"step-6", and @ref step_22
"step-22" (to name just a few of the
532 * more typical programs):
536 * void InformationDensityMeshRefinement<dim>::setup_system ()
538 *
std::cout <<
"Setting up the linear system for the inverse problem..."
541 * dof_handler.distribute_dofs (fe);
544 * hanging_node_constraints.clear ();
546 * hanging_node_constraints);
547 * hanging_node_constraints.close();
549 * std::cout <<
" Number of active cells: "
552 * std::cout <<
" Number of degrees of freedom: "
553 * << dof_handler.n_dofs()
556 *
const std::vector<types::global_dof_index> dofs_per_component =
560 * hanging_node_constraints.condense(c_sparsity);
561 * sparsity_pattern.copy_from(c_sparsity);
563 * system_matrix.reinit (sparsity_pattern);
565 * solution.reinit (dofs_per_component);
566 * system_rhs.reinit (dofs_per_component);
572 *
void InformationDensityMeshRefinement<dim>::assemble_system ()
574 * std::cout <<
"Assembling the linear system for the inverse problem..."
583 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
584 *
const unsigned int n_q_points = quadrature_formula.size();
589 * std::vector<unsigned int> local_dof_indices (dofs_per_cell);
593 *
for (
const auto &cell : dof_handler.active_cell_iterators())
595 * fe_values.
reinit (cell);
599 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
600 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
602 *
const Tensor<1,dim> grad_phi_i = fe_values[c].gradient (i,q_point);
605 *
const double phi_i = fe_values[c].value (i,q_point);
606 *
const double psi_i = fe_values[
lambda].value (i,q_point);
607 *
const double chi_i = fe_values[f].value (i,q_point);
609 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
611 *
const Tensor<1,dim> grad_phi_j = fe_values[c].gradient (j,q_point);
614 *
const double phi_j = fe_values[c].value (j,q_point);
615 *
const double psi_j= fe_values[
lambda].value (j,q_point);
616 *
const double chi_j = fe_values[f].value (j,q_point);
619 * ((grad_phi_i * grad_phi_j
621 * phi_i * (velocity * grad_phi_j)
625 * grad_psi_i * grad_psi_j
627 * psi_i * (velocity * grad_psi_j)
631 * regularization_parameter * chi_i * chi_j
633 * fe_values.JxW (q_point));
635 *
for (
unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
636 *
if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
638 *
cell_matrix(i,j) += psi_i * phi_j / noise_level[n] / noise_level[n];
642 *
for (
unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
643 *
if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
644 * cell_rhs(i) += psi_i * measurement_values[n] / noise_level[n] / noise_level[n];
647 * cell->get_dof_indices (local_dof_indices);
648 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
650 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
651 * system_matrix.add (local_dof_indices[i],
652 * local_dof_indices[j],
655 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
659 * hanging_node_constraints.condense (system_matrix);
660 * hanging_node_constraints.condense (system_rhs);
662 * std::map<unsigned int,double> boundary_values;
663 * std::vector<bool> component_mask (3);
664 * component_mask[0] = component_mask[1] =
true;
665 * component_mask[2] =
false;
676 * std::cout << std::endl;
682 *
void InformationDensityMeshRefinement<dim>::solve ()
684 * std::cout <<
"Solving the linear system for the inverse problem..."
688 * solution = system_rhs;
689 * A_direct.
solve(system_matrix, solution);
691 * hanging_node_constraints.distribute (solution);
693 * std::cout << std::endl;
700 * This is really the only interesting function of
this program. It
701 * computes the
functions @f$h_K = A^{-1} s_K@f$
for each source function
702 * (corresponding to each cell of the mesh). To
do so, it
first
704 *
class to build an LU decomposition
for this matrix. Then it loops
705 * over all cells @f$K@f$ and computes the corresponding @f$h_K@f$ by applying
706 * the LU decomposition to a right hand side vector
for each @f$s_K@f$.
710 * The actual information content is then computed by evaluating these
711 * functions @f$h_K@f$ at measurement locations.
715 *
void InformationDensityMeshRefinement<dim>::compute_information_content ()
717 * std::cout <<
"Computing the information content..."
720 * information_content.reinit (
triangulation.n_active_cells());
722 *
const FE_Q<dim> information_fe (fe.base_element(0).degree);
724 * information_dof_handler.distribute_dofs (information_fe);
728 * constraints.close();
731 * information_dof_handler.max_couplings_between_dofs());
733 * constraints.condense (sparsity);
734 * sparsity.compress ();
740 *
const unsigned int dofs_per_cell = information_fe.dofs_per_cell;
741 *
const unsigned int n_q_points = quadrature_formula.size();
745 * First build the forward
operator
749 *
FEValues<dim> fe_values (information_fe, quadrature_formula,
754 * std::vector<unsigned int> local_dof_indices (dofs_per_cell);
756 *
for (
const auto &cell : information_dof_handler.active_cell_iterators())
758 * fe_values.
reinit (cell);
761 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
762 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
763 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
764 *
cell_matrix(i,j) += (fe_values.shape_grad (i,q_point) *
765 * fe_values.shape_grad(j,q_point)
767 * fe_values.shape_value(i,q_point) *
768 * (velocity * fe_values.shape_grad(j,q_point))) *
769 * fe_values.JxW(q_point);
771 * cell->distribute_local_to_global (cell_matrix,
775 * constraints.condense (system_matrix);
777 * std::map<unsigned int, double> boundary_values;
799 * Now compute the solutions corresponding to the possible
800 * sources. Each source is active on exactly one cell.
804 * As mentioned in the paper,
this is a trivially
parallel job, so
805 * we send the computations
for each of these cells onto a separate
806 * task and let the OS schedule them onto individual processor
817 * std::vector<unsigned int> local_dof_indices (dofs_per_cell);
820 * cell = information_dof_handler.begin_active();
821 * std::advance (cell, K);
823 *
FEValues<dim> fe_values (information_fe, quadrature_formula,
827 * fe_values.reinit (cell);
830 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
831 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
832 * cell_rhs(i) += fe_values.shape_value (i,q_point) *
833 * fe_values.JxW(q_point);
835 * cell->distribute_local_to_global (cell_rhs,
838 * constraints.condense (rhs);
840 * A_inverse.solve(rhs);
842 * constraints.distribute (rhs);
846 * Having computed the forward solutions
847 * corresponding to
this source term, evaluate its
848 * contribution to the information content on all
849 * cells of the mesh by taking into account the
850 * detector locations. We add these into global
851 * objects, so we have to guard access to the
855 *
static std::mutex m;
856 * std::lock_guard<std::mutex> g(m);
859 * information_content(K) = regularization_parameter * cell->measure() * cell->measure();
860 * std::vector<double> local_h_K_values (n_q_points);
861 *
for (
const auto &cell : information_dof_handler.active_cell_iterators())
863 * fe_values.
reinit (cell);
864 * fe_values.get_function_values (rhs, local_h_K_values);
866 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
867 *
for (
unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
868 *
if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
869 * information_content(K) += local_h_K_values[q_point]
870 * * local_h_K_values[q_point]
871 * / noise_level[n] / noise_level[n];
883 * std::cout << std::endl;
890 * Create graphical output
for all of the principal variables of the
891 * problem (c,lambda,f) as well as
for the information content and
892 * density. Then also output the various blocks of the
matrix so we
897 *
void InformationDensityMeshRefinement<dim>::output_results (
const unsigned int cycle)
const
899 * std::cout <<
"Outputting solutions..." << std::flush;
903 * std::vector<std::string> names;
904 * names.push_back (
"forward_solution");
905 * names.push_back (
"adjoint_solution");
906 * names.push_back (
"recovered_parameter");
909 * data_out.add_data_vector (solution, names);
910 * data_out.add_data_vector (information_content,
"information_content");
913 *
for (
const auto &cell :
triangulation.active_cell_iterators())
914 * information_density(cell->active_cell_index())
915 * =
std::
sqrt(information_content(cell->active_cell_index())) / cell->measure();
916 * data_out.add_data_vector (information_density,
"information_density");
918 * data_out.build_patches ();
920 * std::string filename =
"solution-";
921 * filename += (
'0'+cycle);
922 * filename +=
".vtu";
924 * std::ofstream output (filename.c_str());
925 * data_out.write_vtu (output);
930 * Now output the individual blocks of the
matrix into files.
933 *
auto write_block = [&](
const unsigned int block_i,
934 *
const unsigned int block_j,
935 *
const std::string &filename)
937 * std::ofstream o(filename);
938 * system_matrix.block(block_i,block_j).print (o);
940 * write_block(0,0,
"matrix-" + std::to_string(cycle) +
"-A.txt");
941 * write_block(0,2,
"matrix-" + std::to_string(cycle) +
"-B.txt");
942 * write_block(1,0,
"matrix-" + std::to_string(cycle) +
"-C.txt");
943 * write_block(2,2,
"matrix-" + std::to_string(cycle) +
"-M.txt");
945 * std::cout << std::endl;
952 * The following is then a function that refines the mesh based on the
953 * refinement criteria described in the paper. Which criterion to use
954 * is determined by which
value the `refinement_criterion` variable
959 *
void InformationDensityMeshRefinement<dim>::refine_grid ()
961 * std::cout <<
"Refining the mesh..." << std::endl;
963 *
enum RefinementCriterion
966 * information_content,
970 *
const RefinementCriterion refinement_criterion = information_content;
972 *
switch (refinement_criterion)
980 *
case information_content:
983 * this->information_content,
998 * std::vector<double> lambda_values (quadrature.size());
999 * std::vector<double> f_values (quadrature.size());
1001 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1003 * fe_values.
reinit (cell);
1004 * fe_values[
lambda].get_function_values (solution, lambda_values);
1005 * fe_values[f].get_function_values (solution, f_values);
1007 *
for (
unsigned int q=0; q<quadrature.size(); ++q)
1008 * refinement_indicators(cell->active_cell_index())
1009 * += (std::fabs (regularization_parameter * f_values[q]
1012 * * fe_values.JxW(q));
1016 * refinement_indicators,
1029 * refinement_indicators,
1033 * and
scale it to obtain an error indicator.
1036 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1037 * refinement_indicators[cell->active_cell_index()] *=
1042 * refinement_indicators,
1049 *
Assert (
false, ExcInternalError());
1052 * bounce_measurement_points_to_cell_centers ();
1055 * std::cout << std::endl;
1061 *
template <
int dim>
1062 *
void InformationDensityMeshRefinement<dim>::run ()
1064 * std::cout <<
"Solving problem in " << dim <<
" space dimensions." << std::endl;
1066 * compute_synthetic_measurements ();
1067 * bounce_measurement_points_to_cell_centers ();
1069 *
for (
unsigned int cycle=0; cycle<7; ++cycle)
1071 * std::cout <<
"---------- Cycle " << cycle <<
" ------------" << std::endl;
1074 * assemble_system ();
1076 * compute_information_content ();
1077 * output_results (cycle);
1090 * InformationDensityMeshRefinement<2> information_density_mesh_refinement;
1091 * information_density_mesh_refinement.run ();
1093 *
catch (std::exception &exc)
1095 * std::cerr << std::endl << std::endl
1096 * <<
"----------------------------------------------------"
1098 * std::cerr <<
"Exception on processing: " << std::endl
1099 * << exc.what() << std::endl
1100 * <<
"Aborting!" << std::endl
1101 * <<
"----------------------------------------------------"
1108 * std::cerr << std::endl << std::endl
1109 * <<
"----------------------------------------------------"
1111 * std::cerr <<
"Unknown exception!" << std::endl
1112 * <<
"Aborting!" << std::endl
1113 * <<
"----------------------------------------------------"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
unsigned int depth_console(const unsigned int n)
void build_patches(const Particles::ParticleHandler< dim, spacedim > &particles, const std::vector< std::string > &data_component_names={}, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretations={})
void solve(Vector< double > &rhs_and_solution, const bool transpose=false) const
void factorize(const Matrix &matrix)
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Task< RT > new_task(const std::function< RT()> &function)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void approximate_gradient(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
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.
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.)
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 reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
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 > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
const ::Triangulation< dim, spacedim > & tria