Reference documentation for deal.II version 9.5.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
information_based_mesh_refinement.h
Go to the documentation of this file.
1
144 *  
145 *  
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>
172 *  
173 *   #include <deal.II/numerics/data_out.h>
174 *   #include <fstream>
175 *   #include <iostream>
176 *  
177 *   #include <deal.II/base/logstream.h>
178 *   #include <deal.II/grid/grid_refinement.h>
179 *  
180 *   using namespace dealii;
181 *  
182 *  
183 * @endcode
184 *
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
189 * in the paper.
190 *
191 * @code
192 *   template <int dim>
193 *   class InformationDensityMeshRefinement
194 *   {
195 *   public:
196 *   InformationDensityMeshRefinement ();
197 *   void run ();
198 *  
199 *   private:
200 *   void compute_synthetic_measurements();
201 *   void bounce_measurement_points_to_cell_centers ();
202 *   void setup_system();
203 *   void assemble_system ();
204 *   void solve ();
205 *   void compute_information_content ();
206 *   void output_results (const unsigned int cycle) const;
207 *   void refine_grid ();
208 *  
209 *   const Point<dim> source_location;
210 *   const double source_radius;
211 *  
212 *   std::vector<Point<dim>> detector_locations;
213 *  
214 *   const double regularization_parameter;
215 *   Tensor<1,dim> velocity;
216 *  
218 *   FESystem<dim> fe;
219 *   DoFHandler<dim> dof_handler;
220 *  
221 *   AffineConstraints<double> hanging_node_constraints;
222 *  
223 *   BlockSparsityPattern sparsity_pattern;
224 *   BlockSparseMatrix<double> system_matrix;
225 *  
226 *   BlockVector<double> solution;
227 *   BlockVector<double> system_rhs;
228 *  
229 *   Vector<double> information_content;
230 *  
231 *   std::vector<Point<dim>> detector_locations_on_mesh;
232 *   std::vector<double> measurement_values;
233 *   std::vector<double> noise_level;
234 *   };
235 *  
236 *  
237 *  
238 *  
239 *   template <int dim>
240 *   InformationDensityMeshRefinement<dim>::InformationDensityMeshRefinement ()
241 *   :
242 *   source_location (Point<dim>(-0.25,0)),
243 *   source_radius (0.2),
244 *   regularization_parameter (10000),
245 *   fe (FE_Q<dim>(3), 1, // c
246 *   FE_Q<dim>(3), 1, // lambda
247 *   FE_DGQ<dim>(0), 1), // f
248 *   dof_handler (triangulation)
249 *   {
250 *   velocity[0] = 100;
251 *  
252 * @endcode
253 *
254 * We have 50 detector points on an outer ring...
255 *
256 * @code
257 *   for (unsigned int i=0; i<50; ++i)
258 *   {
259 *   const Point<dim> p (0.6 * std::sin(2*numbers::PI * i/50),
260 *   0.6 * std::cos(2*numbers::PI * i/50));
261 *   detector_locations.push_back (p);
262 *   }
263 *  
264 * @endcode
265 *
266 * ...and another 50 detector points on an innner ring:
267 *
268 * @code
269 *   for (unsigned int i=0; i<50; ++i)
270 *   {
271 *   const Point<dim> p (0.2 * std::sin(2*numbers::PI * i/50),
272 *   0.2 * std::cos(2*numbers::PI * i/50));
273 *   detector_locations.push_back (p);
274 *   }
275 *  
276 * @endcode
277 *
278 * Generate the grid we will work on:
279 *
280 * @code
282 *   triangulation.refine_global (4);
283 *  
284 * @endcode
285 *
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.
290 *
291 * @code
292 *   {
295 *   for (const auto &loc : detector_locations)
296 *   {
297 *   Particles::Particle<dim> new_particle;
298 *   new_particle.set_location(loc);
299 * @endcode
300 *
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.
304 *
305 * @code
306 *   particle_handler.insert_particle(new_particle,
307 *   triangulation.begin_active());
308 *   }
309 *  
310 *   Particles::DataOut<dim> particle_out;
311 *   particle_out.build_patches(particle_handler);
312 *   std::ofstream output("detector_locations.vtu");
313 *   particle_out.write_vtu(output);
314 *   }
315 *  
316 * @endcode
317 *
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
320 *
321 * @code
322 *   {
323 *   Particles::ParticleHandler<dim> particle_handler(triangulation,
324 *   StaticMappingQ1<dim>::mapping);
325 *  
326 *   const unsigned int n_points = 1000;
327 *   for (unsigned int i=0; i<n_points; ++i)
328 *   {
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);
332 *  
333 *   Particles::Particle<dim> new_particle;
334 *   new_particle.set_location(loc);
335 *   particle_handler.insert_particle(new_particle,
336 *   triangulation.begin_active());
337 *   }
338 *  
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);
343 *   }
344 *   }
345 *  
346 *  
347 *  
348 * @endcode
349 *
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
353 * inverse crime.
354 *
355 * @code
356 *   template <int dim>
357 *   void InformationDensityMeshRefinement<dim>::compute_synthetic_measurements ()
358 *   {
359 *   std::cout << "Computing synthetic data by solving the forward problem..."
360 *   << std::flush;
361 *  
362 * @endcode
363 *
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
367 * computations.
368 *
369 * @code
370 *   Triangulation<dim> forward_triangulation;
371 *   forward_triangulation.copy_triangulation (triangulation);
372 *   forward_triangulation.refine_global (2);
373 *  
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);
377 *  
378 *   AffineConstraints<double> constraints;
379 *   DoFTools::make_hanging_node_constraints(forward_dof_handler, constraints);
380 *   constraints.close();
381 *  
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 ();
387 *  
388 *   SparseMatrix<double> system_matrix (sparsity);
389 *   Vector<double> system_rhs (forward_dof_handler.n_dofs());
390 *  
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);
395 *  
396 *   const unsigned int dofs_per_cell = forward_fe.dofs_per_cell;
397 *   const unsigned int n_q_points = quadrature_formula.size();
398 *  
399 * @endcode
400 *
401 * First assemble the system matrix and right hand side for the forward
402 * problem:
403 *
404 * @code
405 *   {
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);
409 *  
410 *   for (const auto &cell : forward_dof_handler.active_cell_iterators())
411 *   {
412 *   fe_values.reinit (cell);
413 *   cell_matrix = 0;
414 *   cell_rhs = 0;
415 *  
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)
421 *   +
422 *   fe_values.shape_value(i,q_point) *
423 *   (velocity * fe_values.shape_grad(j,q_point))
424 *   )
425 *   *
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)
429 *   < source_radius)
430 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
431 *   cell_rhs(i) +=
432 *   1.0 *
433 *   fe_values.shape_value (i, q_point) *
434 *   fe_values.JxW(q_point);
435 *  
436 *   cell->get_dof_indices (local_dof_indices);
437 *   constraints.distribute_local_to_global (cell_matrix,
438 *   cell_rhs,
439 *   local_dof_indices,
440 *   system_matrix,
441 *   system_rhs);
442 *   }
443 *  
444 *   std::map<unsigned int, double> boundary_values;
445 *   VectorTools::interpolate_boundary_values (forward_dof_handler,
446 *   0,
447 *   Functions::ZeroFunction<dim>(),
448 *   boundary_values);
449 *   Vector<double> tmp (forward_dof_handler.n_dofs());
450 *   MatrixTools::apply_boundary_values (boundary_values,
451 *   system_matrix,
452 *   tmp,
453 *   system_rhs);
454 *   }
455 *  
456 * @endcode
457 *
458 * Solve the forward problem and output it into its own VTU file :
459 *
460 * @code
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);
465 *  
466 *   const double max_forward_solution = forward_solution.linfty_norm();
467 *  
468 *   {
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);
473 *  
474 *   std::ofstream out ("forward-solution.vtu");
475 *   data_out.write_vtu (out);
476 *   }
477 *  
478 * @endcode
479 *
480 * Now evaluate the forward solution at the measurement points:
481 *
482 * @code
483 *   for (const auto &p : detector_locations)
484 *   {
485 * @endcode
486 *
487 * same 10% noise level for all points
488 *
489 * @code
490 *   noise_level.push_back (0.1 * max_forward_solution);
491 *  
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());
494 *  
495 *   measurement_values.push_back (z_n + eps_n);
496 *   }
497 *  
498 *   std::cout << std::endl;
499 *   }
500 *  
501 *  
502 * @endcode
503 *
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
509 * error to do this.
510 *
511 * @code
512 *   template <int dim>
513 *   void InformationDensityMeshRefinement<dim>::bounce_measurement_points_to_cell_centers ()
514 *   {
515 *   detector_locations_on_mesh = detector_locations;
516 *   for (auto &p : detector_locations_on_mesh)
517 *   {
518 *   for (const auto &cell : triangulation.active_cell_iterators())
519 *   if (cell->point_inside (p))
520 *   {
521 *   p = cell->center();
522 *   break;
523 *   }
524 *   }
525 *   }
526 *  
527 *  
528 * @endcode
529 *
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):
533 *
534 * @code
535 *   template <int dim>
536 *   void InformationDensityMeshRefinement<dim>::setup_system ()
537 *   {
538 *   std::cout << "Setting up the linear system for the inverse problem..."
539 *   << std::endl;
540 *  
541 *   dof_handler.distribute_dofs (fe);
542 *   DoFRenumbering::component_wise (dof_handler);
543 *  
544 *   hanging_node_constraints.clear ();
546 *   hanging_node_constraints);
547 *   hanging_node_constraints.close();
548 *  
549 *   std::cout << " Number of active cells: "
550 *   << triangulation.n_active_cells()
551 *   << std::endl;
552 *   std::cout << " Number of degrees of freedom: "
553 *   << dof_handler.n_dofs()
554 *   << std::endl;
555 *  
556 *   const std::vector<types::global_dof_index> dofs_per_component =
558 *   BlockDynamicSparsityPattern c_sparsity(dofs_per_component,dofs_per_component);
559 *   DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
560 *   hanging_node_constraints.condense(c_sparsity);
561 *   sparsity_pattern.copy_from(c_sparsity);
562 *  
563 *   system_matrix.reinit (sparsity_pattern);
564 *  
565 *   solution.reinit (dofs_per_component);
566 *   system_rhs.reinit (dofs_per_component);
567 *   }
568 *  
569 *  
570 *  
571 *   template <int dim>
572 *   void InformationDensityMeshRefinement<dim>::assemble_system ()
573 *   {
574 *   std::cout << "Assembling the linear system for the inverse problem..."
575 *   << std::flush;
576 *  
577 *   QGauss<dim> quadrature_formula(3);
578 *  
579 *   FEValues<dim> fe_values (fe, quadrature_formula,
582 *  
583 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
584 *   const unsigned int n_q_points = quadrature_formula.size();
585 *  
586 *   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
587 *   Vector<double> cell_rhs (dofs_per_cell);
588 *  
589 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
590 *  
591 *   FEValuesExtractors::Scalar c(0), lambda(1), f(2);
592 *  
593 *   for (const auto &cell : dof_handler.active_cell_iterators())
594 *   {
595 *   fe_values.reinit (cell);
596 *   cell_matrix = 0;
597 *   cell_rhs = 0;
598 *  
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)
601 *   {
602 *   const Tensor<1,dim> grad_phi_i = fe_values[c].gradient (i,q_point);
603 *   const Tensor<1,dim> grad_psi_i = fe_values[lambda].gradient (i,q_point);
604 *  
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);
608 *  
609 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
610 *   {
611 *   const Tensor<1,dim> grad_phi_j = fe_values[c].gradient (j,q_point);
612 *   const Tensor<1,dim> grad_psi_j = fe_values[lambda].gradient (j,q_point);
613 *  
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);
617 *  
618 *   cell_matrix(i,j) +=
619 *   ((grad_phi_i * grad_phi_j
620 *   +
621 *   phi_i * (velocity * grad_phi_j)
622 *   -
623 *   phi_i * chi_j
624 *   +
625 *   grad_psi_i * grad_psi_j
626 *   -
627 *   psi_i * (velocity * grad_psi_j)
628 *   -
629 *   chi_i * psi_j
630 *   +
631 *   regularization_parameter * chi_i * chi_j
632 *   ) *
633 *   fe_values.JxW (q_point));
634 *  
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)
637 *   {
638 *   cell_matrix(i,j) += psi_i * phi_j / noise_level[n] / noise_level[n];
639 *   }
640 *   }
641 *  
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];
645 *   }
646 *  
647 *   cell->get_dof_indices (local_dof_indices);
648 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
649 *   {
650 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
651 *   system_matrix.add (local_dof_indices[i],
652 *   local_dof_indices[j],
653 *   cell_matrix(i,j));
654 *  
655 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
656 *   }
657 *   }
658 *  
659 *   hanging_node_constraints.condense (system_matrix);
660 *   hanging_node_constraints.condense (system_rhs);
661 *  
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;
667 *   0,
669 *   boundary_values,
670 *   component_mask);
671 *   MatrixTools::apply_boundary_values (boundary_values,
672 *   system_matrix,
673 *   solution,
674 *   system_rhs);
675 *  
676 *   std::cout << std::endl;
677 *   }
678 *  
679 *  
680 *  
681 *   template <int dim>
682 *   void InformationDensityMeshRefinement<dim>::solve ()
683 *   {
684 *   std::cout << "Solving the linear system for the inverse problem..."
685 *   << std::flush;
686 *  
687 *   SparseDirectUMFPACK A_direct;
688 *   solution = system_rhs;
689 *   A_direct.solve(system_matrix, solution);
690 *  
691 *   hanging_node_constraints.distribute (solution);
692 *  
693 *   std::cout << std::endl;
694 *   }
695 *  
696 *  
697 *  
698 * @endcode
699 *
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
703 * computes the forward matrix @f$A@f$ and uses the SparseDirectUMFPACK
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$.
707 *
708
709 *
710 * The actual information content is then computed by evaluating these
711 * functions @f$h_K@f$ at measurement locations.
712 *
713 * @code
714 *   template <int dim>
715 *   void InformationDensityMeshRefinement<dim>::compute_information_content ()
716 *   {
717 *   std::cout << "Computing the information content..."
718 *   << std::flush;
719 *  
720 *   information_content.reinit (triangulation.n_active_cells());
721 *  
722 *   const FE_Q<dim> information_fe (fe.base_element(0).degree);
723 *   DoFHandler<dim> information_dof_handler (triangulation);
724 *   information_dof_handler.distribute_dofs (information_fe);
725 *  
726 *   AffineConstraints<double> constraints;
727 *   DoFTools::make_hanging_node_constraints(information_dof_handler, constraints);
728 *   constraints.close();
729 *  
730 *   SparsityPattern sparsity (information_dof_handler.n_dofs(),
731 *   information_dof_handler.max_couplings_between_dofs());
732 *   DoFTools::make_sparsity_pattern (information_dof_handler, sparsity);
733 *   constraints.condense (sparsity);
734 *   sparsity.compress ();
735 *  
736 *   SparseMatrix<double> system_matrix (sparsity);
737 *  
738 *   QGauss<dim> quadrature_formula(3);
739 *  
740 *   const unsigned int dofs_per_cell = information_fe.dofs_per_cell;
741 *   const unsigned int n_q_points = quadrature_formula.size();
742 *  
743 * @endcode
744 *
745 * First build the forward operator
746 *
747 * @code
748 *   {
749 *   FEValues<dim> fe_values (information_fe, quadrature_formula,
752 *  
753 *   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
754 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
755 *  
756 *   for (const auto &cell : information_dof_handler.active_cell_iterators())
757 *   {
758 *   fe_values.reinit (cell);
759 *   cell_matrix = 0;
760 *  
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)
766 *   +
767 *   fe_values.shape_value(i,q_point) *
768 *   (velocity * fe_values.shape_grad(j,q_point))) *
769 *   fe_values.JxW(q_point);
770 *  
771 *   cell->distribute_local_to_global (cell_matrix,
772 *   system_matrix);
773 *   }
774 *  
775 *   constraints.condense (system_matrix);
776 *  
777 *   std::map<unsigned int, double> boundary_values;
778 *   VectorTools::interpolate_boundary_values (information_dof_handler,
779 *   0,
781 *   boundary_values);
782 *   Vector<double> tmp (information_dof_handler.n_dofs());
783 *   MatrixTools::apply_boundary_values (boundary_values,
784 *   system_matrix,
785 *   tmp,
786 *   tmp);
787 *   }
788 *  
789 * @endcode
790 *
791 * Then factorize
792 *
793 * @code
794 *   SparseDirectUMFPACK A_inverse;
795 *   A_inverse.factorize(system_matrix);
796 *  
797 * @endcode
798 *
799 * Now compute the solutions corresponding to the possible
800 * sources. Each source is active on exactly one cell.
801 *
802
803 *
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
807 * cores.
808 *
809 * @code
810 *   Threads::TaskGroup<void> tasks;
811 *   for (unsigned int K=0; K<triangulation.n_active_cells(); ++K)
812 *   tasks +=
813 *   Threads::new_task([&,K]()
814 *   {
815 *   Vector<double> rhs (information_dof_handler.n_dofs());
816 *   Vector<double> cell_rhs (dofs_per_cell);
817 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
818 *  
820 *   cell = information_dof_handler.begin_active();
821 *   std::advance (cell, K);
822 *  
823 *   FEValues<dim> fe_values (information_fe, quadrature_formula,
824 *   update_values |
826 *  
827 *   fe_values.reinit (cell);
828 *   cell_rhs = 0;
829 *  
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);
834 *  
835 *   cell->distribute_local_to_global (cell_rhs,
836 *   rhs);
837 *  
838 *   constraints.condense (rhs);
839 *  
840 *   A_inverse.solve(rhs);
841 *  
842 *   constraints.distribute (rhs);
843 *  
844 * @endcode
845 *
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
852 * global object:
853 *
854 * @code
855 *   static std::mutex m;
856 *   std::lock_guard<std::mutex> g(m);
857 *  
858 *  
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())
862 *   {
863 *   fe_values.reinit (cell);
864 *   fe_values.get_function_values (rhs, local_h_K_values);
865 *  
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];
872 *   }
873 *   }
874 *   );
875 *  
876 * @endcode
877 *
878 * And wait:
879 *
880 * @code
881 *   tasks.join_all();
882 *  
883 *   std::cout << std::endl;
884 *   }
885 *  
886 *  
887 *  
888 * @endcode
889 *
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
893 * can compute the eigenvalues of the H matrix mentioned in the paper.
894 *
895 * @code
896 *   template <int dim>
897 *   void InformationDensityMeshRefinement<dim>::output_results (const unsigned int cycle) const
898 *   {
899 *   std::cout << "Outputting solutions..." << std::flush;
900 *  
901 *   DataOut<dim> data_out;
902 *  
903 *   std::vector<std::string> names;
904 *   names.push_back ("forward_solution");
905 *   names.push_back ("adjoint_solution");
906 *   names.push_back ("recovered_parameter");
907 *  
908 *   data_out.attach_dof_handler (dof_handler);
909 *   data_out.add_data_vector (solution, names);
910 *   data_out.add_data_vector (information_content, "information_content");
911 *  
912 *   Vector<double> information_density (triangulation.n_active_cells());
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");
917 *  
918 *   data_out.build_patches ();
919 *  
920 *   std::string filename = "solution-";
921 *   filename += ('0'+cycle);
922 *   filename += ".vtu";
923 *  
924 *   std::ofstream output (filename.c_str());
925 *   data_out.write_vtu (output);
926 *  
927 *  
928 * @endcode
929 *
930 * Now output the individual blocks of the matrix into files.
931 *
932 * @code
933 *   auto write_block = [&](const unsigned int block_i,
934 *   const unsigned int block_j,
935 *   const std::string &filename)
936 *   {
937 *   std::ofstream o(filename);
938 *   system_matrix.block(block_i,block_j).print (o);
939 *   };
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");
944 *  
945 *   std::cout << std::endl;
946 *   }
947 *  
948 *  
949 *  
950 * @endcode
951 *
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
955 * is set to.
956 *
957 * @code
958 *   template <int dim>
959 *   void InformationDensityMeshRefinement<dim>::refine_grid ()
960 *   {
961 *   std::cout << "Refining the mesh..." << std::endl;
962 *  
963 *   enum RefinementCriterion
964 *   {
965 *   global,
966 *   information_content,
967 *   indicator,
968 *   smoothness
969 *   };
970 *   const RefinementCriterion refinement_criterion = information_content;
971 *  
972 *   switch (refinement_criterion)
973 *   {
974 *   case global:
975 *   {
976 *   triangulation.refine_global();
977 *   break;
978 *   }
979 *  
980 *   case information_content:
981 *   {
983 *   this->information_content,
984 *   0.2, 0.05);
985 *   triangulation.execute_coarsening_and_refinement ();
986 *   break;
987 *   }
988 *  
989 *   case indicator:
990 *   {
991 *   Vector<double> refinement_indicators (triangulation.n_active_cells());
992 *  
993 *   QGauss<dim> quadrature(3);
994 *   FEValues<dim> fe_values (fe, quadrature, update_values | update_JxW_values);
995 *  
997 *  
998 *   std::vector<double> lambda_values (quadrature.size());
999 *   std::vector<double> f_values (quadrature.size());
1000 *  
1001 *   for (const auto &cell : dof_handler.active_cell_iterators())
1002 *   {
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);
1006 *  
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]
1010 *   -
1011 *   lambda_values[q])
1012 *   * fe_values.JxW(q));
1013 *   }
1014 *  
1016 *   refinement_indicators,
1017 *   0.2, 0.05);
1018 *   triangulation.execute_coarsening_and_refinement ();
1019 *   break;
1020 *   }
1021 *  
1022 *  
1023 *   case smoothness:
1024 *   {
1025 *   Vector<float> refinement_indicators (triangulation.n_active_cells());
1026 *  
1028 *   solution,
1029 *   refinement_indicators,
1030 *   /*component=*/2);
1031 * @endcode
1032 *
1033 * and scale it to obtain an error indicator.
1034 *
1035 * @code
1036 *   for (const auto &cell : triangulation.active_cell_iterators())
1037 *   refinement_indicators[cell->active_cell_index()] *=
1038 *   std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
1039 *  
1040 *  
1042 *   refinement_indicators,
1043 *   0.2, 0.05);
1044 *   triangulation.execute_coarsening_and_refinement ();
1045 *   break;
1046 *   }
1047 *  
1048 *   default:
1049 *   Assert (false, ExcInternalError());
1050 *   }
1051 *  
1052 *   bounce_measurement_points_to_cell_centers ();
1053 *  
1054 *  
1055 *   std::cout << std::endl;
1056 *   }
1057 *  
1058 *  
1059 *  
1060 *  
1061 *   template <int dim>
1062 *   void InformationDensityMeshRefinement<dim>::run ()
1063 *   {
1064 *   std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
1065 *  
1066 *   compute_synthetic_measurements ();
1067 *   bounce_measurement_points_to_cell_centers ();
1068 *  
1069 *   for (unsigned int cycle=0; cycle<7; ++cycle)
1070 *   {
1071 *   std::cout << "---------- Cycle " << cycle << " ------------" << std::endl;
1072 *  
1073 *   setup_system ();
1074 *   assemble_system ();
1075 *   solve ();
1076 *   compute_information_content ();
1077 *   output_results (cycle);
1078 *   refine_grid ();
1079 *   }
1080 *   }
1081 *  
1082 *  
1083 *  
1084 *   int main ()
1085 *   {
1086 *   try
1087 *   {
1088 *   deallog.depth_console (0);
1089 *  
1090 *   InformationDensityMeshRefinement<2> information_density_mesh_refinement;
1091 *   information_density_mesh_refinement.run ();
1092 *   }
1093 *   catch (std::exception &exc)
1094 *   {
1095 *   std::cerr << std::endl << std::endl
1096 *   << "----------------------------------------------------"
1097 *   << std::endl;
1098 *   std::cerr << "Exception on processing: " << std::endl
1099 *   << exc.what() << std::endl
1100 *   << "Aborting!" << std::endl
1101 *   << "----------------------------------------------------"
1102 *   << std::endl;
1103 *  
1104 *   return 1;
1105 *   }
1106 *   catch (...)
1107 *   {
1108 *   std::cerr << std::endl << std::endl
1109 *   << "----------------------------------------------------"
1110 *   << std::endl;
1111 *   std::cerr << "Unknown exception!" << std::endl
1112 *   << "Aborting!" << std::endl
1113 *   << "----------------------------------------------------"
1114 *   << std::endl;
1115 *   return 1;
1116 *   }
1117 *  
1118 *   return 0;
1119 *   }
1120 * @endcode
1121
1122
1123*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Definition fe_q.h:551
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:350
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={})
Definition data_out.cc:28
Definition point.h:112
void solve(Vector< double > &rhs_and_solution, const bool transpose=false) const
void factorize(const Matrix &matrix)
Point< 3 > center
Point< 2 > first
Definition grid_out.cc:4615
__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)
LogStream deallog
Definition logstream.cc:37
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 >())
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
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())
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
double diameter(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:88
@ 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.)
Definition advection.h:75
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
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
Definition numbers.h:259
STL namespace.
::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