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
Maxwell-Eigenvalue-hp-Refinement.h
Go to the documentation of this file.
1
223 *  
224 *  
225 *   #include <deal.II/base/function_parser.h>
226 *   #include <deal.II/base/index_set.h>
227 *   #include <deal.II/base/parameter_handler.h>
228 *   #include <deal.II/base/quadrature_lib.h>
229 *   #include <deal.II/base/utilities.h>
230 *  
231 *   #include <deal.II/dofs/dof_handler.h>
232 *   #include <deal.II/dofs/dof_tools.h>
233 *  
234 *   #include <deal.II/fe/fe_nedelec.h>
235 *   #include <deal.II/fe/fe_series.h>
236 *   #include <deal.II/fe/fe_values.h>
237 *  
238 *   #include <deal.II/grid/tria.h>
239 *   #include <deal.II/grid/tria_iterator.h>
240 *  
241 *   #include <deal.II/lac/affine_constraints.h>
242 *   #include <deal.II/lac/full_matrix.h>
243 *   #include <deal.II/lac/petsc_precondition.h>
244 *   #include <deal.II/lac/petsc_sparse_matrix.h>
245 *   #include <deal.II/lac/petsc_vector.h>
246 *   #include <deal.II/lac/slepc_solver.h>
247 *  
248 *   #include <deal.II/numerics/data_out.h>
249 *   #include <deal.II/numerics/vector_tools.h>
250 *  
251 * @endcode
252 *
253 * For parallelization (using WorkStream and Intel TBB)
254 *
255 * @code
256 *   #include <deal.II/base/multithread_info.h>
257 *   #include <deal.II/base/work_stream.h>
258 *  
259 *   #include "petscpc.h"
260 *  
261 * @endcode
262 *
263 * For Error Estimation/Indication and Smoothness Indication
264 *
265 * @code
266 *   #include <deal.II/fe/fe_tools.h>
267 *  
268 *   #include <deal.II/numerics/error_estimator.h>
269 *   #include <deal.II/numerics/smoothness_estimator.h>
270 * @endcode
271 *
272 * For refinement
273 *
274 * @code
275 *   #include <deal.II/grid/grid_refinement.h>
276 *  
277 *   #include <fstream>
278 *   #include <iostream>
279 *   #include <memory>
280 *  
281 *   namespace Operations
282 *   {
283 *  
287 *   double
288 *   curlcurl(const ::FEValues<2> &fe_values,
289 *   const unsigned int & i,
290 *   const unsigned int & j,
291 *   const unsigned int & q_point)
292 *   {
293 *   auto gradu1_x1x2 = fe_values.shape_grad_component(i, q_point, 0);
294 *   auto gradu2_x1x2 = fe_values.shape_grad_component(i, q_point, 1);
295 *  
296 *   auto gradv1_x1x2 = fe_values.shape_grad_component(j, q_point, 0);
297 *   auto gradv2_x1x2 = fe_values.shape_grad_component(j, q_point, 1);
298 *   return (gradu2_x1x2[0] - gradu1_x1x2[1]) *
299 *   (gradv2_x1x2[0] - gradv1_x1x2[1]);
300 *   }
301 *  
302 *  
305 *   template <int dim>
306 *   inline double
307 *   dot_term(const ::FEValues<dim> &fe_values,
308 *   const unsigned int & i,
309 *   const unsigned int & j,
310 *   const unsigned int & q_point)
311 *   {
312 *   double output = 0.0;
313 *   for (unsigned int comp = 0; comp < dim; ++comp)
314 *   {
315 *   output += fe_values.shape_value_component(i, q_point, comp) *
316 *   fe_values.shape_value_component(j, q_point, comp);
317 *   }
318 *   return output;
319 *   }
320 *   } // namespace Operations
321 *  
326 *   namespace Structures
327 *   {
328 *   using namespace dealii;
329 *  
330 *   void
331 *   create_L_waveguide(Triangulation<2> &triangulation, const double &scaling)
332 *   {
333 *   const unsigned int dim = 2;
334 *  
335 *   const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
336 *   {scaling * 0.5, scaling * 0.0},
337 *   {scaling * 0.0, scaling * 0.5},
338 *   {scaling * 0.5, scaling * 0.5},
339 *   {scaling * 0.0, scaling * 1.0},
340 *   {scaling * 0.5, scaling * 1.0},
341 *   {scaling * 1.0, scaling * 0.5},
342 *   {scaling * 1.0, scaling * 1.0}};
343 *  
344 *   const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
345 *   cell_vertices = {{{0, 1, 2, 3}}, {{2, 3, 4, 5}}, {{3, 6, 5, 7}}};
346 *   const unsigned int n_cells = cell_vertices.size();
347 *   std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
348 *   for (unsigned int i = 0; i < n_cells; ++i)
349 *   {
350 *   for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
351 *   cells[i].vertices[j] = cell_vertices[i][j];
352 *   cells[i].material_id = 0;
353 *   }
354 *   triangulation.create_triangulation(vertices, cells, SubCellData());
355 *   triangulation.refine_global(1);
356 *   }
357 *  
358 *  
359 *   void
360 *   create_standard_waveguide(Triangulation<2> &triangulation,
361 *   const double & scaling)
362 *   {
363 *   const unsigned int dim = 2;
364 *  
365 *   const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
366 *   {scaling * 0.6, scaling * 0.0},
367 *   {scaling * 0.0, scaling * 0.3},
368 *   {scaling * 0.6, scaling * 0.3}};
369 *  
370 *   const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
371 *   cell_vertices = {{{0, 1, 2, 3}}};
372 *   const unsigned int n_cells = cell_vertices.size();
373 *   std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
374 *   for (unsigned int i = 0; i < n_cells; ++i)
375 *   {
376 *   for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
377 *   cells[i].vertices[j] = cell_vertices[i][j];
378 *   cells[i].material_id = 0;
379 *   }
380 *   triangulation.create_triangulation(vertices, cells, SubCellData());
381 *   triangulation.refine_global(0);
382 *   }
383 *   } // namespace Structures
384 *  
388 *   namespace Maxwell
389 *   {
390 *   using namespace dealii;
391 *  
392 *   /*
393 *   The "Base" class provides the universal functionality of any eigensolver,
394 *   namely the parameters for the problem, an underlying triangulation, and
395 *   functionality for setting the refinement cycle and to output the solution.
396 *  
397 *   In this case, and for any future class, the use of raw pointers (as opposed to
398 *   "smart" pointers) indicates a lack of ownership. Specifically, the
399 *   triangulation raw pointer is pointing to a triangulation that is owned (and
400 *   created) elsewhere.
401 *   */
402 *   template <int dim>
403 *   class Base
404 *   {
405 *   public:
406 *   Base(const std::string &prm_file, Triangulation<dim> &coarse_grid);
407 *  
408 *   virtual unsigned int
409 *   solve_problem() = 0; // Implemented by a derived class
410 *   virtual void
411 *   set_refinement_cycle(const unsigned int cycle);
412 *  
413 *   virtual void
414 *   output_solution() = 0; // Implemented by a derived class
415 *  
416 *  
417 *   protected:
418 *   Triangulation<dim> * triangulation;
419 *   unsigned int refinement_cycle = 0;
420 *   std::unique_ptr<ParameterHandler> parameters;
421 *   unsigned int n_eigenpairs = 1;
422 *   double target = 0.0;
423 *   unsigned int eigenpair_selection_scheme;
424 *   unsigned int max_cycles = 0;
425 *   ompi_communicator_t * mpi_communicator = PETSC_COMM_SELF;
426 *   };
427 *  
428 *  
431 *   template <int dim>
432 *   Base<dim>::Base(const std::string &prm_file, Triangulation<dim> &coarse_grid)
433 *   : triangulation(&coarse_grid)
434 *   , parameters(std::make_unique<ParameterHandler>())
435 *   {
436 *   parameters->declare_entry(
437 *   "Eigenpair selection scheme",
438 *   "1",
439 *   Patterns::Integer(0, 1),
440 *   "The type of eigenpairs to find (0 - smallest, 1 - target)");
441 *   parameters->declare_entry("Number of eigenvalues/eigenfunctions",
442 *   "1",
443 *   Patterns::Integer(0, 100),
444 *   "The number of eigenvalues/eigenfunctions "
445 *   "to be computed.");
446 *   parameters->declare_entry("Target eigenvalue",
447 *   "1",
448 *   Patterns::Anything(),
449 *   "The target eigenvalue (if scheme == 1)");
450 *  
451 *   parameters->declare_entry("Cycles number",
452 *   "1",
453 *   Patterns::Integer(0, 1500),
454 *   "The number of cycles in refinement");
455 *   parameters->parse_input(prm_file);
456 *  
457 *   eigenpair_selection_scheme =
458 *   parameters->get_integer("Eigenpair selection scheme");
459 *  
460 * @endcode
461 *
462 * The project currently only supports selection by a target eigenvalue.
463 * Furthermore, only one eigenpair can be computed at a time.
464 *
465 * @code
466 *   assert(eigenpair_selection_scheme == 1 &&
467 *   "Selection by a target is the only currently supported option!");
468 *   n_eigenpairs =
469 *   parameters->get_integer("Number of eigenvalues/eigenfunctions");
470 *   assert(
471 *   n_eigenpairs == 1 &&
472 *   "Only the computation of a single eigenpair is currently supported!");
473 *  
474 *   target = parameters->get_double("Target eigenvalue");
475 *   max_cycles = parameters->get_integer("Cycles number");
476 *   if (eigenpair_selection_scheme == 1)
477 *   n_eigenpairs = 1;
478 *   }
479 *  
480 *   template <int dim>
481 *   void
482 *   Base<dim>::set_refinement_cycle(const unsigned int cycle)
483 *   {
484 *   refinement_cycle = cycle;
485 *   }
486 *  
487 *  
492 *   template <int dim>
493 *   class EigenSolver : public virtual Base<dim>
494 *   {
495 *   public:
496 *   EigenSolver(const std::string & prm_file,
497 *   Triangulation<dim> &coarse_grid,
498 *   const unsigned int &minimum_degree,
499 *   const unsigned int &maximum_degree,
500 *   const unsigned int &starting_degree);
501 *  
502 *   virtual unsigned int
503 *   solve_problem() override;
504 *  
505 *   virtual unsigned int
506 *   n_dofs() const;
507 *  
508 *   template <class SolverType>
509 *   void
510 *   initialize_eigensolver(SolverType &eigensolver);
511 *  
512 *   virtual void
513 *   setup_system();
514 *  
515 *   virtual void
516 *   assemble_system();
517 *  
518 *   protected:
519 *   const std::unique_ptr<hp::FECollection<dim>> fe_collection;
520 *   std::unique_ptr<hp::QCollection<dim>> quadrature_collection;
521 *   std::unique_ptr<hp::QCollection<dim - 1>> face_quadrature_collection;
522 *   DoFHandler<dim> dof_handler;
523 *   const unsigned int max_degree, min_degree;
524 * @endcode
525 *
526 * for the actual solution
527 *
528 * @code
529 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> eigenfunctions;
530 *   std::unique_ptr<std::vector<double>> eigenvalues;
531 *   Vector<double> solution;
532 *  
533 *   double *
534 *   get_lambda_h();
535 *  
536 *   Vector<double> *
537 *   get_solution();
538 *  
539 *   void
540 *   convert_solution();
541 *  
542 *   private:
543 *   AffineConstraints<double> constraints;
544 *   PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
545 *   };
546 *  
547 *  
551 *   template <int dim>
552 *   EigenSolver<dim>::EigenSolver(const std::string & prm_file,
553 *   Triangulation<dim> &triangulation,
554 *   const unsigned int &minimum_degree,
555 *   const unsigned int &maximum_degree,
556 *   const unsigned int &starting_degree)
557 *   : Base<dim>(prm_file, triangulation)
558 *   , fe_collection(std::make_unique<hp::FECollection<dim>>())
559 *   , quadrature_collection(std::make_unique<hp::QCollection<dim>>())
560 *   , face_quadrature_collection(std::make_unique<hp::QCollection<dim - 1>>())
561 *   , dof_handler(triangulation)
562 *   , max_degree(maximum_degree)
563 *   , min_degree(minimum_degree)
564 *   , eigenfunctions(
565 *   std::make_unique<std::vector<PETScWrappers::MPI::Vector>>())
566 *   , eigenvalues(std::make_unique<std::vector<double>>())
567 *   {
568 *   for (unsigned int degree = min_degree; degree <= max_degree; ++degree)
569 *   {
570 *   fe_collection->push_back(FE_Nedelec<dim>(degree - 1));
571 * @endcode
572 *
573 * Generate quadrature collection with sorted quadrature weights
574 *
575 * @code
576 *   const QGauss<dim> quadrature(degree + 1);
577 *   const QSorted<dim> sorted_quadrature(quadrature);
578 *   quadrature_collection->push_back(sorted_quadrature);
579 *  
580 *   const QGauss<dim - 1> face_quadrature(degree + 1);
581 *   const QSorted<dim - 1> sorted_face_quadrature(face_quadrature);
582 *   face_quadrature_collection->push_back(sorted_face_quadrature);
583 *   }
584 * @endcode
585 *
586 * adjust the discretization
587 *
588 * @code
589 *   if (starting_degree > min_degree && starting_degree <= max_degree)
590 *   {
591 *   const unsigned int start_diff = starting_degree - min_degree;
592 *   typename DoFHandler<dim>::active_cell_iterator
593 *   cell1 = dof_handler.begin_active(),
594 *   endc1 = dof_handler.end();
595 *   for (; cell1 < endc1; ++cell1)
596 *   {
597 *   cell1->set_active_fe_index(start_diff);
598 *   }
599 *   }
600 *   }
601 *  
602 *  
606 *   template <int dim>
607 *   double *
608 *   EigenSolver<dim>::get_lambda_h()
609 *   {
610 *   return &(*eigenvalues)[0];
611 *   }
612 *  
613 *  
617 *   template <int dim>
618 *   Vector<double> *
619 *   EigenSolver<dim>::get_solution()
620 *   {
621 *   return &solution;
622 *   }
623 *  
624 *  
627 *   template <int dim>
628 *   void
629 *   EigenSolver<dim>::convert_solution()
630 *   {
631 *   solution.reinit((*eigenfunctions)[0].size());
632 *   for (unsigned int i = 0; i < solution.size(); ++i)
633 *   solution[i] = (*eigenfunctions)[0][i];
634 *   }
635 *  
636 *  
642 *   template <int dim>
643 *   template <class SolverType>
644 *   void
645 *   EigenSolver<dim>::initialize_eigensolver(SolverType &eigensolver)
646 *   {
647 * @endcode
648 *
649 * From the parameters class, initialize the eigensolver...
650 *
651 * @code
652 *   switch (this->eigenpair_selection_scheme)
653 *   {
654 *   case 1:
655 *   eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
656 * @endcode
657 *
658 * eigensolver.set_target_eigenvalue(this->target);
659 *
660 * @code
661 *   break;
662 *   default:
663 *   eigensolver.set_which_eigenpairs(EPS_SMALLEST_MAGNITUDE);
664 *  
665 *   break;
666 *   }
667 *   eigensolver.set_problem_type(EPS_GHEP);
668 * @endcode
669 *
670 * apply a Shift-Invert spectrum transformation
671 *
672
673 *
674 *
675 * @code
676 *   double shift_scalar = this->parameters->get_double("Target eigenvalue");
677 * @endcode
678 *
679 * //For the shift-and-invert transformation
680 *
681 * @code
682 *   SLEPcWrappers::TransformationShiftInvert::AdditionalData additional_data(
683 *   shift_scalar);
684 *   SLEPcWrappers::TransformationShiftInvert spectral_transformation(
685 *   this->mpi_communicator, additional_data);
686 *  
687 *   eigensolver.set_transformation(spectral_transformation);
688 *   eigensolver.set_target_eigenvalue(this->target);
689 *   }
690 *  
691 *  
695 *   template <int dim>
696 *   unsigned int
697 *   EigenSolver<dim>::solve_problem()
698 *   {
699 *   setup_system();
700 *   assemble_system();
701 *  
702 *   SolverControl solver_control(dof_handler.n_dofs() * 10,
703 *   5.0e-8,
704 *   false,
705 *   false);
706 *   SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control,
707 *   this->mpi_communicator);
708 *  
709 *   initialize_eigensolver(eigensolver);
710 *  
711 * @endcode
712 *
713 * solve the problem
714 *
715 * @code
716 *   eigensolver.solve(stiffness_matrix,
717 *   mass_matrix,
718 *   *eigenvalues,
719 *   *eigenfunctions,
720 *   eigenfunctions->size());
721 *   for (auto &entry : *eigenfunctions)
722 *   {
723 *   constraints.distribute(entry);
724 *   }
725 *   convert_solution();
726 *  
727 *   return solver_control.last_step();
728 *   }
729 *  
730 *   template <int dim>
731 *   unsigned int
732 *   EigenSolver<dim>::n_dofs() const
733 *   {
734 *   return dof_handler.n_dofs();
735 *   }
736 *  
737 *  
742 *   template <int dim>
743 *   void
744 *   EigenSolver<dim>::setup_system()
745 *   {
746 *   dof_handler.distribute_dofs(*fe_collection);
747 *   constraints.clear();
748 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
749 *   DoFTools::make_zero_boundary_constraints(dof_handler, constraints);
750 *   constraints.close();
751 *  
752 *   eigenfunctions->resize(this->n_eigenpairs);
753 *   eigenvalues->resize(this->n_eigenpairs);
754 *  
755 *   IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs();
756 *  
757 *   for (auto &entry : *eigenfunctions)
758 *   {
759 *   entry.reinit(eigenfunction_index_set, MPI_COMM_WORLD);
760 *   }
761 *   }
762 *  
763 *  
766 *   template <int dim>
767 *   void
768 *   EigenSolver<dim>::assemble_system()
769 *   {
770 *   hp::FEValues<dim> hp_fe_values(*fe_collection,
771 *   *quadrature_collection,
772 *   update_values | update_gradients |
773 *   update_quadrature_points |
774 *   update_JxW_values);
775 * @endcode
776 *
777 * Prep the system matrices for the solution
778 *
779 * @code
780 *   stiffness_matrix.reinit(dof_handler.n_dofs(),
781 *   dof_handler.n_dofs(),
782 *   dof_handler.max_couplings_between_dofs());
783 *   mass_matrix.reinit(dof_handler.n_dofs(),
784 *   dof_handler.n_dofs(),
785 *   dof_handler.max_couplings_between_dofs());
786 *  
787 *   FullMatrix<double> cell_stiffness_matrix, cell_mass_matrix;
788 *   std::vector<types::global_dof_index> local_dof_indices;
789 *  
790 *   for (const auto &cell : dof_handler.active_cell_iterators())
791 *   {
792 *   const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
793 *  
794 *   cell_stiffness_matrix.reinit(dofs_per_cell, dofs_per_cell);
795 *   cell_stiffness_matrix = 0;
796 *  
797 *   cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
798 *   cell_mass_matrix = 0;
799 *  
800 *   hp_fe_values.reinit(cell);
801 *  
802 *   const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
803 *  
804 *   for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
805 *   ++q_point)
806 *   {
807 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
808 *   {
809 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
810 *   {
811 * @endcode
812 *
813 * Note that (in general) the Nedelec element is not
814 * primitive, namely that the shape functions are vectorial
815 * with components in more than one direction
816 *
817
818 *
819 *
820 * @code
821 *   cell_stiffness_matrix(i, j) +=
822 *   Operations::curlcurl(fe_values, i, j, q_point) *
823 *   fe_values.JxW(q_point);
824 *  
825 *   cell_mass_matrix(i, j) +=
826 *   (Operations::dot_term(fe_values, i, j, q_point)) *
827 *   fe_values.JxW(q_point);
828 *   }
829 *   }
830 *   local_dof_indices.resize(dofs_per_cell);
831 *   cell->get_dof_indices(local_dof_indices);
832 *   }
833 *  
834 *   constraints.distribute_local_to_global(cell_stiffness_matrix,
835 *   local_dof_indices,
836 *   stiffness_matrix);
837 *   constraints.distribute_local_to_global(cell_mass_matrix,
838 *   local_dof_indices,
839 *   mass_matrix);
840 *   }
841 *   stiffness_matrix.compress(VectorOperation::add);
842 *   mass_matrix.compress(VectorOperation::add);
843 *  
844 *   for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
845 *   if (constraints.is_constrained(i))
846 *   {
847 *   stiffness_matrix.set(i, i, 10000.0);
848 *   mass_matrix.set(i, i, 1);
849 *   }
850 * @endcode
851 *
852 * since we have just set individual elements, we need the following
853 *
854 * @code
855 *   stiffness_matrix.compress(VectorOperation::insert);
856 *   mass_matrix.compress(VectorOperation::insert);
857 *   }
858 *  
859 *  
863 *   template <int dim>
864 *   class PrimalSolver : public EigenSolver<dim>
865 *   {
866 *   public:
867 *   PrimalSolver(const std::string & prm_file,
868 *   Triangulation<dim> &triangulation,
869 *   const unsigned int &min_degree,
870 *   const unsigned int &max_degree,
871 *   const unsigned int &starting_degree);
872 *  
873 *   virtual void
874 *   output_solution()
875 *   override; // Implements the output solution of the base class...
876 *   virtual unsigned int
877 *   n_dofs() const override;
878 *   };
879 *  
880 *   template <int dim>
881 *   PrimalSolver<dim>::PrimalSolver(const std::string & prm_file,
882 *   Triangulation<dim> &triangulation,
883 *   const unsigned int &min_degree,
884 *   const unsigned int &max_degree,
885 *   const unsigned int &starting_degree)
886 *   : Base<dim>(prm_file, triangulation)
887 *   , EigenSolver<dim>(prm_file,
888 *   triangulation,
889 *   min_degree,
890 *   max_degree,
891 *   starting_degree)
892 *   {}
893 *  
894 *  
898 *   template <int dim>
899 *   void
900 *   PrimalSolver<dim>::output_solution()
901 *   {
902 *   DataOut<dim> data_out;
903 *   data_out.attach_dof_handler(this->dof_handler);
904 *   Vector<double> fe_degrees(this->triangulation->n_active_cells());
905 *   for (const auto &cell : this->dof_handler.active_cell_iterators())
906 *   fe_degrees(cell->active_cell_index()) =
907 *   (*this->fe_collection)[cell->active_fe_index()].degree;
908 *   data_out.add_data_vector(fe_degrees, "fe_degree");
909 *   data_out.add_data_vector((*this->eigenfunctions)[0],
910 *   std::string("eigenfunction_no_") +
911 *   Utilities::int_to_string(0));
912 *  
913 *   std::cout << "Eigenvalue: " << (*this->eigenvalues)[0]
914 *   << " NDoFs: " << this->dof_handler.n_dofs() << std::endl;
915 *   std::ofstream eigenvalues_out(
916 *   "eigenvalues-" + std::to_string(this->refinement_cycle) + ".txt");
917 *  
918 *   eigenvalues_out << std::setprecision(20) << (*this->eigenvalues)[0] << " "
919 *   << this->dof_handler.n_dofs() << std::endl;
920 *  
921 *   eigenvalues_out.close();
922 *  
923 *  
924 *   data_out.build_patches();
925 *   std::ofstream output("eigenvectors-" +
926 *   std::to_string(this->refinement_cycle) + ".vtu");
927 *   data_out.write_vtu(output);
928 *   }
929 *  
930 *   template <int dim>
931 *   unsigned int
932 *   PrimalSolver<dim>::n_dofs() const
933 *   {
934 *   return EigenSolver<dim>::n_dofs();
935 *   }
936 *  
937 * @endcode
938 *
939 * Note, that at least for the demonstrated problem (i.e., a Hermitian problem
940 * and eigenvalue QoI), the dual problem is identical to the primal problem;
941 * however, it is convenient to separate them in this manner (e.g., for
942 * considering functionals of the eigenfunction).
943 *
944 * @code
945 *   template <int dim>
946 *   class DualSolver : public EigenSolver<dim>
947 *   {
948 *   public:
949 *   DualSolver(const std::string & prm_file,
950 *   Triangulation<dim> &triangulation,
951 *   const unsigned int &min_degree,
952 *   const unsigned int &max_degree,
953 *   const unsigned int &starting_degree);
954 *   };
955 *  
956 *   template <int dim>
957 *   DualSolver<dim>::DualSolver(const std::string & prm_file,
958 *   Triangulation<dim> &triangulation,
959 *   const unsigned int &min_degree,
960 *   const unsigned int &max_degree,
961 *   const unsigned int &starting_degree)
962 *   : Base<dim>(prm_file, triangulation)
963 *   , EigenSolver<dim>(prm_file,
964 *   triangulation,
965 *   min_degree,
966 *   max_degree,
967 *   starting_degree)
968 *   {}
969 *  
970 *   } // namespace Maxwell
971 *  
975 *   namespace ErrorIndicators
976 *   {
977 *   using namespace Maxwell;
978 *  
979 *  
984 *   template <int dim, bool report_dual>
985 *   class DualWeightedResidual : public PrimalSolver<dim>, public DualSolver<dim>
986 *   {
987 *   public:
988 *   void
989 *   output_eigenvalue_data(std::ofstream &os);
990 *   void
991 *   output_qoi_error_estimates(std::ofstream &os);
992 *  
993 *   std::string
994 *   name() const
995 *   {
996 *   return "DWR";
997 *   }
998 *   DualWeightedResidual(const std::string & prm_file,
999 *   Triangulation<dim> &triangulation,
1000 *   const unsigned int &min_primal_degree,
1001 *   const unsigned int &max_primal_degree,
1002 *   const unsigned int &starting_primal_degree);
1003 *  
1004 *   virtual unsigned int
1005 *   solve_problem() override;
1006 *  
1007 *   virtual void
1008 *   output_solution() override;
1009 *  
1010 *   virtual unsigned int
1011 *   n_dofs() const override;
1012 *  
1013 *   void
1014 *   estimate_error(Vector<double> &error_indicators);
1015 *  
1016 *   DoFHandler<dim> *
1017 *   get_DoFHandler();
1018 *  
1019 *   DoFHandler<dim> *
1020 *   get_primal_DoFHandler();
1021 *  
1022 *   DoFHandler<dim> *
1023 *   get_dual_DoFHandler();
1024 *  
1025 *   hp::FECollection<dim> *
1026 *   get_FECollection();
1027 *  
1028 *   hp::FECollection<dim> *
1029 *   get_primal_FECollection();
1030 *  
1031 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1032 *   get_eigenfunctions();
1033 *  
1034 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1035 *   get_primal_eigenfunctions();
1036 *  
1037 *   std::unique_ptr<std::vector<double>> &
1038 *   get_primal_eigenvalues();
1039 *  
1040 *   std::unique_ptr<std::vector<double>> &
1041 *   get_dual_eigenvalues();
1042 *  
1043 *   void
1044 *   synchronize_discretization();
1045 *  
1046 *   unsigned int
1047 *   get_max_degree()
1048 *   {
1049 *   return PrimalSolver<dim>::fe_collection->max_degree();
1050 *   }
1051 *   double qoi_error_estimate = 0;
1052 *  
1053 *   private:
1054 *   void
1055 *   embed(const DoFHandler<dim> & dof1,
1056 *   const DoFHandler<dim> & dof2,
1057 *   const AffineConstraints<double> &constraints,
1058 *   const Vector<double> & solution,
1059 *   Vector<double> & u2);
1060 *  
1061 *   void
1062 *   extract(const DoFHandler<dim> & dof1,
1063 *   const DoFHandler<dim> & dof2,
1064 *   const AffineConstraints<double> &constraints,
1065 *   const Vector<double> & solution,
1066 *   Vector<double> & u2);
1067 *  
1068 *  
1069 *  
1070 *   /*The following FEValues objects are unique_ptrs to 1) avoid default
1071 *   constructors for these objects, and 2) automate memory management*/
1072 *   std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values;
1073 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values;
1074 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor;
1075 *   std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values;
1076 *  
1077 *   std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values_forward;
1078 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_forward;
1079 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor_forward;
1080 *   std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values_forward;
1081 *   using FaceIntegrals =
1082 *   typename std::map<typename DoFHandler<dim>::face_iterator, double>;
1083 *  
1084 *   unsigned int
1085 *   solve_primal_problem();
1086 *  
1087 *   unsigned int
1088 *   solve_dual_problem();
1089 *  
1090 *   void
1091 *   normalize_solutions(Vector<double> &primal_solution,
1092 *   Vector<double> &dual_weights);
1093 *  
1094 *   double
1095 *   get_global_QoI_error(Vector<double> &dual_solution,
1096 *   Vector<double> &error_indicators);
1097 *  
1098 *   void
1099 *   initialize_error_estimation_data();
1100 *  
1101 *   void
1102 *   estimate_on_one_cell(
1103 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1104 *   const Vector<double> & primal_solution,
1105 *   const Vector<double> & dual_weights,
1106 *   const double & lambda_h,
1107 *   Vector<double> & error_indicators,
1108 *   FaceIntegrals & face_integrals);
1109 *  
1110 *   void
1111 *   integrate_over_cell(
1112 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1113 *   const Vector<double> & primal_solution,
1114 *   const Vector<double> & dual_weights,
1115 *   const double & lambda_h,
1116 *   Vector<double> & error_indicators);
1117 *  
1118 *   void
1119 *   integrate_over_regular_face(
1120 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1121 *   const unsigned int & face_no,
1122 *   const Vector<double> & primal_solution,
1123 *   const Vector<double> & dual_weights,
1124 *   FaceIntegrals & face_integrals);
1125 *  
1126 *   void
1127 *   integrate_over_irregular_face(
1128 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1129 *   const unsigned int & face_no,
1130 *   const Vector<double> & primal_solution,
1131 *   const Vector<double> & dual_weights,
1132 *   FaceIntegrals & face_integrals);
1133 *   };
1134 *  
1135 *  
1139 *   template <int dim, bool report_dual>
1140 *   DualWeightedResidual<dim, report_dual>::DualWeightedResidual(
1141 *   const std::string & prm_file,
1142 *   Triangulation<dim> &triangulation,
1143 *   const unsigned int &min_primal_degree,
1144 *   const unsigned int &max_primal_degree,
1145 *   const unsigned int &starting_primal_degree)
1146 *   : Base<dim>(prm_file, triangulation)
1147 *   , PrimalSolver<dim>(prm_file,
1148 *   triangulation,
1149 *   min_primal_degree,
1150 *   max_primal_degree,
1151 *   starting_primal_degree)
1152 *   , DualSolver<dim>(prm_file,
1153 *   triangulation,
1154 *   min_primal_degree + 1,
1155 *   max_primal_degree + 1,
1156 *   starting_primal_degree + 1)
1157 *   {
1158 *   initialize_error_estimation_data();
1159 *   }
1160 *  
1161 *  
1165 *   template <int dim, bool report_dual>
1166 *   DoFHandler<dim> *
1167 *   DualWeightedResidual<dim, report_dual>::get_DoFHandler()
1168 *   {
1169 *   if (!report_dual)
1170 *   return &(PrimalSolver<dim>::dof_handler);
1171 *   else
1172 *   return &(DualSolver<dim>::dof_handler);
1173 *   }
1174 *  
1175 * @endcode
1176 *
1177 * See above function, but to specifically output the primal DoFHandler...
1178 *
1179 * @code
1180 *   template <int dim, bool report_dual>
1181 *   DoFHandler<dim> *
1182 *   DualWeightedResidual<dim, report_dual>::get_primal_DoFHandler()
1183 *   {
1184 *   return &(PrimalSolver<dim>::dof_handler);
1185 *   }
1186 *  
1187 * @endcode
1188 *
1189 * See above function, but for the FECollection
1190 *
1191 * @code
1192 *   template <int dim, bool report_dual>
1193 *   hp::FECollection<dim> *
1194 *   DualWeightedResidual<dim, report_dual>::get_FECollection()
1195 *   {
1196 *   if (!report_dual)
1197 *   return &*(PrimalSolver<dim>::fe_collection);
1198 *   else
1199 *   return &*(DualSolver<dim>::fe_collection);
1200 *   }
1201 *  
1202 * @endcode
1203 *
1204 * See above function, but for the primal FECollection
1205 *
1206 * @code
1207 *   template <int dim, bool report_dual>
1208 *   hp::FECollection<dim> *
1209 *   DualWeightedResidual<dim, report_dual>::get_primal_FECollection()
1210 *   {
1211 *   return &*(PrimalSolver<dim>::fe_collection);
1212 *   }
1213 *  
1214 *   template <int dim, bool report_dual>
1215 *   DoFHandler<dim> *
1216 *   DualWeightedResidual<dim, report_dual>::get_dual_DoFHandler()
1217 *   {
1218 *   return &(DualSolver<dim>::dof_handler);
1219 *   }
1220 *  
1221 *  
1222 *   template <int dim, bool report_dual>
1223 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1224 *   DualWeightedResidual<dim, report_dual>::get_eigenfunctions()
1225 *   {
1226 *   if (!report_dual)
1227 *   return (PrimalSolver<dim>::eigenfunctions);
1228 *   else
1229 *   return (DualSolver<dim>::eigenfunctions);
1230 *   }
1231 *  
1232 *  
1233 *   template <int dim, bool report_dual>
1234 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1235 *   DualWeightedResidual<dim, report_dual>::get_primal_eigenfunctions()
1236 *   {
1237 *   return (PrimalSolver<dim>::eigenfunctions);
1238 *   }
1239 *  
1240 *  
1241 *   template <int dim, bool report_dual>
1242 *   std::unique_ptr<std::vector<double>> &
1243 *   DualWeightedResidual<dim, report_dual>::get_primal_eigenvalues()
1244 *   {
1245 *   return PrimalSolver<dim>::eigenvalues;
1246 *   }
1247 *  
1248 *  
1249 *   template <int dim, bool report_dual>
1250 *   std::unique_ptr<std::vector<double>> &
1251 *   DualWeightedResidual<dim, report_dual>::get_dual_eigenvalues()
1252 *   {
1253 *   return DualSolver<dim>::eigenvalues;
1254 *   }
1255 *  
1256 *   template <int dim, bool report_dual>
1257 *   void
1258 *   DualWeightedResidual<dim, report_dual>::output_solution()
1259 *   {
1260 *   PrimalSolver<dim>::output_solution();
1261 *   }
1262 *  
1263 * @endcode
1264 *
1265 * Solves the primal problem
1266 *
1267 * @code
1268 *   template <int dim, bool report_dual>
1269 *   unsigned int
1270 *   DualWeightedResidual<dim, report_dual>::solve_primal_problem()
1271 *   {
1272 *   return PrimalSolver<dim>::solve_problem();
1273 *   }
1274 *  
1275 * @endcode
1276 *
1277 * Solves the dual problem
1278 *
1279 * @code
1280 *   template <int dim, bool report_dual>
1281 *   unsigned int
1282 *   DualWeightedResidual<dim, report_dual>::solve_dual_problem()
1283 *   {
1284 *   return DualSolver<dim>::solve_problem();
1285 *   }
1286 *  
1287 *  
1291 *   template <int dim, bool report_dual>
1292 *   unsigned int
1293 *   DualWeightedResidual<dim, report_dual>::solve_problem()
1294 *   {
1295 *   DualWeightedResidual<dim, report_dual>::solve_primal_problem();
1296 *   return DualWeightedResidual<dim, report_dual>::solve_dual_problem();
1297 *   }
1298 *  
1299 *  
1302 *   template <int dim, bool report_dual>
1303 *   unsigned int
1304 *   DualWeightedResidual<dim, report_dual>::n_dofs() const
1305 *   {
1306 *   return PrimalSolver<dim>::n_dofs();
1307 *   }
1308 *  
1309 *  
1315 *   template <int dim, bool report_dual>
1316 *   void
1317 *   DualWeightedResidual<dim, report_dual>::synchronize_discretization()
1318 *   {
1319 *   /*Note: No additional checks need to be made ensuring that these operations
1320 *   are legal as these checks are made prior to entering this function (i.e.,
1321 *   if the primal attains a degree N,
1322 *   then, by construction, a degree of N+1 must be permissible for the
1323 *   dual)*/
1324 *   DoFHandler<dim> *dof1 = &(PrimalSolver<dim>::dof_handler);
1325 *   DoFHandler<dim> *dof2 = &(DualSolver<dim>::dof_handler);
1326 *  
1327 *   if (report_dual)
1328 *   {
1329 * @endcode
1330 *
1331 * In this case, we have modified the polynomial orders for the dual;
1332 * need to update the primal
1333 *
1334 * @code
1335 *   dof1 = &(DualSolver<dim>::dof_handler);
1336 *   dof2 = &(PrimalSolver<dim>::dof_handler);
1337 *   }
1338 *   typename DoFHandler<dim>::active_cell_iterator cell1 = dof1->begin_active(),
1339 *   endc1 = dof1->end();
1340 *   typename DoFHandler<dim>::active_cell_iterator cell2 = dof2->begin_active();
1341 *   for (; cell1 < endc1; ++cell1, ++cell2)
1342 *   {
1343 *   cell2->set_active_fe_index(cell1->active_fe_index());
1344 *   }
1345 *   }
1346 *  
1347 *  
1351 *   template <int dim, bool report_dual>
1352 *   void
1353 *   DualWeightedResidual<dim, report_dual>::initialize_error_estimation_data()
1354 *   {
1355 * @endcode
1356 *
1357 * initialize the cell fe_values...
1358 *
1359 * @code
1360 *   cell_hp_fe_values = std::make_unique<hp::FEValues<dim>>(
1361 *   *DualSolver<dim>::fe_collection,
1362 *   *DualSolver<dim>::quadrature_collection,
1363 *   update_values | update_hessians | update_quadrature_points |
1364 *   update_JxW_values);
1365 *   face_hp_fe_values = std::make_unique<hp::FEFaceValues<dim>>(
1366 *   *DualSolver<dim>::fe_collection,
1367 *   *DualSolver<dim>::face_quadrature_collection,
1368 *   update_values | update_gradients | update_JxW_values |
1369 *   update_normal_vectors);
1370 *   face_hp_fe_values_neighbor = std::make_unique<hp::FEFaceValues<dim>>(
1371 *   *DualSolver<dim>::fe_collection,
1372 *   *DualSolver<dim>::face_quadrature_collection,
1373 *   update_values | update_gradients | update_JxW_values |
1374 *   update_normal_vectors);
1375 *   subface_hp_fe_values = std::make_unique<hp::FESubfaceValues<dim>>(
1376 *   *DualSolver<dim>::fe_collection,
1377 *   *DualSolver<dim>::face_quadrature_collection,
1378 *   update_gradients);
1379 *   }
1380 *  
1381 *  
1386 *   template <int dim, bool report_dual>
1387 *   void
1388 *   DualWeightedResidual<dim, report_dual>::normalize_solutions(
1389 *   Vector<double> &primal_solution,
1390 *   Vector<double> &dual_weights)
1391 *   {
1392 *   double sum_primal = 0.0, sum_dual = 0.0;
1393 *   for (const auto &cell :
1394 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1395 *   {
1396 *   cell_hp_fe_values->reinit(cell);
1397 *  
1398 * @endcode
1399 *
1400 * grab the fe_values object
1401 *
1402 * @code
1403 *   const FEValues<dim> &fe_values =
1404 *   cell_hp_fe_values->get_present_fe_values();
1405 *  
1406 *   std::vector<Vector<double>> cell_primal_values(
1407 *   fe_values.n_quadrature_points, Vector<double>(dim)),
1408 *   cell_dual_values(fe_values.n_quadrature_points, Vector<double>(dim));
1409 *   fe_values.get_function_values(primal_solution, cell_primal_values);
1410 *   fe_values.get_function_values(dual_weights, cell_dual_values);
1411 *  
1412 *  
1413 *   for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1414 *   {
1415 *   sum_primal +=
1416 *   cell_primal_values[p] * cell_primal_values[p] * fe_values.JxW(p);
1417 *   sum_dual +=
1418 *   cell_dual_values[p] * cell_dual_values[p] * fe_values.JxW(p);
1419 *   }
1420 *   }
1421 *  
1422 *   primal_solution /= sqrt(sum_primal);
1423 *   dual_weights /= sqrt(sum_dual);
1424 *   }
1425 *  
1426 *  
1430 *   template <int dim, bool report_dual>
1431 *   void
1432 *   DualWeightedResidual<dim, report_dual>::estimate_error(
1433 *   Vector<double> &error_indicators)
1434 *   {
1435 * @endcode
1436 *
1437 * The constraints could be grabbed directly, but this is simple
1438 *
1439 * @code
1440 *   AffineConstraints<double> primal_hanging_node_constraints;
1441 *   DoFTools::make_hanging_node_constraints(PrimalSolver<dim>::dof_handler,
1442 *   primal_hanging_node_constraints);
1443 *   primal_hanging_node_constraints.close();
1444 *  
1445 *   AffineConstraints<double> dual_hanging_node_constraints;
1446 *   DoFTools::make_hanging_node_constraints(DualSolver<dim>::dof_handler,
1447 *   dual_hanging_node_constraints);
1448 *   dual_hanging_node_constraints.close();
1449 *  
1450 * @endcode
1451 *
1452 * First map the primal solution to the space of the dual solution
1453 * This allows us to use just one set of FEValues objects (rather than one
1454 * set for the primal, one for dual)
1455 *
1456
1457 *
1458 *
1459 * @code
1460 *   Vector<double> primal_solution(DualSolver<dim>::dof_handler.n_dofs());
1461 *  
1462 *   embed(PrimalSolver<dim>::dof_handler,
1463 *   DualSolver<dim>::dof_handler,
1464 *   dual_hanging_node_constraints,
1465 *   *(PrimalSolver<dim>::get_solution()),
1466 *   primal_solution);
1467 *  
1468 *   Vector<double> &dual_solution = *(DualSolver<dim>::get_solution());
1469 *  
1470 *   normalize_solutions(primal_solution, dual_solution);
1471 *  
1472 *   Vector<double> dual_weights(DualSolver<dim>::dof_handler.n_dofs()),
1473 *   dual_weights_interm(PrimalSolver<dim>::dof_handler.n_dofs());
1474 *  
1475 * @endcode
1476 *
1477 * First extract the dual solution to the space of the primal
1478 *
1479 * @code
1480 *   extract(DualSolver<dim>::dof_handler,
1481 *   PrimalSolver<dim>::dof_handler,
1482 *   primal_hanging_node_constraints,
1483 *   *(DualSolver<dim>::get_solution()),
1484 *   dual_weights_interm);
1485 *  
1486 * @endcode
1487 *
1488 * Now embed this back to the space of the dual solution
1489 *
1490 * @code
1491 *   embed(PrimalSolver<dim>::dof_handler,
1492 *   DualSolver<dim>::dof_handler,
1493 *   dual_hanging_node_constraints,
1494 *   dual_weights_interm,
1495 *   dual_weights);
1496 *  
1497 *  
1498 * @endcode
1499 *
1500 * Subtract this from the full dual solution
1501 *
1502 * @code
1503 *   dual_weights -= *(DualSolver<dim>::get_solution());
1504 *   dual_weights *= -1.0;
1505 *  
1506 *   *(DualSolver<dim>::get_solution()) -= primal_solution;
1507 *  
1508 *   FaceIntegrals face_integrals;
1509 *   for (const auto &cell :
1510 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1511 *   for (const auto &face : cell->face_iterators())
1512 *   face_integrals[face] = -1e20;
1513 *  
1514 *  
1515 *   for (const auto &cell :
1516 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1517 *   {
1518 *   estimate_on_one_cell(cell,
1519 *   primal_solution,
1520 *   dual_weights,
1521 *   *(PrimalSolver<dim>::get_lambda_h()),
1522 *   error_indicators,
1523 *   face_integrals);
1524 *   }
1525 *   unsigned int present_cell = 0;
1526 *   for (const auto &cell :
1527 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1528 *   {
1529 *   for (const auto &face : cell->face_iterators())
1530 *   {
1531 *   Assert(face_integrals.find(face) != face_integrals.end(),
1532 *   ExcInternalError());
1533 *   error_indicators(present_cell) -= 0.5 * face_integrals[face];
1534 *   }
1535 *   ++present_cell;
1536 *   }
1537 *  
1538 * @endcode
1539 *
1540 * Now, with the error indicators computed, let us produce the
1541 * estimate of the QoI error
1542 *
1543 * @code
1544 *   this->qoi_error_estimate =
1545 *   this->get_global_QoI_error(*(DualSolver<dim>::get_solution()),
1546 *   error_indicators);
1547 *   std::cout << "Estimated QoI error: " << std::setprecision(20)
1548 *   << qoi_error_estimate << std::endl;
1549 *   }
1550 *  
1551 *  
1552 *  
1555 *   template <int dim, bool report_dual>
1556 *   void
1557 *   DualWeightedResidual<dim, report_dual>::estimate_on_one_cell(
1558 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1559 *   const Vector<double> & primal_solution,
1560 *   const Vector<double> & dual_weights,
1561 *   const double & lambda_h,
1562 *   Vector<double> & error_indicators,
1563 *   FaceIntegrals & face_integrals)
1564 *   {
1565 *   integrate_over_cell(
1566 *   cell, primal_solution, dual_weights, lambda_h, error_indicators);
1567 *   for (unsigned int face_no : GeometryInfo<dim>::face_indices())
1568 *   {
1569 *   if (cell->face(face_no)->at_boundary())
1570 *   {
1571 *   face_integrals[cell->face(face_no)] = 0.0;
1572 *   continue;
1573 *   }
1574 *   if ((cell->neighbor(face_no)->has_children() == false) &&
1575 *   (cell->neighbor(face_no)->level() == cell->level()) &&
1576 *   (cell->neighbor(face_no)->index() < cell->index()))
1577 *   continue;
1578 *   if (cell->at_boundary(face_no) == false)
1579 *   if (cell->neighbor(face_no)->level() < cell->level())
1580 *   continue;
1581 *   if (cell->face(face_no)->has_children() == false)
1582 *   integrate_over_regular_face(
1583 *   cell, face_no, primal_solution, dual_weights, face_integrals);
1584 *   else
1585 *   integrate_over_irregular_face(
1586 *   cell, face_no, primal_solution, dual_weights, face_integrals);
1587 *   }
1588 *   }
1589 *  
1590 *  
1593 *   template <int dim, bool report_dual>
1594 *   void
1595 *   DualWeightedResidual<dim, report_dual>::integrate_over_cell(
1596 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1597 *   const Vector<double> & primal_solution,
1598 *   const Vector<double> & dual_weights,
1599 *   const double & lambda_h,
1600 *   Vector<double> & error_indicators)
1601 *   {
1602 *   cell_hp_fe_values->reinit(cell);
1603 * @endcode
1604 *
1605 * Grab the fe_values object
1606 *
1607 * @code
1608 *   const FEValues<dim> &fe_values = cell_hp_fe_values->get_present_fe_values();
1609 *   std::vector<std::vector<Tensor<2, dim, double>>> cell_hessians(
1610 *   fe_values.n_quadrature_points, std::vector<Tensor<2, dim, double>>(dim));
1611 *   std::vector<Vector<double>> cell_primal_values(
1612 *   fe_values.n_quadrature_points, Vector<double>(dim)),
1613 *   cell_dual_values(fe_values.n_quadrature_points, Vector<double>(dim));
1614 *   fe_values.get_function_values(primal_solution, cell_primal_values);
1615 *   fe_values.get_function_hessians(primal_solution, cell_hessians);
1616 *   fe_values.get_function_values(dual_weights, cell_dual_values);
1617 *  
1618 *  
1619 *  
1620 *   double sum = 0.0;
1621 *   for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1622 *   {
1623 *   sum +=
1624 *   (/*x-component*/ (cell_hessians[p][1][1][0] -
1625 *   cell_hessians[p][0][1][1]) *
1626 *   (cell_dual_values[p](0)) +
1627 *   /*y-component*/
1628 *   (cell_hessians[p][0][0][1] - cell_hessians[p][1][0][0]) *
1629 *   (cell_dual_values[p](1)) -
1630 *   lambda_h * (cell_primal_values[p](0) * cell_dual_values[p](0) +
1631 *   cell_primal_values[p](1) * cell_dual_values[p](1))) *
1632 *   fe_values.JxW(p);
1633 *   }
1634 *  
1635 *   error_indicators(cell->active_cell_index()) += sum;
1636 *   }
1637 *  
1638 *  
1641 *   template <int dim, bool report_dual>
1642 *   void
1643 *   DualWeightedResidual<dim, report_dual>::integrate_over_regular_face(
1644 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1645 *   const unsigned int & face_no,
1646 *   const Vector<double> & primal_solution,
1647 *   const Vector<double> & dual_weights,
1648 *   FaceIntegrals & face_integrals)
1649 *   {
1650 *   Assert(cell->neighbor(face_no).state() == IteratorState::valid,
1651 *   ExcInternalError());
1652 *   const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1653 *   const auto neighbor = cell->neighbor(face_no);
1654 *  
1655 *   const unsigned int quadrature_index =
1656 *   std::max(cell->active_fe_index(), neighbor->active_fe_index());
1657 *   face_hp_fe_values->reinit(cell, face_no, quadrature_index);
1658 *   const FEFaceValues<dim> &fe_face_values_cell =
1659 *   face_hp_fe_values->get_present_fe_values();
1660 *   std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1661 *   fe_face_values_cell.n_quadrature_points,
1662 *   std::vector<Tensor<1, dim, double>>(dim)),
1663 *   neighbor_primal_grads(fe_face_values_cell.n_quadrature_points,
1664 *   std::vector<Tensor<1, dim, double>>(dim));
1665 *   fe_face_values_cell.get_function_gradients(primal_solution,
1666 *   cell_primal_grads);
1667 *  
1668 *   face_hp_fe_values_neighbor->reinit(neighbor,
1669 *   neighbor_neighbor,
1670 *   quadrature_index);
1671 *   const FEFaceValues<dim> &fe_face_values_cell_neighbor =
1672 *   face_hp_fe_values_neighbor->get_present_fe_values();
1673 *   fe_face_values_cell_neighbor.get_function_gradients(primal_solution,
1674 *   neighbor_primal_grads);
1675 *   const unsigned int n_q_points = fe_face_values_cell.n_quadrature_points;
1676 *   double face_integral = 0.0;
1677 *   std::vector<Vector<double>> cell_dual_values(n_q_points,
1678 *   Vector<double>(dim));
1679 *   fe_face_values_cell.get_function_values(dual_weights, cell_dual_values);
1680 *   for (unsigned int p = 0; p < n_q_points; ++p)
1681 *   {
1682 *   auto face_normal = fe_face_values_cell.normal_vector(p);
1683 *  
1684 *   face_integral +=
1685 *   (cell_primal_grads[p][1][0] - cell_primal_grads[p][0][1] -
1686 *   neighbor_primal_grads[p][1][0] + neighbor_primal_grads[p][0][1]) *
1687 *   (cell_dual_values[p][0] * face_normal[1] -
1688 *   cell_dual_values[p][1] * face_normal[0]) *
1689 *   fe_face_values_cell.JxW(p);
1690 *   }
1691 *   Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
1692 *   ExcInternalError());
1693 *   Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
1694 *   face_integrals[cell->face(face_no)] = face_integral;
1695 *   }
1696 *  
1697 *  
1700 *   template <int dim, bool report_dual>
1701 *   void
1702 *   DualWeightedResidual<dim, report_dual>::integrate_over_irregular_face(
1703 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1704 *   const unsigned int & face_no,
1705 *   const Vector<double> & primal_solution,
1706 *   const Vector<double> & dual_weights,
1707 *   FaceIntegrals & face_integrals)
1708 *   {
1709 *   const typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
1710 *   const typename DoFHandler<dim>::cell_iterator neighbor =
1711 *   cell->neighbor(face_no);
1712 *  
1713 *   Assert(neighbor.state() == IteratorState::valid, ExcInternalError());
1714 *   Assert(neighbor->has_children(), ExcInternalError());
1715 *   (void)neighbor;
1716 *   const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1717 *   for (unsigned int subface_no = 0; subface_no < face->n_children();
1718 *   ++subface_no)
1719 *   {
1720 *   const typename DoFHandler<dim>::active_cell_iterator neighbor_child =
1721 *   cell->neighbor_child_on_subface(face_no, subface_no);
1722 *   Assert(neighbor_child->face(neighbor_neighbor) ==
1723 *   cell->face(face_no)->child(subface_no),
1724 *   ExcInternalError());
1725 *   const unsigned int quadrature_index =
1726 *   std::max(cell->active_fe_index(), neighbor_child->active_fe_index());
1727 * @endcode
1728 *
1729 * initialize fe_subface values_cell
1730 *
1731 * @code
1732 *   subface_hp_fe_values->reinit(cell,
1733 *   face_no,
1734 *   subface_no,
1735 *   quadrature_index);
1736 *   const FESubfaceValues<dim> &subface_fe_values_cell =
1737 *   subface_hp_fe_values->get_present_fe_values();
1738 *   std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1739 *   subface_fe_values_cell.n_quadrature_points,
1740 *   std::vector<Tensor<1, dim, double>>(dim)),
1741 *   neighbor_primal_grads(subface_fe_values_cell.n_quadrature_points,
1742 *   std::vector<Tensor<1, dim, double>>(dim));
1743 *   subface_fe_values_cell.get_function_gradients(primal_solution,
1744 *   cell_primal_grads);
1745 * @endcode
1746 *
1747 * initialize fe_face_values_neighbor
1748 *
1749 * @code
1750 *   face_hp_fe_values_neighbor->reinit(neighbor_child,
1751 *   neighbor_neighbor,
1752 *   quadrature_index);
1753 *   const FEFaceValues<dim> &face_fe_values_neighbor =
1754 *   face_hp_fe_values_neighbor->get_present_fe_values();
1755 *   face_fe_values_neighbor.get_function_gradients(primal_solution,
1756 *   neighbor_primal_grads);
1757 *   const unsigned int n_q_points =
1758 *   subface_fe_values_cell.n_quadrature_points;
1759 *   std::vector<Vector<double>> cell_dual_values(n_q_points,
1760 *   Vector<double>(dim));
1761 *   face_fe_values_neighbor.get_function_values(dual_weights,
1762 *   cell_dual_values);
1763 *  
1764 *   double face_integral = 0.0;
1765 *  
1766 *   for (unsigned int p = 0; p < n_q_points; ++p)
1767 *   {
1768 *   auto face_normal = face_fe_values_neighbor.normal_vector(p);
1769 *   face_integral +=
1770 *   (cell_primal_grads[p][0][1] - cell_primal_grads[p][1][0] +
1771 *   neighbor_primal_grads[p][1][0] -
1772 *   neighbor_primal_grads[p][0][1]) *
1773 *   (cell_dual_values[p][0] * face_normal[1] -
1774 *   cell_dual_values[p][1] * face_normal[0]) *
1775 *   face_fe_values_neighbor.JxW(p);
1776 *   }
1777 *   face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
1778 *   }
1779 *   double sum = 0.0;
1780 *   for (unsigned int subface_no = 0; subface_no < face->n_children();
1781 *   ++subface_no)
1782 *   {
1783 *   Assert(face_integrals.find(face->child(subface_no)) !=
1784 *   face_integrals.end(),
1785 *   ExcInternalError());
1786 *   Assert(face_integrals[face->child(subface_no)] != -1e20,
1787 *   ExcInternalError());
1788 *   sum += face_integrals[face->child(subface_no)];
1789 *   }
1790 *   face_integrals[face] = sum;
1791 *   }
1792 *  
1793 *   template <int dim, bool report_dual>
1794 *   double
1795 *   DualWeightedResidual<dim, report_dual>::get_global_QoI_error(
1796 *   Vector<double> &dual_solution,
1797 *   Vector<double> &error_indicators)
1798 *   {
1799 *   auto dual_less_primal =
1800 *   dual_solution; // Note: We have already extracted the primal solution...
1801 *  
1802 *  
1803 *   double scaling_factor = 0.0;
1804 *   for (const auto &cell :
1805 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1806 *   {
1807 *   cell_hp_fe_values->reinit(cell);
1808 * @endcode
1809 *
1810 * grab the fe_values object
1811 *
1812 * @code
1813 *   const FEValues<dim> &fe_values =
1814 *   cell_hp_fe_values->get_present_fe_values();
1815 *  
1816 *   std::vector<Vector<double>> cell_values(fe_values.n_quadrature_points,
1817 *   Vector<double>(dim));
1818 *   fe_values.get_function_values(dual_less_primal, cell_values);
1819 *  
1820 *   for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1821 *   {
1822 *   scaling_factor +=
1823 *   (cell_values[p] * cell_values[p]) * fe_values.JxW(p);
1824 *   }
1825 *   }
1826 *   double global_QoI_error = 0.0;
1827 *   for (const auto &indicator : error_indicators)
1828 *   {
1829 *   global_QoI_error += indicator;
1830 *   }
1831 *  
1832 *   global_QoI_error /= (1 - 0.5 * scaling_factor);
1833 *   return global_QoI_error;
1834 *   }
1835 *  
1836 *  
1837 *   template <int dim, bool report_dual>
1838 *   void
1839 *   DualWeightedResidual<dim, report_dual>::embed(
1840 *   const DoFHandler<dim> & dof1,
1841 *   const DoFHandler<dim> & dof2,
1842 *   const AffineConstraints<double> &constraints,
1843 *   const Vector<double> & solution,
1844 *   Vector<double> & u2)
1845 *   {
1846 *   assert(u2.size() == dof2.n_dofs() && "Incorrect input vector size!");
1847 *  
1848 *   u2 = 0.0;
1849 *  
1850 *   typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
1851 *   endc1 = dof1.end();
1852 *   typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
1853 *  
1854 *   for (; cell1 < endc1; ++cell1, ++cell2)
1855 *   {
1856 *   const auto &fe1 =
1857 *   dynamic_cast<const FE_Nedelec<dim> &>(cell1->get_fe());
1858 *   const auto &fe2 =
1859 *   dynamic_cast<const FE_Nedelec<dim> &>(cell2->get_fe());
1860 *  
1861 *   assert(fe1.degree < fe2.degree && "Incorrect usage of embed!");
1862 *  
1863 * @endcode
1864 *
1865 * Get the embedding_dofs
1866 *
1867
1868 *
1869 *
1870
1871 *
1872 *
1873 * @code
1874 *   std::vector<unsigned int> embedding_dofs =
1875 *   fe2.get_embedding_dofs(fe1.degree);
1876 *   const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1877 *  
1878 *  
1879 *   Vector<double> local_dof_values_1;
1880 *   Vector<double> local_dof_values_2(dofs_per_cell2);
1881 *  
1882 *   local_dof_values_1.reinit(fe1.dofs_per_cell);
1883 *   cell1->get_dof_values(solution, local_dof_values_1);
1884 *  
1885 *   for (unsigned int i = 0; i < local_dof_values_1.size(); ++i)
1886 *   local_dof_values_2[embedding_dofs[i]] = local_dof_values_1[i];
1887 *  
1888 * @endcode
1889 *
1890 * Now set this changes to the global vector
1891 *
1892 * @code
1893 *   cell2->set_dof_values(local_dof_values_2, u2);
1894 *   }
1895 *  
1896 *   u2.compress(VectorOperation::insert);
1897 * @endcode
1898 *
1899 * Applies the constraints of the target finite element space
1900 *
1901 * @code
1902 *   constraints.distribute(u2);
1903 *   }
1904 *  
1905 *   template <int dim, bool report_dual>
1906 *   void
1907 *   DualWeightedResidual<dim, report_dual>::extract(
1908 *   const DoFHandler<dim> & dof1,
1909 *   const DoFHandler<dim> & dof2,
1910 *   const AffineConstraints<double> &constraints,
1911 *   const Vector<double> & solution,
1912 *   Vector<double> & u2)
1913 *   {
1914 * @endcode
1915 *
1916 * Maps from fe1 to fe2
1917 *
1918 * @code
1919 *   assert(u2.size() == dof2.n_dofs() && "Incorrect input vector size!");
1920 *  
1921 *   u2 = 0.0;
1922 *  
1923 *   typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
1924 *   endc1 = dof1.end();
1925 *   typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
1926 *  
1927 *   for (; cell1 < endc1; ++cell1, ++cell2)
1928 *   {
1929 *   const auto &fe1 =
1930 *   dynamic_cast<const FE_Nedelec<dim> &>(cell1->get_fe());
1931 *   const auto &fe2 =
1932 *   dynamic_cast<const FE_Nedelec<dim> &>(cell2->get_fe());
1933 *  
1934 *   assert(fe1.degree > fe2.degree && "Incorrect usage of extract!");
1935 *  
1936 * @endcode
1937 *
1938 * Get the embedding_dofs
1939 *
1940 * @code
1941 *   std::vector<unsigned int> embedding_dofs =
1942 *   fe1.get_embedding_dofs(fe2.degree);
1943 *   const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1944 *  
1945 *  
1946 *   Vector<double> local_dof_values_1;
1947 *   Vector<double> local_dof_values_2(dofs_per_cell2);
1948 *  
1949 *   local_dof_values_1.reinit(fe1.dofs_per_cell);
1950 *   cell1->get_dof_values(solution, local_dof_values_1);
1951 *  
1952 *   for (unsigned int i = 0; i < local_dof_values_2.size(); ++i)
1953 *   local_dof_values_2[i] = local_dof_values_1[embedding_dofs[i]];
1954 *  
1955 * @endcode
1956 *
1957 * Now set this changes to the global vector
1958 *
1959 * @code
1960 *   cell2->set_dof_values(local_dof_values_2, u2);
1961 *   }
1962 *  
1963 *   u2.compress(VectorOperation::insert);
1964 * @endcode
1965 *
1966 * Applies the constraints of the target finite element space
1967 *
1968 * @code
1969 *   constraints.distribute(u2);
1970 *   }
1971 *   template <int dim, bool report_dual>
1972 *   void
1973 *   DualWeightedResidual<dim, report_dual>::output_eigenvalue_data(
1974 *   std::ofstream &os)
1975 *   {
1976 *   os << (*this->get_primal_eigenvalues())[0] << " "
1977 *   << (this->get_primal_DoFHandler())->n_dofs() << " "
1978 *   << (*this->get_dual_eigenvalues())[0] << " "
1979 *   << (this->get_dual_DoFHandler())->n_dofs() << std::endl;
1980 *   }
1981 *   template <int dim, bool report_dual>
1982 *   void
1983 *   DualWeightedResidual<dim, report_dual>::output_qoi_error_estimates(
1984 *   std::ofstream &os)
1985 *   {
1986 *   os << qoi_error_estimate << std::endl;
1987 *   }
1988 *  
1989 *  
1993 *   template <int dim>
1994 *   class KellyErrorIndicator : public PrimalSolver<dim>
1995 *   {
1996 *   public:
1997 *   std::string
1998 *   name() const
1999 *   {
2000 *   return "Kelly";
2001 *   }
2002 *   void
2003 *   output_eigenvalue_data(std::ofstream &os);
2004 *   void
2005 *   output_qoi_error_estimates(std::ofstream &);
2006 *   KellyErrorIndicator(const std::string & prm_file,
2007 *   Triangulation<dim> &coarse_grid,
2008 *   const unsigned int &min_degree,
2009 *   const unsigned int &max_degree,
2010 *   const unsigned int &starting_degree);
2011 *  
2012 *   virtual unsigned int
2013 *   solve_problem() override;
2014 *  
2015 *   virtual void
2016 *   output_solution() override;
2017 *  
2018 *   hp::FECollection<dim> *
2019 *   get_FECollection();
2020 *  
2021 *   hp::FECollection<dim> *
2022 *   get_primal_FECollection();
2023 *  
2024 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2025 *   get_eigenfunctions();
2026 *  
2027 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2028 *   get_primal_eigenfunctions();
2029 *  
2030 *   std::unique_ptr<std::vector<double>> &
2031 *   get_primal_eigenvalues();
2032 *  
2033 *  
2034 *   void
2035 *   synchronize_discretization();
2036 *  
2037 *   DoFHandler<dim> *
2038 *   get_DoFHandler();
2039 *  
2040 *   DoFHandler<dim> *
2041 *   get_primal_DoFHandler();
2042 *  
2043 *   unsigned int
2044 *   get_max_degree()
2045 *   {
2046 *   return PrimalSolver<dim>::fe_collection->max_degree();
2047 *   }
2048 *   double qoi_error_estimate = 0;
2049 *  
2050 *   protected:
2051 *   void
2052 *   estimate_error(Vector<double> &error_indicators);
2053 *  
2054 *   private:
2055 *   void
2056 *   prune_eigenpairs(const double &TOL);
2057 *  
2058 *   #if DEAL_II_VERSION_GTE(9, 6, 0)
2059 *   std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
2060 *   #else
2061 *   std::vector<const PETScWrappers::MPI::Vector *> eigenfunction_ptrs;
2062 *   #endif
2063 *   std::vector<const double *> eigenvalue_ptrs;
2064 *  
2065 *   std::vector<std::shared_ptr<Vector<float>>> errors;
2066 *   };
2067 *  
2068 *   template <int dim>
2069 *   KellyErrorIndicator<dim>::KellyErrorIndicator(
2070 *   const std::string & prm_file,
2071 *   Triangulation<dim> &coarse_grid,
2072 *   const unsigned int &min_degree,
2073 *   const unsigned int &max_degree,
2074 *   const unsigned int &starting_degree)
2075 *   : Base<dim>(prm_file, coarse_grid)
2076 *   , PrimalSolver<dim>(prm_file,
2077 *   coarse_grid,
2078 *   min_degree,
2079 *   max_degree,
2080 *   starting_degree)
2081 *   {}
2082 *  
2083 *   template <int dim>
2084 *   unsigned int
2085 *   KellyErrorIndicator<dim>::solve_problem()
2086 *   {
2087 *   return PrimalSolver<dim>::solve_problem();
2088 *   }
2089 *  
2090 *   template <int dim>
2091 *   hp::FECollection<dim> *
2092 *   KellyErrorIndicator<dim>::get_FECollection()
2093 *   {
2094 *   return &*(PrimalSolver<dim>::fe_collection);
2095 *   }
2096 *  
2097 *   template <int dim>
2098 *   hp::FECollection<dim> *
2099 *   KellyErrorIndicator<dim>::get_primal_FECollection()
2100 *   {
2101 *   return &*(PrimalSolver<dim>::fe_collection);
2102 *   }
2103 *  
2104 *   template <int dim>
2105 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2106 *   KellyErrorIndicator<dim>::get_eigenfunctions()
2107 *   {
2108 *   return (PrimalSolver<dim>::eigenfunctions);
2109 *   }
2110 *  
2111 *   template <int dim>
2112 *   std::unique_ptr<std::vector<double>> &
2113 *   KellyErrorIndicator<dim>::get_primal_eigenvalues()
2114 *   {
2115 *   return PrimalSolver<dim>::eigenvalues;
2116 *   }
2117 *  
2118 *   template <int dim>
2119 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2120 *   KellyErrorIndicator<dim>::get_primal_eigenfunctions()
2121 *   {
2122 *   return (PrimalSolver<dim>::eigenfunctions);
2123 *   }
2124 *  
2125 *   template <int dim>
2126 *   DoFHandler<dim> *
2127 *   KellyErrorIndicator<dim>::get_DoFHandler()
2128 *   {
2129 *   return &(PrimalSolver<dim>::dof_handler);
2130 *   }
2131 *  
2132 *   template <int dim>
2133 *   DoFHandler<dim> *
2134 *   KellyErrorIndicator<dim>::get_primal_DoFHandler()
2135 *   {
2136 *   return &(PrimalSolver<dim>::dof_handler);
2137 *   }
2138 *  
2139 *   template <int dim>
2140 *   void
2141 *   KellyErrorIndicator<dim>::synchronize_discretization()
2142 *   {
2143 * @endcode
2144 *
2145 * This function does nothing for this error indicator
2146 *
2147 * @code
2148 *   return;
2149 *   }
2150 *  
2151 *   template <int dim>
2152 *   void
2153 *   KellyErrorIndicator<dim>::output_solution()
2154 *   {
2155 *   PrimalSolver<dim>::output_solution();
2156 *   }
2157 *  
2158 *   template <int dim>
2159 *   void
2160 *   KellyErrorIndicator<dim>::prune_eigenpairs(const double &TOL)
2161 *   {
2162 *   unsigned int count = 0;
2163 *   for (size_t eigenpair_index = 0;
2164 *   eigenpair_index < this->eigenfunctions->size();
2165 *   ++eigenpair_index)
2166 *   {
2167 *   if (count >= this->n_eigenpairs)
2168 *   break;
2169 *   if (abs((*this->eigenvalues)[eigenpair_index]) < TOL)
2170 *   continue;
2171 *  
2172 *   eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
2173 *   eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
2174 *   }
2175 *   }
2176 *  
2177 *   template <int dim>
2178 *   void
2179 *   KellyErrorIndicator<dim>::estimate_error(Vector<double> &error_indicators)
2180 *   {
2181 *   std::cout << "Marking cells via Kelly indicator..." << std::endl;
2182 *   prune_eigenpairs(1e-9);
2183 * @endcode
2184 *
2185 * deallocate the errors vector
2186 *
2187 * @code
2188 *   errors.clear();
2189 *   for (size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2190 *   {
2191 *   errors.emplace_back(
2192 *   new Vector<float>(this->triangulation->n_active_cells()));
2193 *   }
2194 *   std::vector<Vector<float> *> estimated_error_per_cell(
2195 *   eigenfunction_ptrs.size());
2196 *   for (size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2197 *   {
2198 *   estimated_error_per_cell[i] = errors[i].get();
2199 *   }
2200 *  
2201 *   #if DEAL_II_VERSION_GTE(9, 6, 0)
2202 *   const auto solution_view = make_array_view(eigenfunction_ptrs);
2203 *   auto error_view = make_array_view(estimated_error_per_cell);
2204 *   KellyErrorEstimator<dim>::estimate(this->dof_handler,
2205 *   *this->face_quadrature_collection,
2206 *   {},
2207 *   solution_view,
2208 *   error_view);
2209 *   #else
2210 *   KellyErrorEstimator<dim>::estimate(this->dof_handler,
2211 *   *this->face_quadrature_collection,
2212 *   {},
2213 *   eigenfunction_ptrs,
2214 *   estimated_error_per_cell);
2215 *   #endif
2216 *  
2217 *   for (auto &error_vec : errors)
2218 *   {
2219 *   auto normalized_vec = *error_vec;
2220 *   normalized_vec /= normalized_vec.l1_norm();
2221 *  
2222 *   for (unsigned int i = 0; i < error_indicators.size(); ++i)
2223 *   error_indicators(i) += double(normalized_vec(i));
2224 *   }
2225 *   std::cout << "...Done!" << std::endl;
2226 *   }
2227 *   template <int dim>
2228 *   void
2229 *   KellyErrorIndicator<dim>::output_eigenvalue_data(std::ofstream &os)
2230 *   {
2231 *   os << (*this->get_primal_eigenvalues())[0] << " "
2232 *   << (this->get_primal_DoFHandler())->n_dofs() << std::endl;
2233 *   }
2234 *   template <int dim>
2235 *   void
2236 *   KellyErrorIndicator<dim>::output_qoi_error_estimates(std::ofstream &)
2237 *   {
2238 *   return;
2239 *   }
2240 *  
2241 *   } // namespace ErrorIndicators
2242 *  
2243 *  
2246 *   namespace RegularityIndicators
2247 *   {
2248 *   using namespace dealii;
2249 *  
2250 *   /* For the Legendre smoothness indicator*/
2251 *   /* Adapted from M. Fehling's smoothness_estimator.cc*/
2252 *   template <int dim>
2253 *   class LegendreInfo
2254 *   {};
2255 *  
2256 *   template <>
2257 *   class LegendreInfo<2>
2258 *   {
2259 *   public:
2260 *   std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
2261 *  
2262 *   hp::FECollection<2> *fe_collection = nullptr;
2263 *   DoFHandler<2> * dof_handler = nullptr;
2264 *  
2265 *   void
2266 *   initialization()
2267 *   {
2268 *   assert(fe_collection != nullptr && dof_handler != nullptr &&
2269 *   "A valid FECollection and DoFHandler must be accessible!");
2270 *  
2271 *   legendre_u = std::make_unique<FESeries::Legendre<2>>(
2272 *   SmoothnessEstimator::Legendre::default_fe_series(*fe_collection, 0));
2273 *   legendre_v = std::make_unique<FESeries::Legendre<2>>(
2274 *   SmoothnessEstimator::Legendre::default_fe_series(*fe_collection, 1));
2275 *  
2276 *   legendre_u->precalculate_all_transformation_matrices();
2277 *   legendre_v->precalculate_all_transformation_matrices();
2278 *   }
2279 *  
2280 *   template <class VectorType>
2281 *   void
2282 *   compute_coefficient_decay(const VectorType & eigenfunction,
2283 *   std::vector<double> &smoothness_indicators)
2284 *   {
2285 * @endcode
2286 *
2287 * Compute the coefficients for the u and v components of the solution
2288 * separately,
2289 *
2290 * @code
2291 *   Vector<float> smoothness_u(smoothness_indicators.size()),
2292 *   smoothness_v(smoothness_indicators.size());
2293 *  
2294 *   SmoothnessEstimator::Legendre::coefficient_decay(*legendre_u,
2295 *   *dof_handler,
2296 *   eigenfunction,
2297 *   smoothness_u);
2298 *  
2299 *   SmoothnessEstimator::Legendre::coefficient_decay(*legendre_v,
2300 *   *dof_handler,
2301 *   eigenfunction,
2302 *   smoothness_v);
2303 *  
2304 *   for (unsigned int i = 0; i < smoothness_indicators.size(); ++i)
2305 *   {
2306 *   smoothness_indicators[i] = std::min(smoothness_u[i], smoothness_v[i]);
2307 *   }
2308 *   }
2309 *   };
2310 *  
2311 *  
2314 *   template <int dim>
2315 *   class LegendreIndicator
2316 *   {
2317 *   public:
2318 *   void
2319 *   attach_FE_info_and_initialize(hp::FECollection<dim> *fe_ptr,
2320 *   DoFHandler<dim> * dof_ptr);
2321 *  
2322 *   protected:
2323 *   template <class VectorType>
2324 *   void
2325 *   estimate_smoothness(
2326 *   const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2327 *   const unsigned int & index_of_goal,
2328 *   std::vector<double> & smoothness_indicators);
2329 *  
2330 *   private:
2331 *   LegendreInfo<dim> legendre;
2332 *   };
2333 *  
2334 *   template <int dim>
2335 *   void
2336 *   LegendreIndicator<dim>::attach_FE_info_and_initialize(
2337 *   hp::FECollection<dim> *fe_ptr,
2338 *   DoFHandler<dim> * dof_ptr)
2339 *   {
2340 *   legendre.fe_collection = fe_ptr;
2341 *   legendre.dof_handler = dof_ptr;
2342 *   this->legendre.initialization();
2343 *   }
2344 *  
2345 *   template <int dim>
2346 *   template <class VectorType>
2347 *   void
2348 *   LegendreIndicator<dim>::estimate_smoothness(
2349 *   const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2350 *   const unsigned int & index_of_goal,
2351 *   std::vector<double> & smoothness_indicators)
2352 *   {
2353 *   this->legendre.compute_coefficient_decay((*eigenfunctions)[index_of_goal],
2354 *   smoothness_indicators);
2355 *   }
2356 *   } // namespace RegularityIndicators
2357 *  
2358 *  
2362 *   namespace Refinement
2363 *   {
2364 *   using namespace dealii;
2365 *   using namespace Maxwell;
2366 *  
2367 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2368 *   class Refiner : public ErrorIndicator, public RegularityIndicator
2369 *   {
2370 *   public:
2371 *   Refiner(const std::string & prm_file,
2372 *   Triangulation<dim> &coarse_grid,
2373 *   const unsigned int &min_degree,
2374 *   const unsigned int &max_degree,
2375 *   const unsigned int &starting_degree);
2376 *  
2377 *   void
2378 *   execute_refinement(const double &smoothness_threshold_fraction);
2379 *  
2380 *   virtual void
2381 *   output_solution() override;
2382 *  
2383 *   private:
2384 *   Vector<double> estimated_error_per_cell;
2385 *   std::vector<double> smoothness_indicators;
2386 *   std::ofstream eigenvalues_out;
2387 *   std::ofstream error_estimate_out;
2388 *   };
2389 *  
2390 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2391 *   Refiner<dim, ErrorIndicator, RegularityIndicator>::Refiner(
2392 *   const std::string & prm_file,
2393 *   Triangulation<dim> &coarse_grid,
2394 *   const unsigned int &min_degree,
2395 *   const unsigned int &max_degree,
2396 *   const unsigned int &starting_degree)
2397 *   : Base<dim>(prm_file, coarse_grid)
2398 *   , ErrorIndicator(prm_file,
2399 *   coarse_grid,
2400 *   min_degree,
2401 *   max_degree,
2402 *   starting_degree)
2403 *   , RegularityIndicator()
2404 *   {
2405 *   if (ErrorIndicator::name() == "DWR")
2406 *   {
2407 *   error_estimate_out.open("error_estimate.txt");
2408 *   error_estimate_out << std::setprecision(20);
2409 *   }
2410 *  
2411 *   eigenvalues_out.open("eigenvalues_" + ErrorIndicator::name() + "_out.txt");
2412 *   eigenvalues_out << std::setprecision(20);
2413 *   }
2414 *  
2415 * @endcode
2416 *
2417 * For generating samples of the curl of the electric field
2418 *
2419 * @code
2420 *   template <int dim>
2421 *   class CurlPostprocessor : public DataPostprocessorScalar<dim>
2422 *   {
2423 *   public:
2424 *   CurlPostprocessor()
2425 *   : DataPostprocessorScalar<dim>("Curl", update_gradients)
2426 *   {}
2427 *  
2428 *   virtual void
2429 *   evaluate_vector_field(
2430 *   const DataPostprocessorInputs::Vector<dim> &input_data,
2431 *   std::vector<Vector<double>> &computed_quantities) const override
2432 *   {
2433 *   AssertDimension(input_data.solution_gradients.size(),
2434 *   computed_quantities.size());
2435 *   for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
2436 *   {
2437 *   computed_quantities[p](0) = input_data.solution_gradients[p][1][0] -
2438 *   input_data.solution_gradients[p][0][1];
2439 *   }
2440 *   }
2441 *   };
2442 *  
2443 *  
2450 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2451 *   void
2452 *   Refiner<dim, ErrorIndicator, RegularityIndicator>::output_solution()
2453 *   {
2454 *   CurlPostprocessor<dim> curl_u;
2455 *  
2456 *   DataOut<dim> data_out;
2457 *   auto & output_dof = *(ErrorIndicator::get_primal_DoFHandler());
2458 *   data_out.attach_dof_handler(output_dof);
2459 *   Vector<double> fe_degrees(this->triangulation->n_active_cells());
2460 *   for (const auto &cell : output_dof.active_cell_iterators())
2461 *   fe_degrees(cell->active_cell_index()) =
2462 *   (*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
2463 *   .degree;
2464 *   data_out.add_data_vector(fe_degrees, "fe_degree");
2465 *  
2466 *   data_out.add_data_vector(estimated_error_per_cell, "error");
2467 *   Vector<double> smoothness_out(this->triangulation->n_active_cells());
2468 *   for (const auto &cell : output_dof.active_cell_iterators())
2469 *   {
2470 *   auto i = cell->active_cell_index();
2471 *   if (!cell->refine_flag_set() && !cell->coarsen_flag_set())
2472 *   smoothness_out(i) = -1;
2473 *   else
2474 *   smoothness_out(i) = smoothness_indicators[i];
2475 *   }
2476 *   data_out.add_data_vector(smoothness_out, "smoothness");
2477 *   data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2478 *   std::string("eigenfunction_no_") +
2479 *   Utilities::int_to_string(0));
2480 *   data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2481 *   curl_u);
2482 *  
2483 *   ErrorIndicator::output_eigenvalue_data(eigenvalues_out);
2484 *   ErrorIndicator::output_qoi_error_estimates(error_estimate_out);
2485 *  
2486 *   std::cout << "Number of DoFs: " << (this->get_primal_DoFHandler())->n_dofs()
2487 *   << std::endl;
2488 *  
2489 *  
2490 *   data_out.build_patches();
2491 *   std::ofstream output("eigenvectors-" + ErrorIndicator::name() + "-" +
2492 *   std::to_string(this->refinement_cycle) + +".vtu");
2493 *   data_out.write_vtu(output);
2494 *   }
2495 *  
2496 *  
2497 *  
2502 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2503 *   void
2504 *   Refiner<dim, ErrorIndicator, RegularityIndicator>::execute_refinement(
2505 *   const double &smoothness_threshold_fraction)
2506 *   {
2507 * @endcode
2508 *
2509 * First initialize the RegularityIndicator...
2510 * Depending on the limits set, this may take a while
2511 *
2512 * @code
2513 *   std::cout << "Initializing RegularityIndicator..." << std::endl;
2514 *   std::cout
2515 *   << "(This may take a while if the max expansion order is set too high)"
2516 *   << std::endl;
2517 *   RegularityIndicator::attach_FE_info_and_initialize(
2518 *   ErrorIndicator::get_FECollection(), ErrorIndicator::get_DoFHandler());
2519 *   std::cout << "Done!" << std::endl << "Starting Refinement..." << std::endl;
2520 *  
2521 *   for (unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
2522 *   {
2523 *   this->set_refinement_cycle(cycle);
2524 *   std::cout << "Cycle: " << cycle << std::endl;
2525 *   ErrorIndicator::solve_problem();
2526 *   this->estimated_error_per_cell.reinit(
2527 *   this->triangulation->n_active_cells());
2528 *  
2529 *   ErrorIndicator::estimate_error(estimated_error_per_cell);
2530 *  
2531 * @endcode
2532 *
2533 * Depending on the source of the error estimation/indication, these
2534 * values might be signed, so we address that with the following
2535 *
2536 * @code
2537 *   for (double &error_indicator : estimated_error_per_cell)
2538 *   error_indicator = std::abs(error_indicator);
2539 *  
2540 *  
2541 *   GridRefinement::refine_and_coarsen_fixed_number(
2542 *   *this->triangulation, estimated_error_per_cell, 1. / 5., 0.000);
2543 *  
2544 * @endcode
2545 *
2546 * Now get regularity indicators
2547 * For those elements which must be refined, swap to increasing @f$p@f$
2548 * depending on the regularity threshold...
2549 *
2550
2551 *
2552 *
2553 * @code
2554 *   smoothness_indicators =
2555 *   std::vector<double>(this->triangulation->n_active_cells(),
2556 *   std::numeric_limits<double>::max());
2557 *   if (ErrorIndicator::PrimalSolver::min_degree !=
2558 *   ErrorIndicator::PrimalSolver::max_degree)
2559 *   RegularityIndicator::estimate_smoothness(
2560 *   ErrorIndicator::get_eigenfunctions(), 0, smoothness_indicators);
2561 * @endcode
2562 *
2563 * save data
2564 *
2565 * @code
2566 *   this->output_solution();
2567 *   const double threshold_smoothness = smoothness_threshold_fraction;
2568 *   unsigned int num_refined = 0, num_coarsened = 0;
2569 *   if (ErrorIndicator::PrimalSolver::min_degree !=
2570 *   ErrorIndicator::PrimalSolver::max_degree)
2571 *   {
2572 *   for (const auto &cell :
2573 *   ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2574 *   {
2575 *   if (cell->refine_flag_set())
2576 *   ++num_refined;
2577 *   if (cell->coarsen_flag_set())
2578 *   ++num_coarsened;
2579 *   if (cell->refine_flag_set() &&
2580 *   smoothness_indicators[cell->active_cell_index()] >
2581 *   threshold_smoothness &&
2582 *   static_cast<unsigned int>(cell->active_fe_index() + 1) <
2583 *   ErrorIndicator::get_FECollection()->size())
2584 *   {
2585 *   cell->clear_refine_flag();
2586 *   cell->set_active_fe_index(cell->active_fe_index() + 1);
2587 *   }
2588 *   else if (cell->coarsen_flag_set() &&
2589 *   smoothness_indicators[cell->active_cell_index()] <
2590 *   threshold_smoothness &&
2591 *   cell->active_fe_index() != 0)
2592 *   {
2593 *   cell->clear_coarsen_flag();
2594 *  
2595 *   cell->set_active_fe_index(cell->active_fe_index() - 1);
2596 *   }
2597 * @endcode
2598 *
2599 * Here we also impose a limit on how small the cells can become
2600 *
2601 * @code
2602 *   else if (cell->refine_flag_set() && cell->diameter() < 5.0e-6)
2603 *   {
2604 *   cell->clear_refine_flag();
2605 *   if (static_cast<unsigned int>(cell->active_fe_index() + 1) <
2606 *   ErrorIndicator::get_FECollection()->size())
2607 *   cell->set_active_fe_index(cell->active_fe_index() + 1);
2608 *   }
2609 *   }
2610 *   }
2611 *  
2612 * @endcode
2613 *
2614 * Check what the smallest diameter is
2615 *
2616 * @code
2617 *   double min_diameter = std::numeric_limits<double>::max();
2618 *   for (const auto &cell :
2619 *   ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2620 *   if (cell->diameter() < min_diameter)
2621 *   min_diameter = cell->diameter();
2622 *  
2623 *   std::cout << "Min diameter: " << min_diameter << std::endl;
2624 *  
2625 *   ErrorIndicator::synchronize_discretization();
2626 *  
2627 *   (this->triangulation)->execute_coarsening_and_refinement();
2628 *   }
2629 *   }
2630 *   } // namespace Refinement
2631 *  
2632 *   int
2633 *   main(int argc, char **argv)
2634 *   {
2635 *   try
2636 *   {
2637 *   using namespace dealii;
2638 *   using namespace Maxwell;
2639 *   using namespace Refinement;
2640 *   using namespace ErrorIndicators;
2641 *   using namespace RegularityIndicators;
2642 *  
2643 *  
2644 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
2645 *  
2646 *  
2647 *   AssertThrow(
2648 *   Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
2649 *   ExcMessage("This program can only be run in serial, use ./maxwell-hp"));
2650 *  
2651 *   Triangulation<2> triangulation_DWR, triangulation_Kelly;
2652 *   Structures::create_L_waveguide(triangulation_DWR, 2.0);
2653 *   Structures::create_L_waveguide(triangulation_Kelly, 2.0);
2654 *  
2655 *   Refiner<2, KellyErrorIndicator<2>, LegendreIndicator<2>> problem_Kelly(
2656 *   "maxwell-hp.prm",
2657 *   triangulation_Kelly,
2658 *   /*Minimum Degree*/ 2,
2659 *   /*Maximum Degree*/ 5,
2660 *   /*Starting Degree*/ 2);
2661 *  
2662 *   Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
2663 *   problem_DWR("maxwell-hp.prm",
2664 *   triangulation_DWR,
2665 *   /*Minimum Degree*/ 2,
2666 *   /*Maximum Degree*/ 5,
2667 *   /*Starting Degree*/ 2);
2668 *  
2669 * @endcode
2670 *
2671 * The threshold for the hp-decision: too small -> not enough
2672 * @f$h@f$-refinement, too large -> not enough @f$p@f$-refinement
2673 *
2674 * @code
2675 *   double smoothness_threshold = 0.75;
2676 *  
2677 *   std::cout << "Executing refinement for the Kelly strategy!" << std::endl;
2678 *   problem_Kelly.execute_refinement(smoothness_threshold);
2679 *   std::cout << "...Done with Kelly refinement strategy!" << std::endl;
2680 *   std::cout << "Executing refinement for the DWR strategy!" << std::endl;
2681 *   problem_DWR.execute_refinement(smoothness_threshold);
2682 *   std::cout << "...Done with DWR refinement strategy!" << std::endl;
2683 *   }
2684 *  
2685 *   catch (std::exception &exc)
2686 *   {
2687 *   std::cerr << std::endl
2688 *   << std::endl
2689 *   << "----------------------------------------------------"
2690 *   << std::endl;
2691 *   std::cerr << "Exception on processing: " << std::endl
2692 *   << exc.what() << std::endl
2693 *   << "Aborting!" << std::endl
2694 *   << "----------------------------------------------------"
2695 *   << std::endl;
2696 *  
2697 *   return 1;
2698 *   }
2699 *   catch (...)
2700 *   {
2701 *   std::cerr << std::endl
2702 *   << std::endl
2703 *   << "----------------------------------------------------"
2704 *   << std::endl;
2705 *   std::cerr << "Unknown exception!" << std::endl
2706 *   << "Aborting!" << std::endl
2707 *   << "----------------------------------------------------"
2708 *   << std::endl;
2709 *   return 1;
2710 *   }
2711 *  
2712 *   std::cout << std::endl << " Job done." << std::endl;
2713 *  
2714 *   return 0;
2715 *   }
2716 * @endcode
2717
2718
2719*/
make_array_view
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
DataPostprocessor::evaluate_vector_field
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
Definition data_postprocessor.cc:59
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
FEFaceValues::get_present_fe_values
const FEFaceValues< dim, spacedim > & get_present_fe_values() const
FESubfaceValues::get_present_fe_values
const FESubfaceValues< dim, spacedim > & get_present_fe_values() const
FEValuesBase::get_function_gradients
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
Definition fe_values_base.cc:851
FEValues::get_present_fe_values
const FEValues< dim, spacedim > & get_present_fe_values() const
FE_Nedelec
Definition fe_nedelec.h:401
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)
Utilities::MPI::MPI_InitFinalize
Definition mpi.h:1042
Vector::reinit
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
DEAL_II_VERSION_GTE
#define DEAL_II_VERSION_GTE(major, minor, subminor)
Definition config.h:329
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1780
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition dof_handler.h:468
DoFHandler::face_iterator
typename ActiveSelector::face_iterator face_iterator
Definition dof_handler.h:486
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:440
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_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_hessians
@ update_hessians
Second derivatives of shape functions.
Definition fe_update_flags.h:87
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:75
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition fe_update_flags.h:141
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
data
std::vector< index_type > data
Definition mpi.cc:735
size
std::size_t size
Definition mpi.cc:734
GridRefinement::refine_and_coarsen_fixed_number
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
Definition grid_refinement.cc:318
GridTools::diameter
double diameter(const Triangulation< dim, spacedim > &tria)
Definition grid_tools_geometry.cc:43
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition tria_iterator_base.h:37
LocalIntegrators::L2::mass_matrix
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Portable::internal::apply
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
Definition portable_tensor_product_kernels.h:251
SmoothnessEstimator::Legendre::default_fe_series
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
Definition smoothness_estimator.cc:288
SmoothnessEstimator::Legendre::coefficient_decay
void coefficient_decay(FESeries::Legendre< dim, spacedim > &fe_legendre, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const VectorTools::NormType regression_strategy=VectorTools::Linfty_norm, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
Definition smoothness_estimator.cc:101
TensorAccessors::extract
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
Definition tensor_accessors.h:218
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition trilinos_sparse_matrix.cc:64
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm mpi_communicator)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
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
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >()), const bool project_to_boundary_first=false)
hp
Definition hp.h:117
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14871
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
std_cxx17::legendre
double legendre(unsigned int l, double x)
Definition cmath.h:65
std
STL namespace.
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6869
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
swap
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
Definition observer_pointer.h:665
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition p4est_wrappers.cc:68
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Generated by doxygen 1.9.8