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