Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
information_based_mesh_refinement.h
Go to the documentation of this file.
1
137 *  
138 *   #include <deal.II/base/numbers.h>
139 *   #include <deal.II/base/thread_management.h>
140 *   #include <deal.II/grid/tria.h>
141 *   #include <deal.II/dofs/dof_handler.h>
142 *   #include <deal.II/grid/grid_generator.h>
143 *   #include <deal.II/grid/tria_accessor.h>
144 *   #include <deal.II/grid/tria_iterator.h>
145 *   #include <deal.II/dofs/dof_accessor.h>
146 *   #include <deal.II/dofs/dof_renumbering.h>
147 *   #include <deal.II/fe/fe_q.h>
148 *   #include <deal.II/fe/fe_dgq.h>
149 *   #include <deal.II/fe/fe_system.h>
150 *   #include <deal.II/dofs/dof_tools.h>
151 *   #include <deal.II/fe/fe_values.h>
152 *   #include <deal.II/base/quadrature_lib.h>
153 *   #include <deal.II/base/function.h>
154 *   #include <deal.II/numerics/vector_tools.h>
155 *   #include <deal.II/numerics/matrix_tools.h>
156 *   #include <deal.II/numerics/derivative_approximation.h>
157 *   #include <deal.II/lac/block_vector.h>
158 *   #include <deal.II/lac/full_matrix.h>
159 *   #include <deal.II/lac/block_sparse_matrix.h>
160 *   #include <deal.II/lac/block_sparsity_pattern.h>
161 *   #include <deal.II/lac/sparse_direct.h>
162 *   #include <deal.II/particles/particle_handler.h>
163 *   #include <deal.II/particles/data_out.h>
164 *  
165 *   #include <deal.II/numerics/data_out.h>
166 *   #include <fstream>
167 *   #include <iostream>
168 *  
169 *   #include <deal.II/base/logstream.h>
170 *   #include <deal.II/grid/grid_refinement.h>
171 *  
172 *   using namespace dealii;
173 *  
174 *  
175 * @endcode
176 *
177 * The following is the main class. It resembles a variation of the @ref step_6 "step-6"
178 * principal class, with the addition of information-specific stuff. It also
179 * has to deal with solving a vector-valued problem for (c,lambda,f) as
180 * primal variable, dual variable, and right hand side, as explained
181 * in the paper.
182 *
183 * @code
184 *   template <int dim>
185 *   class InformationDensityMeshRefinement
186 *   {
187 *   public:
188 *   InformationDensityMeshRefinement ();
189 *   void run ();
190 *  
191 *   private:
192 *   void compute_synthetic_measurements();
193 *   void bounce_measurement_points_to_cell_centers ();
194 *   void setup_system();
195 *   void assemble_system ();
196 *   void solve ();
197 *   void compute_information_content ();
198 *   void output_results (const unsigned int cycle) const;
199 *   void refine_grid ();
200 *  
201 *   const Point<dim> source_location;
202 *   const double source_radius;
203 *  
204 *   std::vector<Point<dim>> detector_locations;
205 *  
206 *   const double regularization_parameter;
207 *   Tensor<1,dim> velocity;
208 *  
210 *   FESystem<dim> fe;
211 *   DoFHandler<dim> dof_handler;
212 *  
213 *   AffineConstraints<double> hanging_node_constraints;
214 *  
215 *   BlockSparsityPattern sparsity_pattern;
216 *   BlockSparseMatrix<double> system_matrix;
217 *  
218 *   BlockVector<double> solution;
219 *   BlockVector<double> system_rhs;
220 *  
221 *   Vector<double> information_content;
222 *  
223 *   std::vector<Point<dim>> detector_locations_on_mesh;
224 *   std::vector<double> measurement_values;
225 *   std::vector<double> noise_level;
226 *   };
227 *  
228 *  
229 *  
230 *  
231 *   template <int dim>
232 *   InformationDensityMeshRefinement<dim>::InformationDensityMeshRefinement ()
233 *   :
234 *   source_location (Point<dim>(-0.25,0)),
235 *   source_radius (0.2),
236 *   regularization_parameter (10000),
237 *   fe (FE_Q<dim>(3), 1, // c
238 *   FE_Q<dim>(3), 1, // lambda
239 *   FE_DGQ<dim>(0), 1), // f
240 *   dof_handler (triangulation)
241 *   {
242 *   velocity[0] = 100;
243 *  
244 * @endcode
245 *
246 * We have 50 detector points on an outer ring...
247 *
248 * @code
249 *   for (unsigned int i=0; i<50; ++i)
250 *   {
251 *   const Point<dim> p (0.6 * std::sin(2*numbers::PI * i/50),
252 *   0.6 * std::cos(2*numbers::PI * i/50));
253 *   detector_locations.push_back (p);
254 *   }
255 *  
256 * @endcode
257 *
258 * ...and another 50 detector points on an innner ring:
259 *
260 * @code
261 *   for (unsigned int i=0; i<50; ++i)
262 *   {
263 *   const Point<dim> p (0.2 * std::sin(2*numbers::PI * i/50),
264 *   0.2 * std::cos(2*numbers::PI * i/50));
265 *   detector_locations.push_back (p);
266 *   }
267 *  
268 * @endcode
269 *
270 * Generate the grid we will work on:
271 *
272 * @code
275 *  
276 * @endcode
277 *
278 * The detector locations are static, so we can already here
279 * generate a file that contains their locations. We use the
280 * particle framework to do this, using detector locations as
281 * particle locations.
282 *
283 * @code
284 *   {
287 *   for (const auto &loc : detector_locations)
288 *   {
289 *   Particles::Particle<dim> new_particle;
290 *   new_particle.set_location(loc);
291 * @endcode
292 *
293 * Insert the particle. It is a lie that the particle is in
294 * the first cell, but nothing we do actually cares about the
295 * cell a particle is in.
296 *
297 * @code
298 *   particle_handler.insert_particle(new_particle,
300 *   }
301 *  
302 *   Particles::DataOut<dim> particle_out;
303 *   particle_out.build_patches(particle_handler);
304 *   std::ofstream output("detector_locations.vtu");
305 *   particle_out.write_vtu(output);
306 *   }
307 *  
308 * @endcode
309 *
310 * While we're generating output, also output the source location. Do this
311 * by outputting many (1000) points that indicate the perimeter of the source
312 *
313 * @code
314 *   {
315 *   Particles::ParticleHandler<dim> particle_handler(triangulation,
316 *   StaticMappingQ1<dim>::mapping);
317 *  
318 *   const unsigned int n_points = 1000;
319 *   for (unsigned int i=0; i<n_points; ++i)
320 *   {
321 *   Point<dim> loc = source_location;
322 *   loc[0] += source_radius * std::cos(2*numbers::PI*i/n_points);
323 *   loc[1] += source_radius * std::sin(2*numbers::PI*i/n_points);
324 *  
325 *   Particles::Particle<dim> new_particle;
326 *   new_particle.set_location(loc);
327 *   particle_handler.insert_particle(new_particle,
328 *   triangulation.begin_active());
329 *   }
330 *  
331 *   Particles::DataOut<dim> particle_out;
332 *   particle_out.build_patches(particle_handler);
333 *   std::ofstream output("source_locations.vtu");
334 *   particle_out.write_vtu(output);
335 *   }
336 *   }
337 *  
338 *  
339 *  
340 * @endcode
341 *
342 * The following function solves a forward problem on a twice
343 * refined mesh to compute "synthetic data". Refining the mesh
344 * beyond the mesh used for the inverse problem avoids an
345 * inverse crime.
346 *
347 * @code
348 *   template <int dim>
349 *   void InformationDensityMeshRefinement<dim>::compute_synthetic_measurements ()
350 *   {
351 *   std::cout << "Computing synthetic data by solving the forward problem..."
352 *   << std::flush;
353 *  
354 * @endcode
355 *
356 * Create a triangulation and DoFHandler that corresponds to a
357 * twice-refined mesh so that we obtain the synthetic data with
358 * higher accuracy than we do on the regular mesh used for all other
359 * computations.
360 *
361 * @code
362 *   Triangulation<dim> forward_triangulation;
363 *   forward_triangulation.copy_triangulation (triangulation);
364 *   forward_triangulation.refine_global (2);
365 *  
366 *   const FE_Q<dim> forward_fe (fe.base_element(0).degree);
367 *   DoFHandler<dim> forward_dof_handler (forward_triangulation);
368 *   forward_dof_handler.distribute_dofs (forward_fe);
369 *  
370 *   AffineConstraints<double> constraints;
371 *   DoFTools::make_hanging_node_constraints(forward_dof_handler, constraints);
372 *   constraints.close();
373 *  
374 *   SparsityPattern sparsity (forward_dof_handler.n_dofs(),
375 *   forward_dof_handler.max_couplings_between_dofs());
376 *   DoFTools::make_sparsity_pattern (forward_dof_handler, sparsity);
377 *   constraints.condense (sparsity);
378 *   sparsity.compress ();
379 *  
380 *   SparseMatrix<double> system_matrix (sparsity);
381 *   Vector<double> system_rhs (forward_dof_handler.n_dofs());
382 *  
383 *   QGauss<dim> quadrature_formula(3);
384 *   FEValues<dim> fe_values (forward_fe, quadrature_formula,
385 *   update_values | update_gradients |
386 *   update_quadrature_points | update_JxW_values);
387 *  
388 *   const unsigned int dofs_per_cell = forward_fe.dofs_per_cell;
389 *   const unsigned int n_q_points = quadrature_formula.size();
390 *  
391 * @endcode
392 *
393 * First assemble the system matrix and right hand side for the forward
394 * problem:
395 *
396 * @code
397 *   {
398 *   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
399 *   Vector<double> cell_rhs (dofs_per_cell);
400 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
401 *  
402 *   for (const auto &cell : forward_dof_handler.active_cell_iterators())
403 *   {
404 *   fe_values.reinit (cell);
405 *   cell_matrix = 0;
406 *   cell_rhs = 0;
407 *  
408 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
409 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
410 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
411 *   cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
412 *   fe_values.shape_grad(j,q_point)
413 *   +
414 *   fe_values.shape_value(i,q_point) *
415 *   (velocity * fe_values.shape_grad(j,q_point))
416 *   )
417 *   *
418 *   fe_values.JxW(q_point);
419 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
420 *   if (fe_values.quadrature_point(q_point).distance (source_location)
421 *   < source_radius)
422 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
423 *   cell_rhs(i) +=
424 *   1.0 *
425 *   fe_values.shape_value (i, q_point) *
426 *   fe_values.JxW(q_point);
427 *  
428 *   cell->get_dof_indices (local_dof_indices);
429 *   constraints.distribute_local_to_global (cell_matrix,
430 *   cell_rhs,
431 *   local_dof_indices,
432 *   system_matrix,
433 *   system_rhs);
434 *   }
435 *  
436 *   std::map<unsigned int, double> boundary_values;
437 *   VectorTools::interpolate_boundary_values (forward_dof_handler,
438 *   0,
439 *   Functions::ZeroFunction<dim>(),
440 *   boundary_values);
441 *   Vector<double> tmp (forward_dof_handler.n_dofs());
442 *   MatrixTools::apply_boundary_values (boundary_values,
443 *   system_matrix,
444 *   tmp,
445 *   system_rhs);
446 *   }
447 *  
448 * @endcode
449 *
450 * Solve the forward problem and output it into its own VTU file :
451 *
452 * @code
453 *   SparseDirectUMFPACK A_inverse;
454 *   Vector<double> forward_solution (forward_dof_handler.n_dofs());
455 *   forward_solution = system_rhs;
456 *   A_inverse.solve(system_matrix, forward_solution);
457 *  
458 *   const double max_forward_solution = forward_solution.linfty_norm();
459 *  
460 *   {
461 *   DataOut<dim> data_out;
462 *   data_out.attach_dof_handler (forward_dof_handler);
463 *   data_out.add_data_vector (forward_solution, "c");
464 *   data_out.build_patches (4);
465 *  
466 *   std::ofstream out ("forward-solution.vtu");
467 *   data_out.write_vtu (out);
468 *   }
469 *  
470 * @endcode
471 *
472 * Now evaluate the forward solution at the measurement points:
473 *
474 * @code
475 *   for (const auto &p : detector_locations)
476 *   {
477 * @endcode
478 *
479 * same 10% noise level for all points
480 *
481 * @code
482 *   noise_level.push_back (0.1 * max_forward_solution);
483 *  
484 *   const double z_n = VectorTools::point_value(forward_dof_handler, forward_solution, p);
485 *   const double eps_n = Utilities::generate_normal_random_number(0, noise_level.back());
486 *  
487 *   measurement_values.push_back (z_n + eps_n);
488 *   }
489 *  
490 *   std::cout << std::endl;
491 *   }
492 *  
493 *  
494 * @endcode
495 *
496 * It will make our lives easier if we can always assume that detector
497 * locations are at cell centers, because then we can evaluate the
498 * solution there using a quadrature formula whose sole quadrature
499 * point lies at the center of a cell. That's of course not where the
500 * "real" detector locations are, but it does not introduce a large
501 * error to do this.
502 *
503 * @code
504 *   template <int dim>
505 *   void InformationDensityMeshRefinement<dim>::bounce_measurement_points_to_cell_centers ()
506 *   {
507 *   detector_locations_on_mesh = detector_locations;
508 *   for (auto &p : detector_locations_on_mesh)
509 *   {
510 *   for (const auto &cell : triangulation.active_cell_iterators())
511 *   if (cell->point_inside (p))
512 *   {
513 *   p = cell->center();
514 *   break;
515 *   }
516 *   }
517 *   }
518 *  
519 *  
520 * @endcode
521 *
522 * The following functions are all quite standard by what we have
523 * 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
524 * more typical programs):
525 *
526 * @code
527 *   template <int dim>
528 *   void InformationDensityMeshRefinement<dim>::setup_system ()
529 *   {
530 *   std::cout << "Setting up the linear system for the inverse problem..."
531 *   << std::endl;
532 *  
533 *   dof_handler.distribute_dofs (fe);
534 *   DoFRenumbering::component_wise (dof_handler);
535 *  
536 *   hanging_node_constraints.clear ();
538 *   hanging_node_constraints);
539 *   hanging_node_constraints.close();
540 *  
541 *   std::cout << " Number of active cells: "
543 *   << std::endl;
544 *   std::cout << " Number of degrees of freedom: "
545 *   << dof_handler.n_dofs()
546 *   << std::endl;
547 *  
548 *   const std::vector<types::global_dof_index> dofs_per_component =
550 *   BlockDynamicSparsityPattern c_sparsity(dofs_per_component,dofs_per_component);
551 *   DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
552 *   hanging_node_constraints.condense(c_sparsity);
553 *   sparsity_pattern.copy_from(c_sparsity);
554 *  
555 *   system_matrix.reinit (sparsity_pattern);
556 *  
557 *   solution.reinit (dofs_per_component);
558 *   system_rhs.reinit (dofs_per_component);
559 *   }
560 *  
561 *  
562 *  
563 *   template <int dim>
564 *   void InformationDensityMeshRefinement<dim>::assemble_system ()
565 *   {
566 *   std::cout << "Assembling the linear system for the inverse problem..."
567 *   << std::flush;
568 *  
569 *   QGauss<dim> quadrature_formula(3);
570 *  
571 *   FEValues<dim> fe_values (fe, quadrature_formula,
574 *  
575 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
576 *   const unsigned int n_q_points = quadrature_formula.size();
577 *  
578 *   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
579 *   Vector<double> cell_rhs (dofs_per_cell);
580 *  
581 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
582 *  
583 *   const FEValuesExtractors::Scalar c(0), lambda(1), f(2);
584 *  
585 *   for (const auto &cell : dof_handler.active_cell_iterators())
586 *   {
587 *   fe_values.reinit (cell);
588 *   cell_matrix = 0;
589 *   cell_rhs = 0;
590 *  
591 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
592 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
593 *   {
594 *   const Tensor<1,dim> grad_phi_i = fe_values[c].gradient (i,q_point);
595 *   const Tensor<1,dim> grad_psi_i = fe_values[lambda].gradient (i,q_point);
596 *  
597 *   const double phi_i = fe_values[c].value (i,q_point);
598 *   const double psi_i = fe_values[lambda].value (i,q_point);
599 *   const double chi_i = fe_values[f].value (i,q_point);
600 *  
601 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
602 *   {
603 *   const Tensor<1,dim> grad_phi_j = fe_values[c].gradient (j,q_point);
604 *   const Tensor<1,dim> grad_psi_j = fe_values[lambda].gradient (j,q_point);
605 *  
606 *   const double phi_j = fe_values[c].value (j,q_point);
607 *   const double psi_j= fe_values[lambda].value (j,q_point);
608 *   const double chi_j = fe_values[f].value (j,q_point);
609 *  
610 *   cell_matrix(i,j) +=
611 *   ((grad_phi_i * grad_phi_j
612 *   +
613 *   phi_i * (velocity * grad_phi_j)
614 *   -
615 *   phi_i * chi_j
616 *   +
617 *   grad_psi_i * grad_psi_j
618 *   -
619 *   psi_i * (velocity * grad_psi_j)
620 *   -
621 *   chi_i * psi_j
622 *   +
623 *   regularization_parameter * chi_i * chi_j
624 *   ) *
625 *   fe_values.JxW (q_point));
626 *  
627 *   for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
628 *   if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
629 *   {
630 *   cell_matrix(i,j) += psi_i * phi_j / noise_level[n] / noise_level[n];
631 *   }
632 *   }
633 *  
634 *   for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
635 *   if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
636 *   cell_rhs(i) += psi_i * measurement_values[n] / noise_level[n] / noise_level[n];
637 *   }
638 *  
639 *   cell->get_dof_indices (local_dof_indices);
640 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
641 *   {
642 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
643 *   system_matrix.add (local_dof_indices[i],
644 *   local_dof_indices[j],
645 *   cell_matrix(i,j));
646 *  
647 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
648 *   }
649 *   }
650 *  
651 *   hanging_node_constraints.condense (system_matrix);
652 *   hanging_node_constraints.condense (system_rhs);
653 *  
654 *   std::map<unsigned int,double> boundary_values;
655 *   std::vector<bool> component_mask (3);
656 *   component_mask[0] = component_mask[1] = true;
657 *   component_mask[2] = false;
659 *   0,
661 *   boundary_values,
662 *   component_mask);
663 *   MatrixTools::apply_boundary_values (boundary_values,
664 *   system_matrix,
665 *   solution,
666 *   system_rhs);
667 *  
668 *   std::cout << std::endl;
669 *   }
670 *  
671 *  
672 *  
673 *   template <int dim>
674 *   void InformationDensityMeshRefinement<dim>::solve ()
675 *   {
676 *   std::cout << "Solving the linear system for the inverse problem..."
677 *   << std::flush;
678 *  
679 *   SparseDirectUMFPACK A_direct;
680 *   solution = system_rhs;
681 *   A_direct.solve(system_matrix, solution);
682 *  
683 *   hanging_node_constraints.distribute (solution);
684 *  
685 *   std::cout << std::endl;
686 *   }
687 *  
688 *  
689 *  
690 * @endcode
691 *
692 * This is really the only interesting function of this program. It
693 * computes the functions @f$h_K = A^{-1} s_K@f$ for each source function
694 * (corresponding to each cell of the mesh). To do so, it first
695 * computes the forward matrix @f$A@f$ and uses the SparseDirectUMFPACK
696 * class to build an LU decomposition for this matrix. Then it loops
697 * over all cells @f$K@f$ and computes the corresponding @f$h_K@f$ by applying
698 * the LU decomposition to a right hand side vector for each @f$s_K@f$.
699 *
700
701 *
702 * The actual information content is then computed by evaluating these
703 * functions @f$h_K@f$ at measurement locations.
704 *
705 * @code
706 *   template <int dim>
707 *   void InformationDensityMeshRefinement<dim>::compute_information_content ()
708 *   {
709 *   std::cout << "Computing the information content..."
710 *   << std::flush;
711 *  
712 *   information_content.reinit (triangulation.n_active_cells());
713 *  
714 *   const FE_Q<dim> information_fe (fe.base_element(0).degree);
715 *   DoFHandler<dim> information_dof_handler (triangulation);
716 *   information_dof_handler.distribute_dofs (information_fe);
717 *  
718 *   AffineConstraints<double> constraints;
719 *   DoFTools::make_hanging_node_constraints(information_dof_handler, constraints);
720 *   constraints.close();
721 *  
722 *   SparsityPattern sparsity (information_dof_handler.n_dofs(),
723 *   information_dof_handler.max_couplings_between_dofs());
724 *   DoFTools::make_sparsity_pattern (information_dof_handler, sparsity);
725 *   constraints.condense (sparsity);
726 *   sparsity.compress ();
727 *  
728 *   SparseMatrix<double> system_matrix (sparsity);
729 *  
730 *   QGauss<dim> quadrature_formula(3);
731 *  
732 *   const unsigned int dofs_per_cell = information_fe.dofs_per_cell;
733 *   const unsigned int n_q_points = quadrature_formula.size();
734 *  
735 * @endcode
736 *
737 * First build the forward operator
738 *
739 * @code
740 *   {
741 *   FEValues<dim> fe_values (information_fe, quadrature_formula,
744 *  
745 *   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
746 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
747 *  
748 *   for (const auto &cell : information_dof_handler.active_cell_iterators())
749 *   {
750 *   fe_values.reinit (cell);
751 *   cell_matrix = 0;
752 *  
753 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
754 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
755 *   for (unsigned int j=0; j<dofs_per_cell; ++j)
756 *   cell_matrix(i,j) += (fe_values.shape_grad (i,q_point) *
757 *   fe_values.shape_grad(j,q_point)
758 *   +
759 *   fe_values.shape_value(i,q_point) *
760 *   (velocity * fe_values.shape_grad(j,q_point))) *
761 *   fe_values.JxW(q_point);
762 *  
763 *   cell->distribute_local_to_global (cell_matrix,
764 *   system_matrix);
765 *   }
766 *  
767 *   constraints.condense (system_matrix);
768 *  
769 *   std::map<unsigned int, double> boundary_values;
770 *   VectorTools::interpolate_boundary_values (information_dof_handler,
771 *   0,
773 *   boundary_values);
774 *   Vector<double> tmp (information_dof_handler.n_dofs());
775 *   MatrixTools::apply_boundary_values (boundary_values,
776 *   system_matrix,
777 *   tmp,
778 *   tmp);
779 *   }
780 *  
781 * @endcode
782 *
783 * Then factorize
784 *
785 * @code
786 *   SparseDirectUMFPACK A_inverse;
787 *   A_inverse.factorize(system_matrix);
788 *  
789 * @endcode
790 *
791 * Now compute the solutions corresponding to the possible
792 * sources. Each source is active on exactly one cell.
793 *
794
795 *
796 * As mentioned in the paper, this is a trivially parallel job, so
797 * we send the computations for each of these cells onto a separate
798 * task and let the OS schedule them onto individual processor
799 * cores.
800 *
801 * @code
802 *   Threads::TaskGroup<void> tasks;
803 *   for (unsigned int K=0; K<triangulation.n_active_cells(); ++K)
804 *   tasks +=
805 *   Threads::new_task([&,K]()
806 *   {
807 *   Vector<double> rhs (information_dof_handler.n_dofs());
808 *   Vector<double> cell_rhs (dofs_per_cell);
809 *   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
810 *  
812 *   cell = information_dof_handler.begin_active();
813 *   std::advance (cell, K);
814 *  
815 *   FEValues<dim> fe_values (information_fe, quadrature_formula,
816 *   update_values |
818 *  
819 *   fe_values.reinit (cell);
820 *   cell_rhs = 0;
821 *  
822 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
823 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
824 *   cell_rhs(i) += fe_values.shape_value (i,q_point) *
825 *   fe_values.JxW(q_point);
826 *  
827 *   cell->distribute_local_to_global (cell_rhs,
828 *   rhs);
829 *  
830 *   constraints.condense (rhs);
831 *  
832 *   A_inverse.solve(rhs);
833 *  
834 *   constraints.distribute (rhs);
835 *  
836 * @endcode
837 *
838 * Having computed the forward solutions
839 * corresponding to this source term, evaluate its
840 * contribution to the information content on all
841 * cells of the mesh by taking into account the
842 * detector locations. We add these into global
843 * objects, so we have to guard access to the
844 * global object:
845 *
846 * @code
847 *   static std::mutex m;
848 *   std::lock_guard<std::mutex> g(m);
849 *  
850 *  
851 *   information_content(K) = regularization_parameter * cell->measure() * cell->measure();
852 *   std::vector<double> local_h_K_values (n_q_points);
853 *   for (const auto &cell : information_dof_handler.active_cell_iterators())
854 *   {
855 *   fe_values.reinit (cell);
856 *   fe_values.get_function_values (rhs, local_h_K_values);
857 *  
858 *   for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
859 *   for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n)
860 *   if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12)
861 *   information_content(K) += local_h_K_values[q_point]
862 *   * local_h_K_values[q_point]
863 *   / noise_level[n] / noise_level[n];
864 *   }
865 *   }
866 *   );
867 *  
868 * @endcode
869 *
870 * And wait:
871 *
872 * @code
873 *   tasks.join_all();
874 *  
875 *   std::cout << std::endl;
876 *   }
877 *  
878 *  
879 *  
880 * @endcode
881 *
882 * Create graphical output for all of the principal variables of the
883 * problem (c,lambda,f) as well as for the information content and
884 * density. Then also output the various blocks of the matrix so we
885 * can compute the eigenvalues of the H matrix mentioned in the paper.
886 *
887 * @code
888 *   template <int dim>
889 *   void InformationDensityMeshRefinement<dim>::output_results (const unsigned int cycle) const
890 *   {
891 *   std::cout << "Outputting solutions..." << std::flush;
892 *  
893 *   DataOut<dim> data_out;
894 *  
895 *   std::vector<std::string> names;
896 *   names.push_back ("forward_solution");
897 *   names.push_back ("adjoint_solution");
898 *   names.push_back ("recovered_parameter");
899 *  
900 *   data_out.attach_dof_handler (dof_handler);
901 *   data_out.add_data_vector (solution, names);
902 *   data_out.add_data_vector (information_content, "information_content");
903 *  
904 *   Vector<double> information_density (triangulation.n_active_cells());
905 *   for (const auto &cell : triangulation.active_cell_iterators())
906 *   information_density(cell->active_cell_index())
907 *   = std::sqrt(information_content(cell->active_cell_index())) / cell->measure();
908 *   data_out.add_data_vector (information_density, "information_density");
909 *  
910 *   data_out.build_patches ();
911 *  
912 *   std::string filename = "solution-";
913 *   filename += ('0'+cycle);
914 *   filename += ".vtu";
915 *  
916 *   std::ofstream output (filename.c_str());
917 *   data_out.write_vtu (output);
918 *  
919 *  
920 * @endcode
921 *
922 * Now output the individual blocks of the matrix into files.
923 *
924 * @code
925 *   auto write_block = [&](const unsigned int block_i,
926 *   const unsigned int block_j,
927 *   const std::string &filename)
928 *   {
929 *   std::ofstream o(filename);
930 *   system_matrix.block(block_i,block_j).print (o);
931 *   };
932 *   write_block(0,0, "matrix-" + std::to_string(cycle) + "-A.txt");
933 *   write_block(0,2, "matrix-" + std::to_string(cycle) + "-B.txt");
934 *   write_block(1,0, "matrix-" + std::to_string(cycle) + "-C.txt");
935 *   write_block(2,2, "matrix-" + std::to_string(cycle) + "-M.txt");
936 *  
937 *   std::cout << std::endl;
938 *   }
939 *  
940 *  
941 *  
942 * @endcode
943 *
944 * The following is then a function that refines the mesh based on the
945 * refinement criteria described in the paper. Which criterion to use
946 * is determined by which value the `refinement_criterion` variable
947 * is set to.
948 *
949 * @code
950 *   template <int dim>
951 *   void InformationDensityMeshRefinement<dim>::refine_grid ()
952 *   {
953 *   std::cout << "Refining the mesh..." << std::endl;
954 *  
955 *   enum RefinementCriterion
956 *   {
957 *   global,
958 *   information_content,
959 *   indicator,
960 *   smoothness
961 *   };
962 *   const RefinementCriterion refinement_criterion = information_content;
963 *  
964 *   switch (refinement_criterion)
965 *   {
966 *   case global:
967 *   {
969 *   break;
970 *   }
971 *  
972 *   case information_content:
973 *   {
975 *   this->information_content,
976 *   0.2, 0.05);
978 *   break;
979 *   }
980 *  
981 *   case indicator:
982 *   {
983 *   Vector<double> refinement_indicators (triangulation.n_active_cells());
984 *  
985 *   QGauss<dim> quadrature(3);
986 *   FEValues<dim> fe_values (fe, quadrature, update_values | update_JxW_values);
987 *  
988 *   const FEValuesExtractors::Scalar lambda(1), f(2);
989 *  
990 *   std::vector<double> lambda_values (quadrature.size());
991 *   std::vector<double> f_values (quadrature.size());
992 *  
993 *   for (const auto &cell : dof_handler.active_cell_iterators())
994 *   {
995 *   fe_values.reinit (cell);
996 *   fe_values[lambda].get_function_values (solution, lambda_values);
997 *   fe_values[f].get_function_values (solution, f_values);
998 *  
999 *   for (unsigned int q=0; q<quadrature.size(); ++q)
1000 *   refinement_indicators(cell->active_cell_index())
1001 *   += (std::fabs (regularization_parameter * f_values[q]
1002 *   -
1003 *   lambda_values[q])
1004 *   * fe_values.JxW(q));
1005 *   }
1006 *  
1008 *   refinement_indicators,
1009 *   0.2, 0.05);
1011 *   break;
1012 *   }
1013 *  
1014 *  
1015 *   case smoothness:
1016 *   {
1017 *   Vector<float> refinement_indicators (triangulation.n_active_cells());
1018 *  
1020 *   solution,
1021 *   refinement_indicators,
1022 *   /*component=*/2);
1023 * @endcode
1024 *
1025 * and scale it to obtain an error indicator.
1026 *
1027 * @code
1028 *   for (const auto &cell : triangulation.active_cell_iterators())
1029 *   refinement_indicators[cell->active_cell_index()] *=
1030 *   std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
1031 *  
1032 *  
1034 *   refinement_indicators,
1035 *   0.2, 0.05);
1037 *   break;
1038 *   }
1039 *  
1040 *   default:
1041 *   Assert (false, ExcInternalError());
1042 *   }
1043 *  
1044 *   bounce_measurement_points_to_cell_centers ();
1045 *  
1046 *  
1047 *   std::cout << std::endl;
1048 *   }
1049 *  
1050 *  
1051 *  
1052 *  
1053 *   template <int dim>
1054 *   void InformationDensityMeshRefinement<dim>::run ()
1055 *   {
1056 *   std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
1057 *  
1058 *   compute_synthetic_measurements ();
1059 *   bounce_measurement_points_to_cell_centers ();
1060 *  
1061 *   for (unsigned int cycle=0; cycle<7; ++cycle)
1062 *   {
1063 *   std::cout << "---------- Cycle " << cycle << " ------------" << std::endl;
1064 *  
1065 *   setup_system ();
1066 *   assemble_system ();
1067 *   solve ();
1068 *   compute_information_content ();
1069 *   output_results (cycle);
1070 *   refine_grid ();
1071 *   }
1072 *   }
1073 *  
1074 *  
1075 *  
1076 *   int main ()
1077 *   {
1078 *   try
1079 *   {
1080 *   deallog.depth_console (0);
1081 *  
1082 *   InformationDensityMeshRefinement<2> information_density_mesh_refinement;
1083 *   information_density_mesh_refinement.run ();
1084 *   }
1085 *   catch (std::exception &exc)
1086 *   {
1087 *   std::cerr << std::endl << std::endl
1088 *   << "----------------------------------------------------"
1089 *   << std::endl;
1090 *   std::cerr << "Exception on processing: " << std::endl
1091 *   << exc.what() << std::endl
1092 *   << "Aborting!" << std::endl
1093 *   << "----------------------------------------------------"
1094 *   << std::endl;
1095 *  
1096 *   return 1;
1097 *   }
1098 *   catch (...)
1099 *   {
1100 *   std::cerr << std::endl << std::endl
1101 *   << "----------------------------------------------------"
1102 *   << std::endl;
1103 *   std::cerr << "Unknown exception!" << std::endl
1104 *   << "Aborting!" << std::endl
1105 *   << "----------------------------------------------------"
1106 *   << std::endl;
1107 *   return 1;
1108 *   }
1109 *  
1110 *   return 0;
1111 *   }
1112 * @endcode
1113
1114
1115*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Definition fe_q.h:554
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:349
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:27
Definition point.h:111
void solve(Vector< double > &rhs_and_solution, const bool transpose=false) const
void factorize(const Matrix &matrix)
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
active_cell_iterator begin_active(const unsigned int level=0) const
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
Point< 3 > center
Point< 2 > first
Definition grid_out.cc:4623
__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={}, 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:36
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)
@ 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:74
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={})
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)