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