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
cdr.h
Go to the documentation of this file.
1
127 *  
128 *   #ifndef dealii__cdr_parameters_h
129 *   #define dealii__cdr_parameters_h
130 *  
131 *   #include <deal.II/base/parameter_handler.h>
132 *  
133 *   #include <string>
134 *  
135 * @endcode
136 *
137 * I prefer to use the ParameterHandler class in a slightly different way than
138 * usual: The class Parameters creates, uses, and then destroys a
139 * ParameterHandler inside the <code>read_parameter_file</code> method instead
140 * of keeping it around. This is nice because now all of the run time
141 * parameters are contained in a simple class and it can be copied or passed
142 * around very easily.
143 *
144 * @code
145 *   namespace CDR
146 *   {
147 *   using namespace dealii;
148 *  
149 *   class Parameters
150 *   {
151 *   public:
152 *   double inner_radius;
153 *   double outer_radius;
154 *  
155 *   double diffusion_coefficient;
156 *   double reaction_coefficient;
157 *   bool time_dependent_forcing;
158 *  
159 *   unsigned int initial_refinement_level;
160 *   unsigned int max_refinement_level;
161 *   unsigned int fe_order;
162 *  
163 *   double start_time;
164 *   double stop_time;
165 *   unsigned int n_time_steps;
166 *  
167 *   unsigned int save_interval;
168 *   unsigned int patch_level;
169 *  
170 *   void
171 *   read_parameter_file(const std::string &file_name);
172 *  
173 *   private:
174 *   void
175 *   configure_parameter_handler(ParameterHandler &parameter_handler);
176 *   };
177 *   } // namespace CDR
178 *   #endif
179 * @endcode
180
181
182<a name="ann-common/include/deal.II-cdr/system_matrix.h"></a>
183<h1>Annotated version of common/include/deal.II-cdr/system_matrix.h</h1>
184 *
185 *
186 *
187 *
188 * @code
189 *   /* -----------------------------------------------------------------------------
190 *   *
191 *   * SPDX-License-Identifier: LGPL-2.1-or-later
192 *   * Copyright (C) 2015 by David Wells
193 *   *
194 *   * This file is part of the deal.II code gallery.
195 *   *
196 *   * -----------------------------------------------------------------------------
197 *   */
198 *  
199 *   #ifndef dealii__cdr_system_matrix_h
200 *   #define dealii__cdr_system_matrix_h
201 *   #include <deal.II/base/quadrature_lib.h>
202 *  
203 *   #include <deal.II/dofs/dof_handler.h>
204 *  
205 *   #include <deal.II/lac/affine_constraints.h>
206 *  
207 *   #include <deal.II-cdr/parameters.h>
208 *  
209 *   #include <functional>
210 *  
211 * @endcode
212 *
213 * One of the goals I had in writing this entry was to split up functions into
214 * different compilation units instead of using one large file. This is the
215 * header file for a pair of functions (only one of which I ultimately use)
216 * which build the system matrix.
217 *
218 * @code
219 *   namespace CDR
220 *   {
221 *   using namespace dealii;
222 *  
223 *   template <int dim, typename MatrixType>
224 *   void
225 *   create_system_matrix(
226 *   const DoFHandler<dim> & dof_handler,
227 *   const QGauss<dim> & quad,
228 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
229 *   const CDR::Parameters & parameters,
230 *   const double time_step,
231 *   MatrixType & system_matrix);
232 *  
233 *   template <int dim, typename MatrixType>
234 *   void
235 *   create_system_matrix(
236 *   const DoFHandler<dim> & dof_handler,
237 *   const QGauss<dim> & quad,
238 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
239 *   const CDR::Parameters & parameters,
240 *   const double time_step,
241 *   const AffineConstraints<double> & constraints,
242 *   MatrixType & system_matrix);
243 *   } // namespace CDR
244 *   #endif
245 * @endcode
246
247
248<a name="ann-common/include/deal.II-cdr/system_matrix.templates.h"></a>
249<h1>Annotated version of common/include/deal.II-cdr/system_matrix.templates.h</h1>
250 *
251 *
252 *
253 *
254 * @code
255 *   /* -----------------------------------------------------------------------------
256 *   *
257 *   * SPDX-License-Identifier: LGPL-2.1-or-later
258 *   * Copyright (C) 2015 by David Wells
259 *   *
260 *   * This file is part of the deal.II code gallery.
261 *   *
262 *   * -----------------------------------------------------------------------------
263 *   */
264 *  
265 *   #ifndef dealii__cdr_system_matrix_templates_h
266 *   #define dealii__cdr_system_matrix_templates_h
267 *   #include <deal.II/base/quadrature_lib.h>
268 *  
269 *   #include <deal.II/dofs/dof_handler.h>
270 *  
271 *   #include <deal.II/fe/fe_q.h>
272 *   #include <deal.II/fe/fe_values.h>
273 *  
274 *   #include <deal.II/lac/affine_constraints.h>
275 *  
276 *   #include <deal.II-cdr/parameters.h>
277 *   #include <deal.II-cdr/system_matrix.h>
278 *  
279 *   #include <functional>
280 *   #include <vector>
281 *  
282 *   namespace CDR
283 *   {
284 *   using namespace dealii;
285 *  
286 * @endcode
287 *
288 * This is the actual implementation of the <code>create_system_matrix</code>
289 * function described in the header file. It is similar to the system matrix
290 * assembly routine in @ref step_40 "step-40".
291 *
292 * @code
293 *   template <int dim, typename UpdateFunction>
294 *   void
295 *   internal_create_system_matrix(
296 *   const DoFHandler<dim> & dof_handler,
297 *   const QGauss<dim> & quad,
298 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
299 *   const CDR::Parameters & parameters,
300 *   const double time_step,
301 *   UpdateFunction update_system_matrix)
302 *   {
303 *   auto & fe = dof_handler.get_fe();
304 *   const auto dofs_per_cell = fe.dofs_per_cell;
305 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
306 *   FEValues<dim> fe_values(fe,
307 *   quad,
308 *   update_values | update_gradients |
309 *   update_quadrature_points | update_JxW_values);
310 *  
311 *   std::vector<types::global_dof_index> local_indices(dofs_per_cell);
312 *  
313 *   for (const auto &cell : dof_handler.active_cell_iterators())
314 *   {
315 *   if (cell->is_locally_owned())
316 *   {
317 *   fe_values.reinit(cell);
318 *   cell_matrix = 0.0;
319 *   cell->get_dof_indices(local_indices);
320 *   for (unsigned int q = 0; q < quad.size(); ++q)
321 *   {
322 *   const auto current_convection =
323 *   convection_function(fe_values.quadrature_point(q));
324 *  
325 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
326 *   {
327 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
328 *   {
329 *   const auto convection_contribution =
330 *   current_convection * fe_values.shape_grad(j, q);
331 *   cell_matrix(i, j) +=
332 *   fe_values.JxW(q) *
333 * @endcode
334 *
335 * Here are the time step, mass, and reaction parts:
336 *
337 * @code
338 *   ((1.0 +
339 *   time_step / 2.0 * parameters.reaction_coefficient) *
340 *   fe_values.shape_value(i, q) *
341 *   fe_values.shape_value(j, q) +
342 *   time_step / 2.0 *
343 * @endcode
344 *
345 * and the convection part:
346 *
347 * @code
348 *   (fe_values.shape_value(i, q) *
349 *   convection_contribution
350 * @endcode
351 *
352 * and, finally, the diffusion part:
353 *
354 * @code
355 *   + parameters.diffusion_coefficient *
356 *   (fe_values.shape_grad(i, q) *
357 *   fe_values.shape_grad(j, q))));
358 *   }
359 *   }
360 *   }
361 *   update_system_matrix(local_indices, cell_matrix);
362 *   }
363 *   }
364 *   }
365 *  
366 *   template <int dim, typename MatrixType>
367 *   void
368 *   create_system_matrix(
369 *   const DoFHandler<dim> & dof_handler,
370 *   const QGauss<dim> & quad,
371 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
372 *   const CDR::Parameters & parameters,
373 *   const double time_step,
374 *   const AffineConstraints<double> & constraints,
375 *   MatrixType & system_matrix)
376 *   {
377 *   internal_create_system_matrix<dim>(
378 *   dof_handler,
379 *   quad,
380 *   convection_function,
381 *   parameters,
382 *   time_step,
383 *   [&constraints, &system_matrix](
384 *   const std::vector<types::global_dof_index> &local_indices,
385 *   const FullMatrix<double> & cell_matrix) {
386 *   constraints.distribute_local_to_global(cell_matrix,
387 *   local_indices,
388 *   system_matrix);
389 *   });
390 *   }
391 *  
392 *   template <int dim, typename MatrixType>
393 *   void
394 *   create_system_matrix(
395 *   const DoFHandler<dim> & dof_handler,
396 *   const QGauss<dim> & quad,
397 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
398 *   const CDR::Parameters & parameters,
399 *   const double time_step,
400 *   MatrixType & system_matrix)
401 *   {
402 *   internal_create_system_matrix<dim>(
403 *   dof_handler,
404 *   quad,
405 *   convection_function,
406 *   parameters,
407 *   time_step,
408 *   [&system_matrix](
409 *   const std::vector<types::global_dof_index> &local_indices,
410 *   const FullMatrix<double> & cell_matrix) {
411 *   system_matrix.add(local_indices, cell_matrix);
412 *   });
413 *   }
414 *   } // namespace CDR
415 *   #endif
416 * @endcode
417
418
419<a name="ann-common/include/deal.II-cdr/system_rhs.h"></a>
420<h1>Annotated version of common/include/deal.II-cdr/system_rhs.h</h1>
421 *
422 *
423 *
424 *
425 * @code
426 *   /* -----------------------------------------------------------------------------
427 *   *
428 *   * SPDX-License-Identifier: LGPL-2.1-or-later
429 *   * Copyright (C) 2015 by David Wells
430 *   *
431 *   * This file is part of the deal.II code gallery.
432 *   *
433 *   * -----------------------------------------------------------------------------
434 *   */
435 *  
436 *   #ifndef dealii__cdr_system_rhs_h
437 *   #define dealii__cdr_system_rhs_h
438 *   #include <deal.II/base/point.h>
439 *   #include <deal.II/base/quadrature_lib.h>
440 *   #include <deal.II/base/tensor.h>
441 *  
442 *   #include <deal.II/dofs/dof_handler.h>
443 *  
444 *   #include <deal.II/lac/affine_constraints.h>
445 *  
446 *   #include <deal.II-cdr/parameters.h>
447 *  
448 *   #include <functional>
449 *  
450 * @endcode
451 *
452 * Similarly to <code>create_system_matrix</code>, I wrote a separate function
453 * to compute the right hand side.
454 *
455 * @code
456 *   namespace CDR
457 *   {
458 *   using namespace dealii;
459 *  
460 *   template <int dim, typename VectorType>
461 *   void
462 *   create_system_rhs(
463 *   const DoFHandler<dim> & dof_handler,
464 *   const QGauss<dim> & quad,
465 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
466 *   const std::function<double(double, const Point<dim>)> &forcing_function,
467 *   const CDR::Parameters & parameters,
468 *   const VectorType & previous_solution,
469 *   const AffineConstraints<double> & constraints,
470 *   const double current_time,
471 *   VectorType & system_rhs);
472 *   } // namespace CDR
473 *   #endif
474 * @endcode
475
476
477<a name="ann-common/include/deal.II-cdr/system_rhs.templates.h"></a>
478<h1>Annotated version of common/include/deal.II-cdr/system_rhs.templates.h</h1>
479 *
480 *
481 *
482 *
483 * @code
484 *   /* -----------------------------------------------------------------------------
485 *   *
486 *   * SPDX-License-Identifier: LGPL-2.1-or-later
487 *   * Copyright (C) 2015 by David Wells
488 *   *
489 *   * This file is part of the deal.II code gallery.
490 *   *
491 *   * -----------------------------------------------------------------------------
492 *   */
493 *  
494 *   #ifndef dealii__cdr_system_rhs_templates_h
495 *   #define dealii__cdr_system_rhs_templates_h
496 *   #include <deal.II/base/point.h>
497 *   #include <deal.II/base/quadrature_lib.h>
498 *   #include <deal.II/base/tensor.h>
499 *  
500 *   #include <deal.II/dofs/dof_handler.h>
501 *  
502 *   #include <deal.II/fe/fe_q.h>
503 *   #include <deal.II/fe/fe_values.h>
504 *  
505 *   #include <deal.II/lac/affine_constraints.h>
506 *   #include <deal.II/lac/vector.h>
507 *  
508 *   #include <deal.II-cdr/parameters.h>
509 *   #include <deal.II-cdr/system_rhs.h>
510 *  
511 *   #include <functional>
512 *   #include <vector>
513 *  
514 *   namespace CDR
515 *   {
516 *   using namespace dealii;
517 *  
518 *   template <int dim, typename VectorType>
519 *   void
520 *   create_system_rhs(
521 *   const DoFHandler<dim> & dof_handler,
522 *   const QGauss<dim> & quad,
523 *   const std::function<Tensor<1, dim>(const Point<dim>)> &convection_function,
524 *   const std::function<double(double, const Point<dim>)> &forcing_function,
525 *   const CDR::Parameters & parameters,
526 *   const VectorType & previous_solution,
527 *   const AffineConstraints<double> & constraints,
528 *   const double current_time,
529 *   VectorType & system_rhs)
530 *   {
531 *   auto & fe = dof_handler.get_fe();
532 *   const auto dofs_per_cell = fe.dofs_per_cell;
533 *   const double time_step =
534 *   (parameters.stop_time - parameters.start_time) / parameters.n_time_steps;
535 *   FEValues<dim> fe_values(fe,
536 *   quad,
537 *   update_values | update_gradients |
538 *   update_quadrature_points | update_JxW_values);
539 *  
540 *   Vector<double> cell_rhs(dofs_per_cell);
541 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
542 *  
543 *   Vector<double> current_fe_coefficients(dofs_per_cell);
544 *   std::vector<types::global_dof_index> local_indices(dofs_per_cell);
545 *  
546 *   const double previous_time{current_time - time_step};
547 *  
548 *   for (const auto &cell : dof_handler.active_cell_iterators())
549 *   {
550 *   if (cell->is_locally_owned())
551 *   {
552 *   fe_values.reinit(cell);
553 *   cell_rhs = 0.0;
554 *   cell->get_dof_indices(local_indices);
555 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
556 *   {
557 *   current_fe_coefficients[i] =
558 *   previous_solution[local_indices[i]];
559 *   }
560 *  
561 *   for (unsigned int q = 0; q < quad.size(); ++q)
562 *   {
563 *   const auto current_convection =
564 *   convection_function(fe_values.quadrature_point(q));
565 *  
566 *   const double current_forcing =
567 *   forcing_function(current_time, fe_values.quadrature_point(q));
568 *   const double previous_forcing =
569 *   forcing_function(previous_time,
570 *   fe_values.quadrature_point(q));
571 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
572 *   {
573 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
574 *   {
575 *   const auto convection_contribution =
576 *   current_convection * fe_values.shape_grad(j, q);
577 *  
578 *   cell_rhs(i) +=
579 *   fe_values.JxW(q) *
580 * @endcode
581 *
582 * Here are the mass and reaction part:
583 *
584 * @code
585 *   (((1.0 - time_step / 2.0 *
586 *   parameters.reaction_coefficient) *
587 *   fe_values.shape_value(i, q) *
588 *   fe_values.shape_value(j, q) -
589 *   time_step / 2.0 *
590 * @endcode
591 *
592 * the convection part:
593 *
594 * @code
595 *   (fe_values.shape_value(i, q) *
596 *   convection_contribution
597 * @endcode
598 *
599 * the diffusion part:
600 *
601 * @code
602 *   + parameters.diffusion_coefficient *
603 *   (fe_values.shape_grad(i, q) *
604 *   fe_values.shape_grad(j, q)))) *
605 *   current_fe_coefficients[j]
606 * @endcode
607 *
608 * and, finally, the forcing function part:
609 *
610 * @code
611 *   + time_step / 2.0 *
612 *   (current_forcing + previous_forcing) *
613 *   fe_values.shape_value(i, q));
614 *   }
615 *   }
616 *   }
617 *   constraints.distribute_local_to_global(cell_rhs,
618 *   local_indices,
619 *   system_rhs);
620 *   }
621 *   }
622 *   }
623 *   } // namespace CDR
624 *   #endif
625 * @endcode
626
627
628<a name="ann-common/include/deal.II-cdr/write_pvtu_output.h"></a>
629<h1>Annotated version of common/include/deal.II-cdr/write_pvtu_output.h</h1>
630 *
631 *
632 *
633 *
634 * @code
635 *   /* -----------------------------------------------------------------------------
636 *   *
637 *   * SPDX-License-Identifier: LGPL-2.1-or-later
638 *   * Copyright (C) 2015 by David Wells
639 *   *
640 *   * This file is part of the deal.II code gallery.
641 *   *
642 *   * -----------------------------------------------------------------------------
643 *   */
644 *  
645 *   #ifndef dealii__cdr_write_pvtu_output_h
646 *   #define dealii__cdr_write_pvtu_output_h
647 *   #include <deal.II/dofs/dof_handler.h>
648 *  
649 * @endcode
650 *
651 * This is a small class which handles PVTU output.
652 *
653 * @code
654 *   namespace CDR
655 *   {
656 *   using namespace dealii;
657 *  
658 *   class WritePVTUOutput
659 *   {
660 *   public:
661 *   WritePVTUOutput(const unsigned int patch_level);
662 *  
663 *   template <int dim, typename VectorType>
664 *   void
665 *   write_output(const DoFHandler<dim> &dof_handler,
666 *   const VectorType & solution,
667 *   const unsigned int time_step_n,
668 *   const double current_time);
669 *  
670 *   private:
671 *   const unsigned int patch_level;
672 *   const unsigned int this_mpi_process;
673 *   };
674 *   } // namespace CDR
675 *   #endif
676 * @endcode
677
678
679<a name="ann-common/include/deal.II-cdr/write_pvtu_output.templates.h"></a>
680<h1>Annotated version of common/include/deal.II-cdr/write_pvtu_output.templates.h</h1>
681 *
682 *
683 *
684 *
685 * @code
686 *   /* -----------------------------------------------------------------------------
687 *   *
688 *   * SPDX-License-Identifier: LGPL-2.1-or-later
689 *   * Copyright (C) 2015 by David Wells
690 *   *
691 *   * This file is part of the deal.II code gallery.
692 *   *
693 *   * -----------------------------------------------------------------------------
694 *   */
695 *  
696 *   #ifndef dealii__cdr_write_pvtu_output_templates_h
697 *   #define dealii__cdr_write_pvtu_output_templates_h
698 *   #include <deal.II/base/data_out_base.h>
699 *   #include <deal.II/base/utilities.h>
700 *  
701 *   #include <deal.II/dofs/dof_handler.h>
702 *  
703 *   #include <deal.II/lac/vector.h>
704 *  
705 *   #include <deal.II/numerics/data_out.h>
706 *  
707 *   #include <deal.II-cdr/write_pvtu_output.h>
708 *  
709 *   #include <fstream>
710 *   #include <string>
711 *   #include <vector>
712 *  
713 * @endcode
714 *
715 * Here is the implementation of the important function. This is similar to
716 * what is presented in @ref step_40 "step-40".
717 *
718 * @code
719 *   namespace CDR
720 *   {
721 *   using namespace dealii;
722 *  
723 *   template <int dim, typename VectorType>
724 *   void
725 *   WritePVTUOutput::write_output(const DoFHandler<dim> &dof_handler,
726 *   const VectorType & solution,
727 *   const unsigned int time_step_n,
728 *   const double current_time)
729 *   {
730 *   DataOut<dim> data_out;
731 *   data_out.attach_dof_handler(dof_handler);
732 *   data_out.add_data_vector(solution, "u");
733 *  
734 *   const auto & triangulation = dof_handler.get_triangulation();
735 *   Vector<float> subdomain(triangulation.n_active_cells());
736 *   for (auto &domain : subdomain)
737 *   {
738 *   domain = triangulation.locally_owned_subdomain();
739 *   }
740 *   data_out.add_data_vector(subdomain, "subdomain");
741 *   data_out.build_patches(patch_level);
742 *  
743 *   DataOutBase::VtkFlags flags;
744 *   flags.time = current_time;
745 * @endcode
746 *
747 * While the default flag is for the best compression level, using
748 * <code>best_speed</code> makes this function much faster.
749 *
750 * @code
751 *   flags.compression_level =
752 *   DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
753 *   data_out.set_flags(flags);
754 *  
755 *   unsigned int subdomain_n;
756 *   if (Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1)
757 *   {
758 *   subdomain_n = 0;
759 *   }
760 *   else
761 *   {
762 *   subdomain_n = triangulation.locally_owned_subdomain();
763 *   }
764 *  
765 *   std::ofstream output("solution-" + Utilities::int_to_string(time_step_n) +
766 *   "." + Utilities::int_to_string(subdomain_n, 4) +
767 *   ".vtu");
768 *  
769 *   data_out.write_vtu(output);
770 *  
771 *   if (this_mpi_process == 0)
772 *   {
773 *   std::vector<std::string> filenames;
774 *   for (unsigned int i = 0;
775 *   i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
776 *   ++i)
777 *   filenames.push_back("solution-" +
778 *   Utilities::int_to_string(time_step_n) + "." +
779 *   Utilities::int_to_string(i, 4) + ".vtu");
780 *   std::ofstream master_output(
781 *   "solution-" + Utilities::int_to_string(time_step_n) + ".pvtu");
782 *   data_out.write_pvtu_record(master_output, filenames);
783 *   }
784 *   }
785 *   } // namespace CDR
786 *   #endif
787 * @endcode
788
789
790<a name="ann-common/source/parameters.cc"></a>
791<h1>Annotated version of common/source/parameters.cc</h1>
792 *
793 *
794 *
795 * -----------------------------------------------------------------------------
796 *
797
798 *
799 * SPDX-License-Identifier: LGPL-2.1-or-later
800 * Copyright (C) 2015 by David Wells
801 *
802
803 *
804 * This file is part of the deal.II code gallery.
805 *
806
807 *
808 * -----------------------------------------------------------------------------
809 *
810
811 *
812 *
813 * @code
814 *   #include <deal.II-cdr/parameters.h>
815 *  
816 *   #include <fstream>
817 *   #include <string>
818 *  
819 *   namespace CDR
820 *   {
821 *   void
822 *   Parameters::configure_parameter_handler(ParameterHandler &parameter_handler)
823 *   {
824 *   parameter_handler.enter_subsection("Geometry");
825 *   {
826 *   parameter_handler.declare_entry("inner_radius",
827 *   "1.0",
828 *   Patterns::Double(0.0),
829 *   "Inner radius.");
830 *   parameter_handler.declare_entry("outer_radius",
831 *   "2.0",
832 *   Patterns::Double(0.0),
833 *   "Outer radius.");
834 *   }
835 *   parameter_handler.leave_subsection();
836 *  
837 *   parameter_handler.enter_subsection("Physical Parameters");
838 *   {
839 *   parameter_handler.declare_entry("diffusion_coefficient",
840 *   "1.0",
841 *   Patterns::Double(0.0),
842 *   "Diffusion coefficient.");
843 *   parameter_handler.declare_entry("reaction_coefficient",
844 *   "1.0",
845 *   Patterns::Double(0.0),
846 *   "Reaction coefficient.");
847 *   parameter_handler.declare_entry("time_dependent_forcing",
848 *   "true",
849 *   Patterns::Bool(),
850 *   "Whether or not "
851 *   "the forcing function depends on time.");
852 *   }
853 *   parameter_handler.leave_subsection();
854 *  
855 *   parameter_handler.enter_subsection("Finite Element");
856 *   {
857 *   parameter_handler.declare_entry("initial_refinement_level",
858 *   "1",
859 *   Patterns::Integer(1),
860 *   "Initial number of levels in the mesh.");
861 *   parameter_handler.declare_entry("max_refinement_level",
862 *   "1",
863 *   Patterns::Integer(1),
864 *   "Maximum number of levels in the mesh.");
865 *   parameter_handler.declare_entry("fe_order",
866 *   "1",
867 *   Patterns::Integer(1),
868 *   "Finite element order.");
869 *   }
870 *   parameter_handler.leave_subsection();
871 *  
872 *   parameter_handler.enter_subsection("Time Step");
873 *   {
874 *   parameter_handler.declare_entry("start_time",
875 *   "0.0",
876 *   Patterns::Double(0.0),
877 *   "Start time.");
878 *   parameter_handler.declare_entry("stop_time",
879 *   "1.0",
880 *   Patterns::Double(1.0),
881 *   "Stop time.");
882 *   parameter_handler.declare_entry("n_time_steps",
883 *   "1",
884 *   Patterns::Integer(1),
885 *   "Number of time steps.");
886 *   }
887 *   parameter_handler.leave_subsection();
888 *  
889 *   parameter_handler.enter_subsection("Output");
890 *   {
891 *   parameter_handler.declare_entry("save_interval",
892 *   "10",
893 *   Patterns::Integer(1),
894 *   "Save interval.");
895 *   parameter_handler.declare_entry("patch_level",
896 *   "2",
897 *   Patterns::Integer(0),
898 *   "Patch level.");
899 *   }
900 *   parameter_handler.leave_subsection();
901 *   }
902 *  
903 *   void
904 *   Parameters::read_parameter_file(const std::string &file_name)
905 *   {
906 *   ParameterHandler parameter_handler;
907 *   {
908 *   std::ifstream file(file_name);
909 *   configure_parameter_handler(parameter_handler);
910 *   parameter_handler.parse_input(file);
911 *   }
912 *  
913 *   parameter_handler.enter_subsection("Geometry");
914 *   {
915 *   inner_radius = parameter_handler.get_double("inner_radius");
916 *   outer_radius = parameter_handler.get_double("outer_radius");
917 *   }
918 *   parameter_handler.leave_subsection();
919 *  
920 *   parameter_handler.enter_subsection("Physical Parameters");
921 *   {
922 *   diffusion_coefficient =
923 *   parameter_handler.get_double("diffusion_coefficient");
924 *   reaction_coefficient =
925 *   parameter_handler.get_double("reaction_coefficient");
926 *   time_dependent_forcing =
927 *   parameter_handler.get_bool("time_dependent_forcing");
928 *   }
929 *   parameter_handler.leave_subsection();
930 *  
931 *  
932 *   parameter_handler.enter_subsection("Finite Element");
933 *   {
934 *   initial_refinement_level =
935 *   parameter_handler.get_integer("initial_refinement_level");
936 *   max_refinement_level =
937 *   parameter_handler.get_integer("max_refinement_level");
938 *   fe_order = parameter_handler.get_integer("fe_order");
939 *   }
940 *   parameter_handler.leave_subsection();
941 *  
942 *   parameter_handler.enter_subsection("Time Step");
943 *   {
944 *   start_time = parameter_handler.get_double("start_time");
945 *   stop_time = parameter_handler.get_double("stop_time");
946 *   n_time_steps = parameter_handler.get_integer("n_time_steps");
947 *   }
948 *   parameter_handler.leave_subsection();
949 *  
950 *   parameter_handler.enter_subsection("Output");
951 *   {
952 *   save_interval = parameter_handler.get_integer("save_interval");
953 *   patch_level = parameter_handler.get_integer("patch_level");
954 *   }
955 *   parameter_handler.leave_subsection();
956 *   }
957 *   } // namespace CDR
958 * @endcode
959
960
961<a name="ann-common/source/system_matrix.cc"></a>
962<h1>Annotated version of common/source/system_matrix.cc</h1>
963 *
964 *
965 *
966 *
967 * @code
968 *   /* -----------------------------------------------------------------------------
969 *   *
970 *   * SPDX-License-Identifier: LGPL-2.1-or-later
971 *   * Copyright (C) 2015 by David Wells
972 *   *
973 *   * This file is part of the deal.II code gallery.
974 *   *
975 *   * -----------------------------------------------------------------------------
976 *   */
977 *  
978 *   #include <deal.II/lac/sparse_matrix.h>
979 *   #include <deal.II/lac/trilinos_sparse_matrix.h>
980 *  
981 *   #include <deal.II-cdr/parameters.h>
982 *   #include <deal.II-cdr/system_matrix.h>
983 *   #include <deal.II-cdr/system_matrix.templates.h>
984 *  
985 * @endcode
986 *
987 * This file exists just to build template specializations of
988 * <code>create_system_matrix</code>. Even though the solver is run in
989 * parallel with Trilinos objects, other serial solvers can use the same
990 * function without recompilation by compiling everything here just one time.
991 *
992 * @code
993 *   namespace CDR
994 *   {
995 *   using namespace dealii;
996 *  
997 *   template void
998 *   create_system_matrix<2, SparseMatrix<double>>(
999 *   const DoFHandler<2> & dof_handler,
1000 *   const QGauss<2> & quad,
1001 *   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
1002 *   const CDR::Parameters & parameters,
1003 *   const double time_step,
1004 *   SparseMatrix<double> & system_matrix);
1005 *  
1006 *   template void
1007 *   create_system_matrix<3, SparseMatrix<double>>(
1008 *   const DoFHandler<3> & dof_handler,
1009 *   const QGauss<3> & quad,
1010 *   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
1011 *   const CDR::Parameters & parameters,
1012 *   const double time_step,
1013 *   SparseMatrix<double> & system_matrix);
1014 *  
1015 *   template void
1016 *   create_system_matrix<2, SparseMatrix<double>>(
1017 *   const DoFHandler<2> & dof_handler,
1018 *   const QGauss<2> & quad,
1019 *   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
1020 *   const CDR::Parameters & parameters,
1021 *   const double time_step,
1022 *   const AffineConstraints<double> & constraints,
1023 *   SparseMatrix<double> & system_matrix);
1024 *  
1025 *   template void
1026 *   create_system_matrix<3, SparseMatrix<double>>(
1027 *   const DoFHandler<3> & dof_handler,
1028 *   const QGauss<3> & quad,
1029 *   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
1030 *   const CDR::Parameters & parameters,
1031 *   const double time_step,
1032 *   const AffineConstraints<double> & constraints,
1033 *   SparseMatrix<double> & system_matrix);
1034 *  
1035 *   template void
1036 *   create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
1037 *   const DoFHandler<2> & dof_handler,
1038 *   const QGauss<2> & quad,
1039 *   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
1040 *   const CDR::Parameters & parameters,
1041 *   const double time_step,
1042 *   TrilinosWrappers::SparseMatrix & system_matrix);
1043 *  
1044 *   template void
1045 *   create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
1046 *   const DoFHandler<3> & dof_handler,
1047 *   const QGauss<3> & quad,
1048 *   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
1049 *   const CDR::Parameters & parameters,
1050 *   const double time_step,
1051 *   TrilinosWrappers::SparseMatrix & system_matrix);
1052 *  
1053 *   template void
1054 *   create_system_matrix<2, TrilinosWrappers::SparseMatrix>(
1055 *   const DoFHandler<2> & dof_handler,
1056 *   const QGauss<2> & quad,
1057 *   const std::function<Tensor<1, 2>(const Point<2>)> &convection_function,
1058 *   const CDR::Parameters & parameters,
1059 *   const double time_step,
1060 *   const AffineConstraints<double> & constraints,
1061 *   TrilinosWrappers::SparseMatrix & system_matrix);
1062 *  
1063 *   template void
1064 *   create_system_matrix<3, TrilinosWrappers::SparseMatrix>(
1065 *   const DoFHandler<3> & dof_handler,
1066 *   const QGauss<3> & quad,
1067 *   const std::function<Tensor<1, 3>(const Point<3>)> &convection_function,
1068 *   const CDR::Parameters & parameters,
1069 *   const double time_step,
1070 *   const AffineConstraints<double> & constraints,
1071 *   TrilinosWrappers::SparseMatrix & system_matrix);
1072 *   } // namespace CDR
1073 * @endcode
1074
1075
1076<a name="ann-common/source/system_rhs.cc"></a>
1077<h1>Annotated version of common/source/system_rhs.cc</h1>
1078 *
1079 *
1080 *
1081 *
1082 * @code
1083 *   /* -----------------------------------------------------------------------------
1084 *   *
1085 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1086 *   * Copyright (C) 2015 by David Wells
1087 *   *
1088 *   * This file is part of the deal.II code gallery.
1089 *   *
1090 *   * -----------------------------------------------------------------------------
1091 *   */
1092 *  
1093 *   #include <deal.II/lac/sparse_matrix.h>
1094 *   #include <deal.II/lac/trilinos_sparse_matrix.h>
1095 *   #include <deal.II/lac/trilinos_vector.h>
1096 *   #include <deal.II/lac/vector.h>
1097 *  
1098 *   #include <deal.II-cdr/system_rhs.templates.h>
1099 *  
1100 * @endcode
1101 *
1102 * Like <code>system_matrix.cc</code>, this file just compiles template
1103 * specializations.
1104 *
1105 * @code
1106 *   namespace CDR
1107 *   {
1108 *   using namespace dealii;
1109 *  
1110 *   template void
1111 *   create_system_rhs<2, Vector<double>>(
1112 *   const DoFHandler<2> & dof_handler,
1113 *   const QGauss<2> & quad,
1114 *   const std::function<Tensor<1, 2>(const Point<2>)> & convection_function,
1115 *   const std::function<double(double, const Point<2>)> &forcing_function,
1116 *   const CDR::Parameters & parameters,
1117 *   const Vector<double> & previous_solution,
1118 *   const AffineConstraints<double> & constraints,
1119 *   const double current_time,
1120 *   Vector<double> & system_rhs);
1121 *  
1122 *   template void
1123 *   create_system_rhs<3, Vector<double>>(
1124 *   const DoFHandler<3> & dof_handler,
1125 *   const QGauss<3> & quad,
1126 *   const std::function<Tensor<1, 3>(const Point<3>)> & convection_function,
1127 *   const std::function<double(double, const Point<3>)> &forcing_function,
1128 *   const CDR::Parameters & parameters,
1129 *   const Vector<double> & previous_solution,
1130 *   const AffineConstraints<double> & constraints,
1131 *   const double current_time,
1132 *   Vector<double> & system_rhs);
1133 *  
1134 *   template void
1135 *   create_system_rhs<2, TrilinosWrappers::MPI::Vector>(
1136 *   const DoFHandler<2> & dof_handler,
1137 *   const QGauss<2> & quad,
1138 *   const std::function<Tensor<1, 2>(const Point<2>)> & convection_function,
1139 *   const std::function<double(double, const Point<2>)> &forcing_function,
1140 *   const CDR::Parameters & parameters,
1141 *   const TrilinosWrappers::MPI::Vector & previous_solution,
1142 *   const AffineConstraints<double> & constraints,
1143 *   const double current_time,
1144 *   TrilinosWrappers::MPI::Vector & system_rhs);
1145 *  
1146 *   template void
1147 *   create_system_rhs<3, TrilinosWrappers::MPI::Vector>(
1148 *   const DoFHandler<3> & dof_handler,
1149 *   const QGauss<3> & quad,
1150 *   const std::function<Tensor<1, 3>(const Point<3>)> & convection_function,
1151 *   const std::function<double(double, const Point<3>)> &forcing_function,
1152 *   const CDR::Parameters & parameters,
1153 *   const TrilinosWrappers::MPI::Vector & previous_solution,
1154 *   const AffineConstraints<double> & constraints,
1155 *   const double current_time,
1156 *   TrilinosWrappers::MPI::Vector & system_rhs);
1157 *   } // namespace CDR
1158 * @endcode
1159
1160
1161<a name="ann-common/source/write_pvtu_output.cc"></a>
1162<h1>Annotated version of common/source/write_pvtu_output.cc</h1>
1163 *
1164 *
1165 *
1166 *
1167 * @code
1168 *   /* -----------------------------------------------------------------------------
1169 *   *
1170 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1171 *   * Copyright (C) 2015 by David Wells
1172 *   *
1173 *   * This file is part of the deal.II code gallery.
1174 *   *
1175 *   * -----------------------------------------------------------------------------
1176 *   */
1177 *  
1178 *   #include <deal.II/lac/trilinos_vector.h>
1179 *   #include <deal.II/lac/vector.h>
1180 *  
1181 *   #include <deal.II-cdr/write_pvtu_output.templates.h>
1182 *  
1183 * @endcode
1184 *
1185 * Again, this file just compiles the constructor and also the templated
1186 * functions.
1187 *
1188 * @code
1189 *   namespace CDR
1190 *   {
1191 *   using namespace dealii;
1192 *  
1193 *   WritePVTUOutput::WritePVTUOutput(const unsigned int patch_level)
1194 *   : patch_level{patch_level}
1195 *   , this_mpi_process{Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)}
1196 *   {}
1197 *  
1198 *   template void
1199 *   WritePVTUOutput::write_output(const DoFHandler<2> & dof_handler,
1200 *   const Vector<double> &solution,
1201 *   const unsigned int time_step_n,
1202 *   const double current_time);
1203 *  
1204 *   template void
1205 *   WritePVTUOutput::write_output(const DoFHandler<3> & dof_handler,
1206 *   const Vector<double> &solution,
1207 *   const unsigned int time_step_n,
1208 *   const double current_time);
1209 *  
1210 *   template void
1211 *   WritePVTUOutput::write_output(const DoFHandler<2> &dof_handler,
1212 *   const TrilinosWrappers::MPI::Vector &solution,
1213 *   const unsigned int time_step_n,
1214 *   const double current_time);
1215 *  
1216 *   template void
1217 *   WritePVTUOutput::write_output(const DoFHandler<3> &dof_handler,
1218 *   const TrilinosWrappers::MPI::Vector &solution,
1219 *   const unsigned int time_step_n,
1220 *   const double current_time);
1221 *   } // namespace CDR
1222 * @endcode
1223
1224
1225<a name="ann-solver/cdr.cc"></a>
1226<h1>Annotated version of solver/cdr.cc</h1>
1227 *
1228 *
1229 *
1230 *
1231 * @code
1232 *   /* -----------------------------------------------------------------------------
1233 *   *
1234 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1235 *   * Copyright (C) 2015 by David Wells
1236 *   *
1237 *   * This file is part of the deal.II code gallery.
1238 *   *
1239 *   * -----------------------------------------------------------------------------
1240 *   */
1241 *  
1242 *   #include <deal.II/base/conditional_ostream.h>
1243 *   #include <deal.II/base/quadrature_lib.h>
1244 *  
1245 *   #include <deal.II/dofs/dof_handler.h>
1246 *   #include <deal.II/dofs/dof_tools.h>
1247 *  
1248 *   #include <deal.II/fe/fe_q.h>
1249 *  
1250 *   #include <deal.II/grid/grid_generator.h>
1251 *   #include <deal.II/grid/manifold_lib.h>
1252 *  
1253 *   #include <deal.II/lac/affine_constraints.h>
1254 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
1255 *  
1256 *   #include <deal.II/numerics/error_estimator.h>
1257 *  
1258 * @endcode
1259 *
1260 * These headers are for distributed computations:
1261 *
1262 * @code
1263 *   #include <deal.II/base/index_set.h>
1264 *   #include <deal.II/base/utilities.h>
1265 *  
1266 *   #include <deal.II/distributed/grid_refinement.h>
1267 *   #include <deal.II/distributed/solution_transfer.h>
1268 *   #include <deal.II/distributed/tria.h>
1269 *  
1270 *   #include <deal.II/lac/sparsity_tools.h>
1271 *   #include <deal.II/lac/trilinos_precondition.h>
1272 *   #include <deal.II/lac/trilinos_solver.h>
1273 *   #include <deal.II/lac/trilinos_sparse_matrix.h>
1274 *   #include <deal.II/lac/trilinos_vector.h>
1275 *  
1276 *   #include <deal.II-cdr/parameters.h>
1277 *   #include <deal.II-cdr/system_matrix.h>
1278 *   #include <deal.II-cdr/system_rhs.h>
1279 *   #include <deal.II-cdr/write_pvtu_output.h>
1280 *  
1281 *   #include <chrono>
1282 *   #include <functional>
1283 *   #include <iostream>
1284 *  
1285 *   using namespace dealii;
1286 *  
1287 *   constexpr int manifold_id{0};
1288 *  
1289 * @endcode
1290 *
1291 * This is the actual solver class which performs time iteration and calls the
1292 * appropriate library functions to do it.
1293 *
1294 * @code
1295 *   template <int dim>
1296 *   class CDRProblem
1297 *   {
1298 *   public:
1299 *   CDRProblem(const CDR::Parameters &parameters);
1300 *   void
1301 *   run();
1302 *  
1303 *   private:
1304 *   const CDR::Parameters parameters;
1305 *   const double time_step;
1306 *   double current_time;
1307 *  
1308 *   MPI_Comm mpi_communicator;
1309 *   const unsigned int n_mpi_processes;
1310 *   const unsigned int this_mpi_process;
1311 *  
1312 *   FE_Q<dim> fe;
1313 *   QGauss<dim> quad;
1314 *   const SphericalManifold<dim> boundary_description;
1315 *   parallel::distributed::Triangulation<dim> triangulation;
1316 *   DoFHandler<dim> dof_handler;
1317 *  
1318 *   const std::function<Tensor<1, dim>(const Point<dim>)> convection_function;
1319 *   const std::function<double(double, const Point<dim>)> forcing_function;
1320 *  
1321 *   IndexSet locally_owned_dofs;
1322 *   IndexSet locally_relevant_dofs;
1323 *  
1324 *   AffineConstraints<double> constraints;
1325 *   bool first_run;
1326 *  
1327 * @endcode
1328 *
1329 * As is usual in parallel programs, I keep two copies of parts of the
1330 * complete solution: <code>locally_relevant_solution</code> contains both
1331 * the locally calculated solution as well as the layer of cells at its
1332 * boundary (the @ref GlossGhostCell "ghost cells") while
1333 * <code>completely_distributed_solution</code> only contains the parts of
1334 * the solution computed on the current @ref GlossMPIProcess "MPI process".
1335 *
1336 * @code
1337 *   TrilinosWrappers::MPI::Vector locally_relevant_solution;
1338 *   TrilinosWrappers::MPI::Vector completely_distributed_solution;
1339 *   TrilinosWrappers::MPI::Vector system_rhs;
1340 *   TrilinosWrappers::SparseMatrix system_matrix;
1341 *   TrilinosWrappers::PreconditionAMG preconditioner;
1342 *  
1343 *   ConditionalOStream pcout;
1344 *  
1345 *   void
1346 *   setup_geometry();
1347 *   void
1348 *   setup_system();
1349 *   void
1350 *   setup_dofs();
1351 *   void
1352 *   refine_mesh();
1353 *   void
1354 *   time_iterate();
1355 *   };
1356 *  
1357 *  
1358 *   template <int dim>
1359 *   CDRProblem<dim>::CDRProblem(const CDR::Parameters &parameters)
1360 *   : parameters(parameters)
1361 *   , time_step{(parameters.stop_time - parameters.start_time) /
1362 *   parameters.n_time_steps}
1363 *   , current_time{parameters.start_time}
1364 *   , mpi_communicator(MPI_COMM_WORLD)
1365 *   , n_mpi_processes{Utilities::MPI::n_mpi_processes(mpi_communicator)}
1366 *   , this_mpi_process{Utilities::MPI::this_mpi_process(mpi_communicator)}
1367 *   , fe(parameters.fe_order)
1368 *   , quad(parameters.fe_order + 2)
1369 *   , boundary_description(Point<dim>())
1370 *   , triangulation(mpi_communicator,
1371 *   typename Triangulation<dim>::MeshSmoothing(
1372 *   Triangulation<dim>::smoothing_on_refinement |
1373 *   Triangulation<dim>::smoothing_on_coarsening))
1374 *   , dof_handler(triangulation)
1375 *   , convection_function{[](const Point<dim> p) -> Tensor<1, dim> {
1376 *   Tensor<1, dim> v;
1377 *   v[0] = -p[1];
1378 *   v[1] = p[0];
1379 *   return v;
1380 *   }}
1381 *   , forcing_function{[](double t, const Point<dim> p) -> double {
1382 *   return std::exp(-8 * t) *
1383 *   std::exp(-40 * Utilities::fixed_power<6>(p[0] - 1.5)) *
1384 *   std::exp(-40 * Utilities::fixed_power<6>(p[1]));
1385 *   }}
1386 *   , first_run{true}
1387 *   , pcout(std::cout, this_mpi_process == 0)
1388 *   {
1389 *   Assert(dim == 2, ExcNotImplemented());
1390 *   }
1391 *  
1392 *  
1393 *   template <int dim>
1394 *   void
1395 *   CDRProblem<dim>::setup_geometry()
1396 *   {
1397 *   const Point<dim> center;
1398 *   GridGenerator::hyper_shell(triangulation,
1399 *   center,
1400 *   parameters.inner_radius,
1401 *   parameters.outer_radius);
1402 *   triangulation.set_manifold(manifold_id, boundary_description);
1403 *   for (const auto &cell : triangulation.active_cell_iterators())
1404 *   {
1405 *   cell->set_all_manifold_ids(manifold_id);
1406 *   }
1407 *   triangulation.refine_global(parameters.initial_refinement_level);
1408 *   }
1409 *  
1410 *  
1411 *   template <int dim>
1412 *   void
1413 *   CDRProblem<dim>::setup_dofs()
1414 *   {
1415 *   dof_handler.distribute_dofs(fe);
1416 *   pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
1417 *   << std::endl;
1418 *   locally_owned_dofs = dof_handler.locally_owned_dofs();
1419 *   locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
1420 *  
1421 *   constraints.clear();
1422 *   constraints.reinit(locally_relevant_dofs);
1423 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1424 *   DoFTools::make_zero_boundary_constraints(dof_handler,
1425 *   manifold_id,
1426 *   constraints);
1427 *   constraints.close();
1428 *  
1429 *   completely_distributed_solution.reinit(locally_owned_dofs, mpi_communicator);
1430 *  
1431 *   locally_relevant_solution.reinit(locally_owned_dofs,
1432 *   locally_relevant_dofs,
1433 *   mpi_communicator);
1434 *   }
1435 *  
1436 *  
1437 *   template <int dim>
1438 *   void
1439 *   CDRProblem<dim>::setup_system()
1440 *   {
1441 *   DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs());
1442 *   DoFTools::make_sparsity_pattern(dof_handler,
1443 *   dynamic_sparsity_pattern,
1444 *   constraints,
1445 *   /*keep_constrained_dofs*/ true);
1446 *   SparsityTools::distribute_sparsity_pattern(dynamic_sparsity_pattern,
1447 *   dof_handler.locally_owned_dofs(),
1448 *   mpi_communicator,
1449 *   locally_relevant_dofs);
1450 *  
1451 *   system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1452 *   system_matrix.reinit(locally_owned_dofs,
1453 *   dynamic_sparsity_pattern,
1454 *   mpi_communicator);
1455 *  
1456 *   CDR::create_system_matrix<dim>(dof_handler,
1457 *   quad,
1458 *   convection_function,
1459 *   parameters,
1460 *   time_step,
1461 *   constraints,
1462 *   system_matrix);
1463 *   system_matrix.compress(VectorOperation::add);
1464 *   preconditioner.initialize(system_matrix);
1465 *   }
1466 *  
1467 *  
1468 *   template <int dim>
1469 *   void
1470 *   CDRProblem<dim>::time_iterate()
1471 *   {
1472 *   double current_time = parameters.start_time;
1473 *   CDR::WritePVTUOutput pvtu_output(parameters.patch_level);
1474 *   for (unsigned int time_step_n = 0; time_step_n < parameters.n_time_steps;
1475 *   ++time_step_n)
1476 *   {
1477 *   current_time += time_step;
1478 *  
1479 *   system_rhs = 0.0;
1480 *   CDR::create_system_rhs<dim>(dof_handler,
1481 *   quad,
1482 *   convection_function,
1483 *   forcing_function,
1484 *   parameters,
1485 *   locally_relevant_solution,
1486 *   constraints,
1487 *   current_time,
1488 *   system_rhs);
1489 *   system_rhs.compress(VectorOperation::add);
1490 *  
1491 *   SolverControl solver_control(dof_handler.n_dofs(),
1492 *   1e-6 * system_rhs.l2_norm(),
1493 *   /*log_history = */ false,
1494 *   /*log_result = */ false);
1495 *   TrilinosWrappers::SolverGMRES solver(solver_control);
1496 *   solver.solve(system_matrix,
1497 *   completely_distributed_solution,
1498 *   system_rhs,
1499 *   preconditioner);
1500 *   constraints.distribute(completely_distributed_solution);
1501 *   locally_relevant_solution = completely_distributed_solution;
1502 *  
1503 *   if (time_step_n % parameters.save_interval == 0)
1504 *   {
1505 *   pvtu_output.write_output(dof_handler,
1506 *   locally_relevant_solution,
1507 *   time_step_n,
1508 *   current_time);
1509 *   }
1510 *  
1511 *   refine_mesh();
1512 *   }
1513 *   }
1514 *  
1515 *  
1516 *   template <int dim>
1517 *   void
1518 *   CDRProblem<dim>::refine_mesh()
1519 *   {
1520 *   using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;
1521 *  
1522 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1523 *   KellyErrorEstimator<dim>::estimate(dof_handler,
1524 *   QGauss<dim - 1>(fe.degree + 1),
1525 *   FunctionMap(),
1526 *   locally_relevant_solution,
1527 *   estimated_error_per_cell);
1528 *  
1529 * @endcode
1530 *
1531 * This solver uses a crude refinement strategy where cells with relatively
1532 * high errors are refined and cells with relatively low errors are
1533 * coarsened. The maximum refinement level is capped to prevent run-away
1534 * refinement.
1535 *
1536 * @code
1537 *   for (const auto &cell : triangulation.active_cell_iterators())
1538 *   {
1539 *   if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) >= 1e-3)
1540 *   {
1541 *   cell->set_refine_flag();
1542 *   }
1543 *   else if (std::abs(estimated_error_per_cell[cell->active_cell_index()]) <=
1544 *   1e-5)
1545 *   {
1546 *   cell->set_coarsen_flag();
1547 *   }
1548 *   }
1549 *  
1550 *   if (triangulation.n_levels() > parameters.max_refinement_level)
1551 *   {
1552 *   for (const auto &cell : triangulation.cell_iterators_on_level(
1553 *   parameters.max_refinement_level))
1554 *   {
1555 *   cell->clear_refine_flag();
1556 *   }
1557 *   }
1558 *  
1559 * @endcode
1560 *
1561 * Transferring the solution between different grids is ultimately just a
1562 * few function calls but they must be made in exactly the right order.
1563 *
1564 * @code
1565 *   parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
1566 *   solution_transfer(dof_handler);
1567 *  
1568 *   triangulation.prepare_coarsening_and_refinement();
1569 *   solution_transfer.prepare_for_coarsening_and_refinement(
1570 *   locally_relevant_solution);
1571 *   triangulation.execute_coarsening_and_refinement();
1572 *  
1573 *   setup_dofs();
1574 *  
1575 * @endcode
1576 *
1577 * The <code>solution_transfer</code> object stores a pointer to
1578 * <code>locally_relevant_solution</code>, so when
1579 * parallel::distributed::SolutionTransfer::interpolate is called it uses
1580 * those values to populate <code>temporary</code>.
1581 *
1582 * @code
1583 *   TrilinosWrappers::MPI::Vector temporary(locally_owned_dofs, mpi_communicator);
1584 *   solution_transfer.interpolate(temporary);
1585 * @endcode
1586 *
1587 * After <code>temporary</code> has the correct value, this call correctly
1588 * populates <code>completely_distributed_solution</code>, which had its
1589 * index set updated above with the call to <code>setup_dofs</code>.
1590 *
1591 * @code
1592 *   completely_distributed_solution = temporary;
1593 * @endcode
1594 *
1595 * Constraints cannot be applied to
1596 * @ref GlossGhostedVector "vectors with ghost entries" since the ghost
1597 * entries are write only, so this first goes through the completely
1598 * distributed vector.
1599 *
1600 * @code
1601 *   constraints.distribute(completely_distributed_solution);
1602 *   locally_relevant_solution = completely_distributed_solution;
1603 *   setup_system();
1604 *   }
1605 *  
1606 *  
1607 *   template <int dim>
1608 *   void
1609 *   CDRProblem<dim>::run()
1610 *   {
1611 *   setup_geometry();
1612 *   setup_dofs();
1613 *   setup_system();
1614 *   time_iterate();
1615 *   }
1616 *  
1617 *  
1618 *   constexpr int dim{2};
1619 *  
1620 *  
1621 *   int
1622 *   main(int argc, char *argv[])
1623 *   {
1624 * @endcode
1625 *
1626 * One of the new features in C++11 is the <code>chrono</code> component of
1627 * the standard library. This gives us an easy way to time the output.
1628 *
1629 * @code
1630 *   auto t0 = std::chrono::high_resolution_clock::now();
1631 *  
1632 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
1633 *   CDR::Parameters parameters;
1634 *   parameters.read_parameter_file("parameters.prm");
1635 *   CDRProblem<dim> cdr_problem(parameters);
1636 *   cdr_problem.run();
1637 *  
1638 *   auto t1 = std::chrono::high_resolution_clock::now();
1639 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
1640 *   {
1641 *   std::cout << "time elapsed: "
1642 *   << std::chrono::duration_cast<std::chrono::milliseconds>(t1 -
1643 *   t0)
1644 *   .count()
1645 *   << " milliseconds." << std::endl;
1646 *   }
1647 *  
1648 *   return 0;
1649 *   }
1650 * @endcode
1651
1652
1653*/
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
FE_Q
Definition fe_q.h:554
KellyErrorEstimator::estimate
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
MPI_Comm
Point
Definition point.h:111
Triangulation::n_active_cells
unsigned int n_active_cells() const
Triangulation::refine_global
void refine_global(const unsigned int times=1)
Triangulation::n_levels
unsigned int n_levels() const
Triangulation::get_triangulation
Triangulation< dim, spacedim > & get_triangulation()
Utilities::MPI::MPI_InitFinalize
Definition mpi.h:1081
parallel::TriangulationBase::locally_owned_subdomain
types::subdomain_id locally_owned_subdomain() const override
Definition tria_base.cc:345
parallel::distributed::SolutionTransfer::interpolate
void interpolate(std::vector< VectorType * > &all_out)
Definition solution_transfer.cc:232
parallel::distributed::Triangulation::execute_coarsening_and_refinement
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
parallel::distributed::Triangulation::prepare_coarsening_and_refinement
virtual bool prepare_coarsening_and_refinement() override
Definition tria.cc:2805
center
Point< 3 > center
Definition data_out_base.cc:267
first
Point< 2 > first
Definition grid_out.cc:4623
level
unsigned int level
Definition grid_out.cc:4626
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
Assert
#define Assert(cond, exc)
Definition exceptions.h:1638
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition dof_tools_constraints.cc:3073
DoFTools::make_sparsity_pattern
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)
Definition dof_tools_sparsity.cc:56
DoFTools::make_zero_boundary_constraints
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
Definition dof_tools_constraints.cc:4712
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:75
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:134
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:81
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition fe_update_flags.h:127
Triangulation::set_manifold
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
AdaptationStrategies::Refinement::split
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
DoFTools::extract_locally_relevant_dofs
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_tools.cc:1164
GridGenerator::hyper_shell
void hyper_shell(Triangulation< dim, spacedim > &tria, const Point< spacedim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:57
LocalIntegrators::Advection::cell_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
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SparsityTools::distribute_sparsity_pattern
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
Definition sparsity_tools.cc:1020
Utilities::MPI::ConsensusAlgorithms::serial
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
WorkStream::internal::tbb_no_coloring::run
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)
Definition work_stream.h:471
internal::p4est::functions
int(& functions)(const void *v1, const void *v2)
Definition p4est_wrappers.cc:338
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition matrix_block.h:617
std
STL namespace.
std::exp
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6829
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
types::manifold_id
unsigned int manifold_id
Definition types.h:156
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition p4est_wrappers.cc:68
Generated by doxygen 1.11.0