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