deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
Maxwell-Eigenvalue-hp-Refinement.h
Go to the documentation of this file.
1
243 *  
244 *  
245 *   #include <deal.II/base/function_parser.h>
246 *   #include <deal.II/base/index_set.h>
247 *   #include <deal.II/base/parameter_handler.h>
248 *   #include <deal.II/base/quadrature_lib.h>
249 *   #include <deal.II/base/utilities.h>
250 *  
251 *   #include <deal.II/dofs/dof_handler.h>
252 *   #include <deal.II/dofs/dof_tools.h>
253 *  
254 *   #include <deal.II/fe/fe_nedelec.h>
255 *   #include <deal.II/fe/fe_series.h>
256 *   #include <deal.II/fe/fe_values.h>
257 *  
258 *   #if DEAL_II_VERSION_GTE(9, 7, 0)
259 *   #include <deal.II/grid/cell_data.h>
260 *   #endif
261 *   #include <deal.II/grid/tria.h>
262 *   #include <deal.II/grid/tria_iterator.h>
263 *  
264 *   #include <deal.II/lac/affine_constraints.h>
265 *   #include <deal.II/lac/full_matrix.h>
268 *   #include <deal.II/lac/petsc_vector.h>
269 *   #include <deal.II/lac/slepc_solver.h>
270 *  
271 *   #include <deal.II/numerics/data_out.h>
272 *   #include <deal.II/numerics/vector_tools.h>
273 *  
274 * @endcode
275 *
277 *
278 * @code
279 *   #include <deal.II/base/multithread_info.h>
280 *   #include <deal.II/base/work_stream.h>
281 *  
282 *   #include "petscpc.h"
283 *  
284 * @endcode
285 *
287 *
288 * @code
289 *   #include <deal.II/fe/fe_tools.h>
290 *  
291 *   #include <deal.II/numerics/error_estimator.h>
292 *   #include <deal.II/numerics/smoothness_estimator.h>
293 * @endcode
294 *
296 *
297 * @code
298 *   #include <deal.II/grid/grid_refinement.h>
299 *  
303 *  
304 *   namespace Operations
305 *   {
306 *  
310 *   double
311 *   curlcurl(const ::FEValues<2> &fe_values,
312 *   const unsigned int & i,
313 *   const unsigned int & j,
314 *   const unsigned int & q_point)
315 *   {
316 *   auto gradu1_x1x2 = fe_values.shape_grad_component(i, q_point, 0);
317 *   auto gradu2_x1x2 = fe_values.shape_grad_component(i, q_point, 1);
318 *  
319 *   auto gradv1_x1x2 = fe_values.shape_grad_component(j, q_point, 0);
320 *   auto gradv2_x1x2 = fe_values.shape_grad_component(j, q_point, 1);
321 *   return (gradu2_x1x2[0] - gradu1_x1x2[1]) *
322 *   (gradv2_x1x2[0] - gradv1_x1x2[1]);
323 *   }
324 *  
325 *  
328 *   template <int dim>
329 *   inline double
330 *   dot_term(const ::FEValues<dim> &fe_values,
331 *   const unsigned int & i,
332 *   const unsigned int & j,
333 *   const unsigned int & q_point)
334 *   {
335 *   double output = 0.0;
336 *   for (unsigned int comp = 0; comp < dim; ++comp)
337 *   {
338 *   output += fe_values.shape_value_component(i, q_point, comp) *
339 *   fe_values.shape_value_component(j, q_point, comp);
340 *   }
341 *   return output;
342 *   }
343 *   } // namespace Operations
344 *  
349 *   namespace Structures
350 *   {
351 *   using namespace dealii;
352 *  
353 *   void
354 *   create_L_waveguide(Triangulation<2> &triangulation, const double &scaling)
355 *   {
356 *   const unsigned int dim = 2;
357 *  
358 *   const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
359 *   {scaling * 0.5, scaling * 0.0},
360 *   {scaling * 0.0, scaling * 0.5},
361 *   {scaling * 0.5, scaling * 0.5},
362 *   {scaling * 0.0, scaling * 1.0},
363 *   {scaling * 0.5, scaling * 1.0},
364 *   {scaling * 1.0, scaling * 0.5},
365 *   {scaling * 1.0, scaling * 1.0}};
366 *  
367 *   const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
368 *   cell_vertices = {{{0, 1, 2, 3}}, {{2, 3, 4, 5}}, {{3, 6, 5, 7}}};
369 *   const unsigned int n_cells = cell_vertices.size();
370 *   std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
371 *   for (unsigned int i = 0; i < n_cells; ++i)
372 *   {
373 *   for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
374 *   cells[i].vertices[j] = cell_vertices[i][j];
375 *   cells[i].material_id = 0;
376 *   }
377 *   triangulation.create_triangulation(vertices, cells, SubCellData());
378 *   triangulation.refine_global(1);
379 *   }
380 *  
381 *  
382 *   void
384 *   const double & scaling)
385 *   {
386 *   const unsigned int dim = 2;
387 *  
388 *   const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
389 *   {scaling * 0.6, scaling * 0.0},
390 *   {scaling * 0.0, scaling * 0.3},
391 *   {scaling * 0.6, scaling * 0.3}};
392 *  
393 *   const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
394 *   cell_vertices = {{{0, 1, 2, 3}}};
395 *   const unsigned int n_cells = cell_vertices.size();
396 *   std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
397 *   for (unsigned int i = 0; i < n_cells; ++i)
398 *   {
399 *   for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
400 *   cells[i].vertices[j] = cell_vertices[i][j];
401 *   cells[i].material_id = 0;
402 *   }
403 *   triangulation.create_triangulation(vertices, cells, SubCellData());
404 *   triangulation.refine_global(0);
405 *   }
406 *   } // namespace Structures
407 *  
411 *   namespace Maxwell
412 *   {
413 *   using namespace dealii;
414 *  
415 *   /*
416 *   The "Base" class provides the universal functionality of any eigensolver,
417 *   namely the parameters for the problem, an underlying triangulation, and
418 *   functionality for setting the refinement cycle and to output the solution.
419 *  
420 *   In this case, and for any future class, the use of raw pointers (as opposed to
421 *   "smart" pointers) indicates a lack of ownership. Specifically, the
422 *   triangulation raw pointer is pointing to a triangulation that is owned (and
423 *   created) elsewhere.
424 *   */
425 *   template <int dim>
426 *   class Base
427 *   {
428 *   public:
429 *   Base(const std::string &prm_file, Triangulation<dim> &coarse_grid);
430 *  
431 *   virtual unsigned int
432 *   solve_problem() = 0; // Implemented by a derived class
433 *   virtual void
434 *   set_refinement_cycle(const unsigned int cycle);
435 *  
436 *   virtual void
437 *   output_solution() = 0; // Implemented by a derived class
438 *  
439 *  
440 *   protected:
442 *   unsigned int refinement_cycle = 0;
443 *   std::unique_ptr<ParameterHandler> parameters;
444 *   unsigned int n_eigenpairs = 1;
445 *   double target = 0.0;
446 *   unsigned int eigenpair_selection_scheme;
447 *   unsigned int max_cycles = 0;
448 *   ompi_communicator_t * mpi_communicator = PETSC_COMM_SELF;
449 *   };
450 *  
451 *  
454 *   template <int dim>
455 *   Base<dim>::Base(const std::string &prm_file, Triangulation<dim> &coarse_grid)
456 *   : triangulation(&coarse_grid)
457 *   , parameters(std::make_unique<ParameterHandler>())
458 *   {
459 *   parameters->declare_entry(
460 *   "Eigenpair selection scheme",
461 *   "1",
462 *   Patterns::Integer(0, 1),
463 *   "The type of eigenpairs to find (0 - smallest, 1 - target)");
464 *   parameters->declare_entry("Number of eigenvalues/eigenfunctions",
465 *   "1",
466 *   Patterns::Integer(0, 100),
467 *   "The number of eigenvalues/eigenfunctions "
468 *   "to be computed.");
469 *   parameters->declare_entry("Target eigenvalue",
470 *   "1",
472 *   "The target eigenvalue (if scheme == 1)");
473 *  
474 *   parameters->declare_entry("Cycles number",
475 *   "1",
476 *   Patterns::Integer(0, 1500),
477 *   "The number of cycles in refinement");
478 *   parameters->parse_input(prm_file);
479 *  
481 *   parameters->get_integer("Eigenpair selection scheme");
482 *  
483 * @endcode
484 *
485 * The project currently only supports selection by a target eigenvalue.
487 *
488 * @code
490 *   "Selection by a target is the only currently supported option!");
491 *   n_eigenpairs =
492 *   parameters->get_integer("Number of eigenvalues/eigenfunctions");
493 *   assert(
494 *   n_eigenpairs == 1 &&
495 *   "Only the computation of a single eigenpair is currently supported!");
496 *  
497 *   target = parameters->get_double("Target eigenvalue");
498 *   max_cycles = parameters->get_integer("Cycles number");
500 *   n_eigenpairs = 1;
501 *   }
502 *  
503 *   template <int dim>
504 *   void
505 *   Base<dim>::set_refinement_cycle(const unsigned int cycle)
506 *   {
507 *   refinement_cycle = cycle;
508 *   }
509 *  
510 *  
515 *   template <int dim>
516 *   class EigenSolver : public virtual Base<dim>
517 *   {
518 *   public:
519 *   EigenSolver(const std::string & prm_file,
520 *   Triangulation<dim> &coarse_grid,
521 *   const unsigned int &minimum_degree,
522 *   const unsigned int &maximum_degree,
523 *   const unsigned int &starting_degree);
524 *  
525 *   virtual unsigned int
526 *   solve_problem() override;
527 *  
528 *   virtual unsigned int
529 *   n_dofs() const;
530 *  
531 *   template <class SolverType>
532 *   void
534 *  
535 *   virtual void
536 *   setup_system();
537 *  
538 *   virtual void
540 *  
541 *   protected:
542 *   const std::unique_ptr<hp::FECollection<dim>> fe_collection;
543 *   std::unique_ptr<hp::QCollection<dim>> quadrature_collection;
544 *   std::unique_ptr<hp::QCollection<dim - 1>> face_quadrature_collection;
545 *   DoFHandler<dim> dof_handler;
546 *   const unsigned int max_degree, min_degree;
547 * @endcode
548 *
549 * for the actual solution
550 *
551 * @code
552 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> eigenfunctions;
553 *   std::unique_ptr<std::vector<double>> eigenvalues;
554 *   Vector<double> solution;
555 *  
556 *   double *
557 *   get_lambda_h();
558 *  
560 *   get_solution();
561 *  
562 *   void
564 *  
565 *   private:
566 *   AffineConstraints<double> constraints;
568 *   };
569 *  
570 *  
574 *   template <int dim>
575 *   EigenSolver<dim>::EigenSolver(const std::string & prm_file,
577 *   const unsigned int &minimum_degree,
578 *   const unsigned int &maximum_degree,
579 *   const unsigned int &starting_degree)
581 *   , fe_collection(std::make_unique<hp::FECollection<dim>>())
582 *   , quadrature_collection(std::make_unique<hp::QCollection<dim>>())
583 *   , face_quadrature_collection(std::make_unique<hp::QCollection<dim - 1>>())
584 *   , dof_handler(triangulation)
585 *   , max_degree(maximum_degree)
586 *   , min_degree(minimum_degree)
587 *   , eigenfunctions(
588 *   std::make_unique<std::vector<PETScWrappers::MPI::Vector>>())
589 *   , eigenvalues(std::make_unique<std::vector<double>>())
590 *   {
591 *   for (unsigned int degree = min_degree; degree <= max_degree; ++degree)
592 *   {
593 *   fe_collection->push_back(FE_Nedelec<dim>(degree - 1));
594 * @endcode
595 *
596 * Generate quadrature collection with sorted quadrature weights
597 *
598 * @code
599 *   const QGauss<dim> quadrature(degree + 1);
600 *   const QSorted<dim> sorted_quadrature(quadrature);
602 *  
603 *   const QGauss<dim - 1> face_quadrature(degree + 1);
604 *   const QSorted<dim - 1> sorted_face_quadrature(face_quadrature);
605 *   face_quadrature_collection->push_back(sorted_face_quadrature);
606 *   }
607 * @endcode
608 *
610 *
611 * @code
612 *   if (starting_degree > min_degree && starting_degree <= max_degree)
613 *   {
614 *   const unsigned int start_diff = starting_degree - min_degree;
616 *   cell1 = dof_handler.begin_active(),
617 *   endc1 = dof_handler.end();
618 *   for (; cell1 < endc1; ++cell1)
619 *   {
620 *   cell1->set_active_fe_index(start_diff);
621 *   }
622 *   }
623 *   }
624 *  
625 *  
629 *   template <int dim>
630 *   double *
632 *   {
633 *   return &(*eigenvalues)[0];
634 *   }
635 *  
636 *  
640 *   template <int dim>
643 *   {
644 *   return &solution;
645 *   }
646 *  
647 *  
650 *   template <int dim>
651 *   void
653 *   {
654 *   solution.reinit((*eigenfunctions)[0].size());
655 *   for (unsigned int i = 0; i < solution.size(); ++i)
656 *   solution[i] = (*eigenfunctions)[0][i];
657 *   }
658 *  
659 *  
665 *   template <int dim>
666 *   template <class SolverType>
667 *   void
669 *   {
670 * @endcode
671 *
672 * From the parameters class, initialize the eigensolver...
673 *
674 * @code
675 *   switch (this->eigenpair_selection_scheme)
676 *   {
677 *   case 1:
678 *   eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
679 * @endcode
680 *
681 * eigensolver.set_target_eigenvalue(this->target);
682 *
683 * @code
684 *   break;
685 *   default:
686 *   eigensolver.set_which_eigenpairs(EPS_SMALLEST_MAGNITUDE);
687 *  
688 *   break;
689 *   }
690 *   eigensolver.set_problem_type(EPS_GHEP);
691 * @endcode
692 *
694 *
695
696 *
697 *
698 * @code
699 *   double shift_scalar = this->parameters->get_double("Target eigenvalue");
700 * @endcode
701 *
702 * //For the shift-and-invert transformation
703 *
704 * @code
706 *   shift_scalar);
708 *   this->mpi_communicator, additional_data);
709 *  
710 *   eigensolver.set_transformation(spectral_transformation);
711 *   eigensolver.set_target_eigenvalue(this->target);
712 *   }
713 *  
714 *  
718 *   template <int dim>
719 *   unsigned int
721 *   {
722 *   setup_system();
724 *  
725 *   SolverControl solver_control(dof_handler.n_dofs() * 10,
726 *   5.0e-8,
727 *   false,
728 *   false);
730 *   this->mpi_communicator);
731 *  
733 *  
734 * @endcode
735 *
736 * solve the problem
737 *
738 * @code
740 *   mass_matrix,
741 *   *eigenvalues,
743 *   eigenfunctions->size());
744 *   for (auto &entry : *eigenfunctions)
745 *   {
746 *   constraints.distribute(entry);
747 *   }
749 *  
750 *   return solver_control.last_step();
751 *   }
752 *  
753 *   template <int dim>
754 *   unsigned int
756 *   {
757 *   return dof_handler.n_dofs();
758 *   }
759 *  
760 *  
765 *   template <int dim>
766 *   void
768 *   {
769 *   dof_handler.distribute_dofs(*fe_collection);
770 *   constraints.clear();
771 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
772 *   DoFTools::make_zero_boundary_constraints(dof_handler, constraints);
773 *   constraints.close();
774 *  
775 *   eigenfunctions->resize(this->n_eigenpairs);
776 *   eigenvalues->resize(this->n_eigenpairs);
777 *  
778 *   IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs();
779 *  
780 *   for (auto &entry : *eigenfunctions)
781 *   {
783 *   }
784 *   }
785 *  
786 *  
789 *   template <int dim>
790 *   void
792 *   {
793 *   hp::FEValues<dim> hp_fe_values(*fe_collection,
798 * @endcode
799 *
800 * Prep the system matrices for the solution
801 *
802 * @code
803 *   stiffness_matrix.reinit(dof_handler.n_dofs(),
804 *   dof_handler.n_dofs(),
805 *   dof_handler.max_couplings_between_dofs());
806 *   mass_matrix.reinit(dof_handler.n_dofs(),
807 *   dof_handler.n_dofs(),
808 *   dof_handler.max_couplings_between_dofs());
809 *  
811 *   std::vector<types::global_dof_index> local_dof_indices;
812 *  
813 *   for (const auto &cell : dof_handler.active_cell_iterators())
814 *   {
815 *   const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
816 *  
817 *   cell_stiffness_matrix.reinit(dofs_per_cell, dofs_per_cell);
819 *  
820 *   cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
821 *   cell_mass_matrix = 0;
822 *  
823 *   hp_fe_values.reinit(cell);
824 *  
825 *   const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
826 *  
827 *   for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
828 *   ++q_point)
829 *   {
830 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
831 *   {
832 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
833 *   {
834 * @endcode
835 *
836 * Note that (in general) the Nedelec element is not
837 * primitive, namely that the shape functions are vectorial
838 * with components in more than one direction
839 *
840
841 *
842 *
843 * @code
844 *   cell_stiffness_matrix(i, j) +=
845 *   Operations::curlcurl(fe_values, i, j, q_point) *
846 *   fe_values.JxW(q_point);
847 *  
848 *   cell_mass_matrix(i, j) +=
849 *   (Operations::dot_term(fe_values, i, j, q_point)) *
850 *   fe_values.JxW(q_point);
851 *   }
852 *   }
853 *   local_dof_indices.resize(dofs_per_cell);
854 *   cell->get_dof_indices(local_dof_indices);
855 *   }
856 *  
857 *   constraints.distribute_local_to_global(cell_stiffness_matrix,
858 *   local_dof_indices,
860 *   constraints.distribute_local_to_global(cell_mass_matrix,
861 *   local_dof_indices,
862 *   mass_matrix);
863 *   }
866 *  
867 *   for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
868 *   if (constraints.is_constrained(i))
869 *   {
870 *   stiffness_matrix.set(i, i, 10000.0);
871 *   mass_matrix.set(i, i, 1);
872 *   }
873 * @endcode
874 *
875 * since we have just set individual elements, we need the following
876 *
877 * @code
880 *   }
881 *  
882 *  
886 *   template <int dim>
887 *   class PrimalSolver : public EigenSolver<dim>
888 *   {
889 *   public:
890 *   PrimalSolver(const std::string & prm_file,
892 *   const unsigned int &min_degree,
893 *   const unsigned int &max_degree,
894 *   const unsigned int &starting_degree);
895 *  
896 *   virtual void
898 *   override; // Implements the output solution of the base class...
899 *   virtual unsigned int
900 *   n_dofs() const override;
901 *   };
902 *  
903 *   template <int dim>
904 *   PrimalSolver<dim>::PrimalSolver(const std::string & prm_file,
906 *   const unsigned int &min_degree,
907 *   const unsigned int &max_degree,
908 *   const unsigned int &starting_degree)
912 *   min_degree,
913 *   max_degree,
915 *   {}
916 *  
917 *  
921 *   template <int dim>
922 *   void
924 *   {
925 *   DataOut<dim> data_out;
926 *   data_out.attach_dof_handler(this->dof_handler);
927 *   Vector<double> fe_degrees(this->triangulation->n_active_cells());
928 *   for (const auto &cell : this->dof_handler.active_cell_iterators())
929 *   fe_degrees(cell->active_cell_index()) =
930 *   (*this->fe_collection)[cell->active_fe_index()].degree;
931 *   data_out.add_data_vector(fe_degrees, "fe_degree");
932 *   data_out.add_data_vector((*this->eigenfunctions)[0],
933 *   std::string("eigenfunction_no_") +
935 *  
936 *   std::cout << "Eigenvalue: " << (*this->eigenvalues)[0]
937 *   << " NDoFs: " << this->dof_handler.n_dofs() << std::endl;
938 *   std::ofstream eigenvalues_out(
939 *   "eigenvalues-" + std::to_string(this->refinement_cycle) + ".txt");
940 *  
941 *   eigenvalues_out << std::setprecision(20) << (*this->eigenvalues)[0] << " "
942 *   << this->dof_handler.n_dofs() << std::endl;
943 *  
944 *   eigenvalues_out.close();
945 *  
946 *  
947 *   data_out.build_patches();
948 *   std::ofstream output("eigenvectors-" +
949 *   std::to_string(this->refinement_cycle) + ".vtu");
950 *   data_out.write_vtu(output);
951 *   }
952 *  
953 *   template <int dim>
954 *   unsigned int
956 *   {
957 *   return EigenSolver<dim>::n_dofs();
958 *   }
959 *  
960 * @endcode
961 *
964 * however, it is convenient to separate them in this manner (e.g., for
966 *
967 * @code
968 *   template <int dim>
969 *   class DualSolver : public EigenSolver<dim>
970 *   {
971 *   public:
972 *   DualSolver(const std::string & prm_file,
974 *   const unsigned int &min_degree,
975 *   const unsigned int &max_degree,
976 *   const unsigned int &starting_degree);
977 *   };
978 *  
979 *   template <int dim>
980 *   DualSolver<dim>::DualSolver(const std::string & prm_file,
982 *   const unsigned int &min_degree,
983 *   const unsigned int &max_degree,
984 *   const unsigned int &starting_degree)
988 *   min_degree,
989 *   max_degree,
991 *   {}
992 *  
993 *   } // namespace Maxwell
994 *  
998 *   namespace ErrorIndicators
999 *   {
1000 *   using namespace Maxwell;
1001 *  
1002 *  
1007 *   template <int dim, bool report_dual>
1008 *   class DualWeightedResidual : public PrimalSolver<dim>, public DualSolver<dim>
1009 *   {
1010 *   public:
1011 *   void
1012 *   output_eigenvalue_data(std::ofstream &os);
1013 *   void
1014 *   output_qoi_error_estimates(std::ofstream &os);
1015 *  
1016 *   std::string
1017 *   name() const
1018 *   {
1019 *   return "DWR";
1020 *   }
1021 *   DualWeightedResidual(const std::string & prm_file,
1023 *   const unsigned int &min_primal_degree,
1024 *   const unsigned int &max_primal_degree,
1025 *   const unsigned int &starting_primal_degree);
1026 *  
1027 *   virtual unsigned int
1028 *   solve_problem() override;
1029 *  
1030 *   virtual void
1031 *   output_solution() override;
1032 *  
1033 *   virtual unsigned int
1034 *   n_dofs() const override;
1035 *  
1036 *   void
1038 *  
1040 *   get_DoFHandler();
1041 *  
1044 *  
1047 *  
1049 *   get_FECollection();
1050 *  
1053 *  
1054 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1056 *  
1057 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1059 *  
1060 *   std::unique_ptr<std::vector<double>> &
1062 *  
1063 *   std::unique_ptr<std::vector<double>> &
1065 *  
1066 *   void
1068 *  
1069 *   unsigned int
1070 *   get_max_degree()
1071 *   {
1072 *   return PrimalSolver<dim>::fe_collection->max_degree();
1073 *   }
1074 *   double qoi_error_estimate = 0;
1075 *  
1076 *   private:
1077 *   void
1078 *   embed(const DoFHandler<dim> & dof1,
1079 *   const DoFHandler<dim> & dof2,
1080 *   const AffineConstraints<double> &constraints,
1081 *   const Vector<double> & solution,
1082 *   Vector<double> & u2);
1083 *  
1084 *   void
1085 *   extract(const DoFHandler<dim> & dof1,
1086 *   const DoFHandler<dim> & dof2,
1087 *   const AffineConstraints<double> &constraints,
1088 *   const Vector<double> & solution,
1089 *   Vector<double> & u2);
1090 *  
1091 *  
1092 *  
1093 *   /*The following FEValues objects are unique_ptrs to 1) avoid default
1094 *   constructors for these objects, and 2) automate memory management*/
1095 *   std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values;
1096 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values;
1097 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor;
1098 *   std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values;
1099 *  
1100 *   std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values_forward;
1101 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_forward;
1102 *   std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor_forward;
1103 *   std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values_forward;
1104 *   using FaceIntegrals =
1105 *   typename std::map<typename DoFHandler<dim>::face_iterator, double>;
1106 *  
1107 *   unsigned int
1108 *   solve_primal_problem();
1109 *  
1110 *   unsigned int
1111 *   solve_dual_problem();
1112 *  
1113 *   void
1116 *  
1117 *   double
1120 *  
1121 *   void
1123 *  
1124 *   void
1126 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1128 *   const Vector<double> & dual_weights,
1129 *   const double & lambda_h,
1132 *  
1133 *   void
1135 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1137 *   const Vector<double> & dual_weights,
1138 *   const double & lambda_h,
1140 *  
1141 *   void
1143 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1144 *   const unsigned int & face_no,
1146 *   const Vector<double> & dual_weights,
1148 *  
1149 *   void
1151 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1152 *   const unsigned int & face_no,
1154 *   const Vector<double> & dual_weights,
1156 *   };
1157 *  
1158 *  
1162 *   template <int dim, bool report_dual>
1164 *   const std::string & prm_file,
1166 *   const unsigned int &min_primal_degree,
1167 *   const unsigned int &max_primal_degree,
1168 *   const unsigned int &starting_primal_degree)
1171 *   triangulation,
1176 *   triangulation,
1177 *   min_primal_degree + 1,
1178 *   max_primal_degree + 1,
1180 *   {
1182 *   }
1183 *  
1184 *  
1188 *   template <int dim, bool report_dual>
1191 *   {
1192 *   if (!report_dual)
1194 *   else
1195 *   return &(DualSolver<dim>::dof_handler);
1196 *   }
1197 *  
1198 * @endcode
1199 *
1200 * See above function, but to specifically output the primal DoFHandler...
1201 *
1202 * @code
1203 *   template <int dim, bool report_dual>
1206 *   {
1208 *   }
1209 *  
1210 * @endcode
1211 *
1212 * See above function, but for the FECollection
1213 *
1214 * @code
1215 *   template <int dim, bool report_dual>
1218 *   {
1219 *   if (!report_dual)
1221 *   else
1223 *   }
1224 *  
1225 * @endcode
1226 *
1227 * See above function, but for the primal FECollection
1228 *
1229 * @code
1230 *   template <int dim, bool report_dual>
1233 *   {
1235 *   }
1236 *  
1237 *   template <int dim, bool report_dual>
1240 *   {
1241 *   return &(DualSolver<dim>::dof_handler);
1242 *   }
1243 *  
1244 *  
1245 *   template <int dim, bool report_dual>
1246 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1248 *   {
1249 *   if (!report_dual)
1251 *   else
1253 *   }
1254 *  
1255 *  
1256 *   template <int dim, bool report_dual>
1257 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1259 *   {
1261 *   }
1262 *  
1263 *  
1264 *   template <int dim, bool report_dual>
1265 *   std::unique_ptr<std::vector<double>> &
1267 *   {
1269 *   }
1270 *  
1271 *  
1272 *   template <int dim, bool report_dual>
1273 *   std::unique_ptr<std::vector<double>> &
1275 *   {
1277 *   }
1278 *  
1279 *   template <int dim, bool report_dual>
1280 *   void
1282 *   {
1284 *   }
1285 *  
1286 * @endcode
1287 *
1289 *
1290 * @code
1291 *   template <int dim, bool report_dual>
1292 *   unsigned int
1294 *   {
1296 *   }
1297 *  
1298 * @endcode
1299 *
1301 *
1302 * @code
1303 *   template <int dim, bool report_dual>
1304 *   unsigned int
1306 *   {
1308 *   }
1309 *  
1310 *  
1314 *   template <int dim, bool report_dual>
1315 *   unsigned int
1317 *   {
1320 *   }
1321 *  
1322 *  
1325 *   template <int dim, bool report_dual>
1326 *   unsigned int
1328 *   {
1329 *   return PrimalSolver<dim>::n_dofs();
1330 *   }
1331 *  
1332 *  
1338 *   template <int dim, bool report_dual>
1339 *   void
1341 *   {
1342 *   /*Note: No additional checks need to be made ensuring that these operations
1343 *   are legal as these checks are made prior to entering this function (i.e.,
1344 *   if the primal attains a degree N,
1345 *   then, by construction, a degree of N+1 must be permissible for the
1346 *   dual)*/
1349 *  
1350 *   if (report_dual)
1351 *   {
1352 * @endcode
1353 *
1354 * In this case, we have modified the polynomial orders for the dual;
1356 *
1357 * @code
1360 *   }
1361 *   typename DoFHandler<dim>::active_cell_iterator cell1 = dof1->begin_active(),
1362 *   endc1 = dof1->end();
1363 *   typename DoFHandler<dim>::active_cell_iterator cell2 = dof2->begin_active();
1364 *   for (; cell1 < endc1; ++cell1, ++cell2)
1365 *   {
1366 *   cell2->set_active_fe_index(cell1->active_fe_index());
1367 *   }
1368 *   }
1369 *  
1370 *  
1374 *   template <int dim, bool report_dual>
1375 *   void
1377 *   {
1378 * @endcode
1379 *
1380 * initialize the cell fe_values...
1381 *
1382 * @code
1383 *   cell_hp_fe_values = std::make_unique<hp::FEValues<dim>>(
1388 *   face_hp_fe_values = std::make_unique<hp::FEFaceValues<dim>>(
1393 *   face_hp_fe_values_neighbor = std::make_unique<hp::FEFaceValues<dim>>(
1398 *   subface_hp_fe_values = std::make_unique<hp::FESubfaceValues<dim>>(
1402 *   }
1403 *  
1404 *  
1409 *   template <int dim, bool report_dual>
1410 *   void
1414 *   {
1415 *   double sum_primal = 0.0, sum_dual = 0.0;
1416 *   for (const auto &cell :
1417 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1418 *   {
1419 *   cell_hp_fe_values->reinit(cell);
1420 *  
1421 * @endcode
1422 *
1423 * grab the fe_values object
1424 *
1425 * @code
1426 *   const FEValues<dim> &fe_values =
1427 *   cell_hp_fe_values->get_present_fe_values();
1428 *  
1429 *   std::vector<Vector<double>> cell_primal_values(
1430 *   fe_values.n_quadrature_points, Vector<double>(dim)),
1431 *   cell_dual_values(fe_values.n_quadrature_points, Vector<double>(dim));
1432 *   fe_values.get_function_values(primal_solution, cell_primal_values);
1433 *   fe_values.get_function_values(dual_weights, cell_dual_values);
1434 *  
1435 *  
1436 *   for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1437 *   {
1438 *   sum_primal +=
1439 *   cell_primal_values[p] * cell_primal_values[p] * fe_values.JxW(p);
1440 *   sum_dual +=
1441 *   cell_dual_values[p] * cell_dual_values[p] * fe_values.JxW(p);
1442 *   }
1443 *   }
1444 *  
1447 *   }
1448 *  
1449 *  
1453 *   template <int dim, bool report_dual>
1454 *   void
1457 *   {
1458 * @endcode
1459 *
1460 * The constraints could be grabbed directly, but this is simple
1461 *
1462 * @code
1467 *  
1472 *  
1473 * @endcode
1474 *
1475 * First map the primal solution to the space of the dual solution
1476 * This allows us to use just one set of FEValues objects (rather than one
1477 * set for the primal, one for dual)
1478 *
1479
1480 *
1481 *
1482 * @code
1484 *  
1490 *  
1492 *  
1494 *  
1497 *  
1498 * @endcode
1499 *
1501 *
1502 * @code
1508 *  
1509 * @endcode
1510 *
1511 * Now embed this back to the space of the dual solution
1512 *
1513 * @code
1518 *   dual_weights);
1519 *  
1520 *  
1521 * @endcode
1522 *
1523 * Subtract this from the full dual solution
1524 *
1525 * @code
1527 *   dual_weights *= -1.0;
1528 *  
1530 *  
1532 *   for (const auto &cell :
1533 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1534 *   for (const auto &face : cell->face_iterators())
1535 *   face_integrals[face] = -1e20;
1536 *  
1537 *  
1538 *   for (const auto &cell :
1539 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1540 *   {
1541 *   estimate_on_one_cell(cell,
1543 *   dual_weights,
1544 *   *(PrimalSolver<dim>::get_lambda_h()),
1546 *   face_integrals);
1547 *   }
1548 *   unsigned int present_cell = 0;
1549 *   for (const auto &cell :
1550 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1551 *   {
1552 *   for (const auto &face : cell->face_iterators())
1553 *   {
1554 *   Assert(face_integrals.find(face) != face_integrals.end(),
1555 *   ExcInternalError());
1556 *   error_indicators(present_cell) -= 0.5 * face_integrals[face];
1557 *   }
1558 *   ++present_cell;
1559 *   }
1560 *  
1561 * @endcode
1562 *
1564 * estimate of the QoI error
1565 *
1566 * @code
1567 *   this->qoi_error_estimate =
1570 *   std::cout << "Estimated QoI error: " << std::setprecision(20)
1571 *   << qoi_error_estimate << std::endl;
1572 *   }
1573 *  
1574 *  
1575 *  
1578 *   template <int dim, bool report_dual>
1579 *   void
1581 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1583 *   const Vector<double> & dual_weights,
1584 *   const double & lambda_h,
1587 *   {
1590 *   for (unsigned int face_no : GeometryInfo<dim>::face_indices())
1591 *   {
1592 *   if (cell->face(face_no)->at_boundary())
1593 *   {
1594 *   face_integrals[cell->face(face_no)] = 0.0;
1595 *   continue;
1596 *   }
1597 *   if ((cell->neighbor(face_no)->has_children() == false) &&
1598 *   (cell->neighbor(face_no)->level() == cell->level()) &&
1599 *   (cell->neighbor(face_no)->index() < cell->index()))
1600 *   continue;
1601 *   if (cell->at_boundary(face_no) == false)
1602 *   if (cell->neighbor(face_no)->level() < cell->level())
1603 *   continue;
1604 *   if (cell->face(face_no)->has_children() == false)
1607 *   else
1610 *   }
1611 *   }
1612 *  
1613 *  
1616 *   template <int dim, bool report_dual>
1617 *   void
1619 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1621 *   const Vector<double> & dual_weights,
1622 *   const double & lambda_h,
1624 *   {
1625 *   cell_hp_fe_values->reinit(cell);
1626 * @endcode
1627 *
1628 * Grab the fe_values object
1629 *
1630 * @code
1631 *   const FEValues<dim> &fe_values = cell_hp_fe_values->get_present_fe_values();
1632 *   std::vector<std::vector<Tensor<2, dim, double>>> cell_hessians(
1633 *   fe_values.n_quadrature_points, std::vector<Tensor<2, dim, double>>(dim));
1634 *   std::vector<Vector<double>> cell_primal_values(
1635 *   fe_values.n_quadrature_points, Vector<double>(dim)),
1636 *   cell_dual_values(fe_values.n_quadrature_points, Vector<double>(dim));
1637 *   fe_values.get_function_values(primal_solution, cell_primal_values);
1638 *   fe_values.get_function_hessians(primal_solution, cell_hessians);
1639 *   fe_values.get_function_values(dual_weights, cell_dual_values);
1640 *  
1641 *  
1642 *  
1643 *   double sum = 0.0;
1644 *   for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1645 *   {
1646 *   sum +=
1647 *   (/*x-component*/ (cell_hessians[p][1][1][0] -
1648 *   cell_hessians[p][0][1][1]) *
1649 *   (cell_dual_values[p](0)) +
1650 *   /*y-component*/
1651 *   (cell_hessians[p][0][0][1] - cell_hessians[p][1][0][0]) *
1652 *   (cell_dual_values[p](1)) -
1653 *   lambda_h * (cell_primal_values[p](0) * cell_dual_values[p](0) +
1654 *   cell_primal_values[p](1) * cell_dual_values[p](1))) *
1655 *   fe_values.JxW(p);
1656 *   }
1657 *  
1658 *   error_indicators(cell->active_cell_index()) += sum;
1659 *   }
1660 *  
1661 *  
1664 *   template <int dim, bool report_dual>
1665 *   void
1667 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1668 *   const unsigned int & face_no,
1670 *   const Vector<double> & dual_weights,
1672 *   {
1673 *   Assert(cell->neighbor(face_no).state() == IteratorState::valid,
1674 *   ExcInternalError());
1675 *   const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1676 *   const auto neighbor = cell->neighbor(face_no);
1677 *  
1678 *   const unsigned int quadrature_index =
1679 *   std::max(cell->active_fe_index(), neighbor->active_fe_index());
1680 *   face_hp_fe_values->reinit(cell, face_no, quadrature_index);
1682 *   face_hp_fe_values->get_present_fe_values();
1683 *   std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1684 *   fe_face_values_cell.n_quadrature_points,
1685 *   std::vector<Tensor<1, dim, double>>(dim)),
1686 *   neighbor_primal_grads(fe_face_values_cell.n_quadrature_points,
1687 *   std::vector<Tensor<1, dim, double>>(dim));
1688 *   fe_face_values_cell.get_function_gradients(primal_solution,
1690 *  
1691 *   face_hp_fe_values_neighbor->reinit(neighbor,
1695 *   face_hp_fe_values_neighbor->get_present_fe_values();
1696 *   fe_face_values_cell_neighbor.get_function_gradients(primal_solution,
1698 *   const unsigned int n_q_points = fe_face_values_cell.n_quadrature_points;
1699 *   double face_integral = 0.0;
1700 *   std::vector<Vector<double>> cell_dual_values(n_q_points,
1701 *   Vector<double>(dim));
1702 *   fe_face_values_cell.get_function_values(dual_weights, cell_dual_values);
1703 *   for (unsigned int p = 0; p < n_q_points; ++p)
1704 *   {
1705 *   auto face_normal = fe_face_values_cell.normal_vector(p);
1706 *  
1707 *   face_integral +=
1708 *   (cell_primal_grads[p][1][0] - cell_primal_grads[p][0][1] -
1709 *   neighbor_primal_grads[p][1][0] + neighbor_primal_grads[p][0][1]) *
1710 *   (cell_dual_values[p][0] * face_normal[1] -
1711 *   cell_dual_values[p][1] * face_normal[0]) *
1712 *   fe_face_values_cell.JxW(p);
1713 *   }
1714 *   Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
1715 *   ExcInternalError());
1716 *   Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
1717 *   face_integrals[cell->face(face_no)] = face_integral;
1718 *   }
1719 *  
1720 *  
1723 *   template <int dim, bool report_dual>
1724 *   void
1726 *   const typename DoFHandler<dim>::active_cell_iterator &cell,
1727 *   const unsigned int & face_no,
1729 *   const Vector<double> & dual_weights,
1731 *   {
1732 *   const typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
1733 *   const typename DoFHandler<dim>::cell_iterator neighbor =
1734 *   cell->neighbor(face_no);
1735 *  
1736 *   Assert(neighbor.state() == IteratorState::valid, ExcInternalError());
1737 *   Assert(neighbor->has_children(), ExcInternalError());
1738 *   (void)neighbor;
1739 *   const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1740 *   for (unsigned int subface_no = 0; subface_no < face->n_children();
1741 *   ++subface_no)
1742 *   {
1744 *   cell->neighbor_child_on_subface(face_no, subface_no);
1746 *   cell->face(face_no)->child(subface_no),
1747 *   ExcInternalError());
1748 *   const unsigned int quadrature_index =
1749 *   std::max(cell->active_fe_index(), neighbor_child->active_fe_index());
1750 * @endcode
1751 *
1752 * initialize fe_subface values_cell
1753 *
1754 * @code
1755 *   subface_hp_fe_values->reinit(cell,
1756 *   face_no,
1757 *   subface_no,
1760 *   subface_hp_fe_values->get_present_fe_values();
1761 *   std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1762 *   subface_fe_values_cell.n_quadrature_points,
1763 *   std::vector<Tensor<1, dim, double>>(dim)),
1764 *   neighbor_primal_grads(subface_fe_values_cell.n_quadrature_points,
1765 *   std::vector<Tensor<1, dim, double>>(dim));
1766 *   subface_fe_values_cell.get_function_gradients(primal_solution,
1768 * @endcode
1769 *
1770 * initialize fe_face_values_neighbor
1771 *
1772 * @code
1777 *   face_hp_fe_values_neighbor->get_present_fe_values();
1778 *   face_fe_values_neighbor.get_function_gradients(primal_solution,
1780 *   const unsigned int n_q_points =
1781 *   subface_fe_values_cell.n_quadrature_points;
1782 *   std::vector<Vector<double>> cell_dual_values(n_q_points,
1783 *   Vector<double>(dim));
1784 *   face_fe_values_neighbor.get_function_values(dual_weights,
1786 *  
1787 *   double face_integral = 0.0;
1788 *  
1789 *   for (unsigned int p = 0; p < n_q_points; ++p)
1790 *   {
1791 *   auto face_normal = face_fe_values_neighbor.normal_vector(p);
1792 *   face_integral +=
1793 *   (cell_primal_grads[p][0][1] - cell_primal_grads[p][1][0] +
1794 *   neighbor_primal_grads[p][1][0] -
1795 *   neighbor_primal_grads[p][0][1]) *
1796 *   (cell_dual_values[p][0] * face_normal[1] -
1797 *   cell_dual_values[p][1] * face_normal[0]) *
1798 *   face_fe_values_neighbor.JxW(p);
1799 *   }
1801 *   }
1802 *   double sum = 0.0;
1803 *   for (unsigned int subface_no = 0; subface_no < face->n_children();
1804 *   ++subface_no)
1805 *   {
1806 *   Assert(face_integrals.find(face->child(subface_no)) !=
1807 *   face_integrals.end(),
1808 *   ExcInternalError());
1809 *   Assert(face_integrals[face->child(subface_no)] != -1e20,
1810 *   ExcInternalError());
1811 *   sum += face_integrals[face->child(subface_no)];
1812 *   }
1813 *   face_integrals[face] = sum;
1814 *   }
1815 *  
1816 *   template <int dim, bool report_dual>
1817 *   double
1821 *   {
1822 *   auto dual_less_primal =
1823 *   dual_solution; // Note: We have already extracted the primal solution...
1824 *  
1825 *  
1826 *   double scaling_factor = 0.0;
1827 *   for (const auto &cell :
1828 *   DualSolver<dim>::dof_handler.active_cell_iterators())
1829 *   {
1830 *   cell_hp_fe_values->reinit(cell);
1831 * @endcode
1832 *
1833 * grab the fe_values object
1834 *
1835 * @code
1836 *   const FEValues<dim> &fe_values =
1837 *   cell_hp_fe_values->get_present_fe_values();
1838 *  
1839 *   std::vector<Vector<double>> cell_values(fe_values.n_quadrature_points,
1840 *   Vector<double>(dim));
1841 *   fe_values.get_function_values(dual_less_primal, cell_values);
1842 *  
1843 *   for (unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1844 *   {
1845 *   scaling_factor +=
1846 *   (cell_values[p] * cell_values[p]) * fe_values.JxW(p);
1847 *   }
1848 *   }
1849 *   double global_QoI_error = 0.0;
1850 *   for (const auto &indicator : error_indicators)
1851 *   {
1853 *   }
1854 *  
1855 *   global_QoI_error /= (1 - 0.5 * scaling_factor);
1856 *   return global_QoI_error;
1857 *   }
1858 *  
1859 *  
1860 *   template <int dim, bool report_dual>
1861 *   void
1863 *   const DoFHandler<dim> & dof1,
1864 *   const DoFHandler<dim> & dof2,
1865 *   const AffineConstraints<double> &constraints,
1866 *   const Vector<double> & solution,
1867 *   Vector<double> & u2)
1868 *   {
1869 *   assert(u2.size() == dof2.n_dofs() && "Incorrect input vector size!");
1870 *  
1871 *   u2 = 0.0;
1872 *  
1873 *   typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
1874 *   endc1 = dof1.end();
1875 *   typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
1876 *  
1877 *   for (; cell1 < endc1; ++cell1, ++cell2)
1878 *   {
1879 *   const auto &fe1 =
1880 *   dynamic_cast<const FE_Nedelec<dim> &>(cell1->get_fe());
1881 *   const auto &fe2 =
1882 *   dynamic_cast<const FE_Nedelec<dim> &>(cell2->get_fe());
1883 *  
1884 *   assert(fe1.degree < fe2.degree && "Incorrect usage of embed!");
1885 *  
1886 * @endcode
1887 *
1889 *
1890
1891 *
1892 *
1893
1894 *
1895 *
1896 * @code
1897 *   std::vector<unsigned int> embedding_dofs =
1898 *   fe2.get_embedding_dofs(fe1.degree);
1899 *   const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1900 *  
1901 *  
1904 *  
1905 *   local_dof_values_1.reinit(fe1.dofs_per_cell);
1906 *   cell1->get_dof_values(solution, local_dof_values_1);
1907 *  
1908 *   for (unsigned int i = 0; i < local_dof_values_1.size(); ++i)
1910 *  
1911 * @endcode
1912 *
1913 * Now set this changes to the global vector
1914 *
1915 * @code
1916 *   cell2->set_dof_values(local_dof_values_2, u2);
1917 *   }
1918 *  
1919 *   u2.compress(VectorOperation::insert);
1920 * @endcode
1921 *
1922 * Applies the constraints of the target finite element space
1923 *
1924 * @code
1925 *   constraints.distribute(u2);
1926 *   }
1927 *  
1928 *   template <int dim, bool report_dual>
1929 *   void
1931 *   const DoFHandler<dim> & dof1,
1932 *   const DoFHandler<dim> & dof2,
1933 *   const AffineConstraints<double> &constraints,
1934 *   const Vector<double> & solution,
1935 *   Vector<double> & u2)
1936 *   {
1937 * @endcode
1938 *
1939 * Maps from fe1 to fe2
1940 *
1941 * @code
1942 *   assert(u2.size() == dof2.n_dofs() && "Incorrect input vector size!");
1943 *  
1944 *   u2 = 0.0;
1945 *  
1946 *   typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
1947 *   endc1 = dof1.end();
1948 *   typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
1949 *  
1950 *   for (; cell1 < endc1; ++cell1, ++cell2)
1951 *   {
1952 *   const auto &fe1 =
1953 *   dynamic_cast<const FE_Nedelec<dim> &>(cell1->get_fe());
1954 *   const auto &fe2 =
1955 *   dynamic_cast<const FE_Nedelec<dim> &>(cell2->get_fe());
1956 *  
1957 *   assert(fe1.degree > fe2.degree && "Incorrect usage of extract!");
1958 *  
1959 * @endcode
1960 *
1962 *
1963 * @code
1964 *   std::vector<unsigned int> embedding_dofs =
1965 *   fe1.get_embedding_dofs(fe2.degree);
1966 *   const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1967 *  
1968 *  
1971 *  
1972 *   local_dof_values_1.reinit(fe1.dofs_per_cell);
1973 *   cell1->get_dof_values(solution, local_dof_values_1);
1974 *  
1975 *   for (unsigned int i = 0; i < local_dof_values_2.size(); ++i)
1977 *  
1978 * @endcode
1979 *
1980 * Now set this changes to the global vector
1981 *
1982 * @code
1983 *   cell2->set_dof_values(local_dof_values_2, u2);
1984 *   }
1985 *  
1986 *   u2.compress(VectorOperation::insert);
1987 * @endcode
1988 *
1989 * Applies the constraints of the target finite element space
1990 *
1991 * @code
1992 *   constraints.distribute(u2);
1993 *   }
1994 *   template <int dim, bool report_dual>
1995 *   void
1997 *   std::ofstream &os)
1998 *   {
1999 *   os << (*this->get_primal_eigenvalues())[0] << " "
2000 *   << (this->get_primal_DoFHandler())->n_dofs() << " "
2001 *   << (*this->get_dual_eigenvalues())[0] << " "
2002 *   << (this->get_dual_DoFHandler())->n_dofs() << std::endl;
2003 *   }
2004 *   template <int dim, bool report_dual>
2005 *   void
2007 *   std::ofstream &os)
2008 *   {
2009 *   os << qoi_error_estimate << std::endl;
2010 *   }
2011 *  
2012 *  
2016 *   template <int dim>
2017 *   class KellyErrorIndicator : public PrimalSolver<dim>
2018 *   {
2019 *   public:
2020 *   std::string
2021 *   name() const
2022 *   {
2023 *   return "Kelly";
2024 *   }
2025 *   void
2026 *   output_eigenvalue_data(std::ofstream &os);
2027 *   void
2028 *   output_qoi_error_estimates(std::ofstream &);
2029 *   KellyErrorIndicator(const std::string & prm_file,
2030 *   Triangulation<dim> &coarse_grid,
2031 *   const unsigned int &min_degree,
2032 *   const unsigned int &max_degree,
2033 *   const unsigned int &starting_degree);
2034 *  
2035 *   virtual unsigned int
2036 *   solve_problem() override;
2037 *  
2038 *   virtual void
2039 *   output_solution() override;
2040 *  
2042 *   get_FECollection();
2043 *  
2046 *  
2047 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2049 *  
2050 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2052 *  
2053 *   std::unique_ptr<std::vector<double>> &
2055 *  
2056 *  
2057 *   void
2059 *  
2061 *   get_DoFHandler();
2062 *  
2065 *  
2066 *   unsigned int
2067 *   get_max_degree()
2068 *   {
2069 *   return PrimalSolver<dim>::fe_collection->max_degree();
2070 *   }
2071 *   double qoi_error_estimate = 0;
2072 *  
2073 *   protected:
2074 *   void
2076 *  
2077 *   private:
2078 *   void
2079 *   prune_eigenpairs(const double &TOL);
2080 *  
2081 *   #if DEAL_II_VERSION_GTE(9, 6, 0)
2082 *   std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
2083 *   #else
2084 *   std::vector<const PETScWrappers::MPI::Vector *> eigenfunction_ptrs;
2085 *   #endif
2086 *   std::vector<const double *> eigenvalue_ptrs;
2087 *  
2088 *   std::vector<std::shared_ptr<Vector<float>>> errors;
2089 *   };
2090 *  
2091 *   template <int dim>
2093 *   const std::string & prm_file,
2094 *   Triangulation<dim> &coarse_grid,
2095 *   const unsigned int &min_degree,
2096 *   const unsigned int &max_degree,
2097 *   const unsigned int &starting_degree)
2098 *   : Base<dim>(prm_file, coarse_grid)
2100 *   coarse_grid,
2101 *   min_degree,
2102 *   max_degree,
2104 *   {}
2105 *  
2106 *   template <int dim>
2107 *   unsigned int
2109 *   {
2111 *   }
2112 *  
2113 *   template <int dim>
2116 *   {
2118 *   }
2119 *  
2120 *   template <int dim>
2123 *   {
2125 *   }
2126 *  
2127 *   template <int dim>
2128 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2130 *   {
2132 *   }
2133 *  
2134 *   template <int dim>
2135 *   std::unique_ptr<std::vector<double>> &
2137 *   {
2139 *   }
2140 *  
2141 *   template <int dim>
2142 *   std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2144 *   {
2146 *   }
2147 *  
2148 *   template <int dim>
2151 *   {
2153 *   }
2154 *  
2155 *   template <int dim>
2158 *   {
2160 *   }
2161 *  
2162 *   template <int dim>
2163 *   void
2165 *   {
2166 * @endcode
2167 *
2168 * This function does nothing for this error indicator
2169 *
2170 * @code
2171 *   return;
2172 *   }
2173 *  
2174 *   template <int dim>
2175 *   void
2177 *   {
2179 *   }
2180 *  
2181 *   template <int dim>
2182 *   void
2184 *   {
2185 *   unsigned int count = 0;
2186 *   for (size_t eigenpair_index = 0;
2187 *   eigenpair_index < this->eigenfunctions->size();
2188 *   ++eigenpair_index)
2189 *   {
2190 *   if (count >= this->n_eigenpairs)
2191 *   break;
2192 *   if (abs((*this->eigenvalues)[eigenpair_index]) < TOL)
2193 *   continue;
2194 *  
2195 *   eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
2196 *   eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
2197 *   }
2198 *   }
2199 *  
2200 *   template <int dim>
2201 *   void
2203 *   {
2204 *   std::cout << "Marking cells via Kelly indicator..." << std::endl;
2205 *   prune_eigenpairs(1e-9);
2206 * @endcode
2207 *
2208 * deallocate the errors vector
2209 *
2210 * @code
2211 *   errors.clear();
2212 *   for (size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2213 *   {
2214 *   errors.emplace_back(
2215 *   new Vector<float>(this->triangulation->n_active_cells()));
2216 *   }
2217 *   std::vector<Vector<float> *> estimated_error_per_cell(
2218 *   eigenfunction_ptrs.size());
2219 *   for (size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2220 *   {
2221 *   estimated_error_per_cell[i] = errors[i].get();
2222 *   }
2223 *  
2224 *   #if DEAL_II_VERSION_GTE(9, 6, 0)
2227 *   KellyErrorEstimator<dim>::estimate(this->dof_handler,
2228 *   *this->face_quadrature_collection,
2229 *   {},
2230 *   solution_view,
2231 *   error_view);
2232 *   #else
2233 *   KellyErrorEstimator<dim>::estimate(this->dof_handler,
2234 *   *this->face_quadrature_collection,
2235 *   {},
2238 *   #endif
2239 *  
2240 *   for (auto &error_vec : errors)
2241 *   {
2244 *  
2245 *   for (unsigned int i = 0; i < error_indicators.size(); ++i)
2246 *   error_indicators(i) += double(normalized_vec(i));
2247 *   }
2248 *   std::cout << "...Done!" << std::endl;
2249 *   }
2250 *   template <int dim>
2251 *   void
2253 *   {
2254 *   os << (*this->get_primal_eigenvalues())[0] << " "
2255 *   << (this->get_primal_DoFHandler())->n_dofs() << std::endl;
2256 *   }
2257 *   template <int dim>
2258 *   void
2260 *   {
2261 *   return;
2262 *   }
2263 *  
2264 *   } // namespace ErrorIndicators
2265 *  
2266 *  
2269 *   namespace RegularityIndicators
2270 *   {
2271 *   using namespace dealii;
2272 *  
2273 *   /* For the Legendre smoothness indicator*/
2274 *   /* Adapted from M. Fehling's smoothness_estimator.cc*/
2275 *   template <int dim>
2276 *   class LegendreInfo
2277 *   {};
2278 *  
2279 *   template <>
2280 *   class LegendreInfo<2>
2281 *   {
2282 *   public:
2283 *   std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
2284 *  
2285 *   hp::FECollection<2> *fe_collection = nullptr;
2286 *   DoFHandler<2> * dof_handler = nullptr;
2287 *  
2288 *   void
2289 *   initialization()
2290 *   {
2291 *   assert(fe_collection != nullptr && dof_handler != nullptr &&
2292 *   "A valid FECollection and DoFHandler must be accessible!");
2293 *  
2294 *   legendre_u = std::make_unique<FESeries::Legendre<2>>(
2296 *   legendre_v = std::make_unique<FESeries::Legendre<2>>(
2298 *  
2299 *   legendre_u->precalculate_all_transformation_matrices();
2300 *   legendre_v->precalculate_all_transformation_matrices();
2301 *   }
2302 *  
2303 *   template <class VectorType>
2304 *   void
2305 *   compute_coefficient_decay(const VectorType & eigenfunction,
2306 *   std::vector<double> &smoothness_indicators)
2307 *   {
2308 * @endcode
2309 *
2310 * Compute the coefficients for the u and v components of the solution
2311 * separately,
2312 *
2313 * @code
2316 *  
2318 *   *dof_handler,
2319 *   eigenfunction,
2320 *   smoothness_u);
2321 *  
2323 *   *dof_handler,
2324 *   eigenfunction,
2325 *   smoothness_v);
2326 *  
2327 *   for (unsigned int i = 0; i < smoothness_indicators.size(); ++i)
2328 *   {
2330 *   }
2331 *   }
2332 *   };
2333 *  
2334 *  
2337 *   template <int dim>
2338 *   class LegendreIndicator
2339 *   {
2340 *   public:
2341 *   void
2344 *  
2345 *   protected:
2346 *   template <class VectorType>
2347 *   void
2349 *   const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2350 *   const unsigned int & index_of_goal,
2351 *   std::vector<double> & smoothness_indicators);
2352 *  
2353 *   private:
2355 *   };
2356 *  
2357 *   template <int dim>
2358 *   void
2362 *   {
2363 *   legendre.fe_collection = fe_ptr;
2364 *   legendre.dof_handler = dof_ptr;
2365 *   this->legendre.initialization();
2366 *   }
2367 *  
2368 *   template <int dim>
2369 *   template <class VectorType>
2370 *   void
2372 *   const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2373 *   const unsigned int & index_of_goal,
2374 *   std::vector<double> & smoothness_indicators)
2375 *   {
2376 *   this->legendre.compute_coefficient_decay((*eigenfunctions)[index_of_goal],
2378 *   }
2379 *   } // namespace RegularityIndicators
2380 *  
2381 *  
2385 *   namespace Refinement
2386 *   {
2387 *   using namespace dealii;
2388 *   using namespace Maxwell;
2389 *  
2390 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2391 *   class Refiner : public ErrorIndicator, public RegularityIndicator
2392 *   {
2393 *   public:
2394 *   Refiner(const std::string & prm_file,
2395 *   Triangulation<dim> &coarse_grid,
2396 *   const unsigned int &min_degree,
2397 *   const unsigned int &max_degree,
2398 *   const unsigned int &starting_degree);
2399 *  
2400 *   void
2401 *   execute_refinement(const double &smoothness_threshold_fraction);
2402 *  
2403 *   virtual void
2404 *   output_solution() override;
2405 *  
2406 *   private:
2408 *   std::vector<double> smoothness_indicators;
2409 *   std::ofstream eigenvalues_out;
2410 *   std::ofstream error_estimate_out;
2411 *   };
2412 *  
2413 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2415 *   const std::string & prm_file,
2416 *   Triangulation<dim> &coarse_grid,
2417 *   const unsigned int &min_degree,
2418 *   const unsigned int &max_degree,
2419 *   const unsigned int &starting_degree)
2420 *   : Base<dim>(prm_file, coarse_grid)
2422 *   coarse_grid,
2423 *   min_degree,
2424 *   max_degree,
2427 *   {
2428 *   if (ErrorIndicator::name() == "DWR")
2429 *   {
2430 *   error_estimate_out.open("error_estimate.txt");
2431 *   error_estimate_out << std::setprecision(20);
2432 *   }
2433 *  
2434 *   eigenvalues_out.open("eigenvalues_" + ErrorIndicator::name() + "_out.txt");
2435 *   eigenvalues_out << std::setprecision(20);
2436 *   }
2437 *  
2438 * @endcode
2439 *
2441 *
2442 * @code
2443 *   template <int dim>
2444 *   class CurlPostprocessor : public DataPostprocessorScalar<dim>
2445 *   {
2446 *   public:
2449 *   {}
2450 *  
2451 *   virtual void
2452 *   evaluate_vector_field(
2454 *   std::vector<Vector<double>> &computed_quantities) const override
2455 *   {
2456 *   AssertDimension(input_data.solution_gradients.size(),
2457 *   computed_quantities.size());
2458 *   for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
2459 *   {
2460 *   computed_quantities[p](0) = input_data.solution_gradients[p][1][0] -
2461 *   input_data.solution_gradients[p][0][1];
2462 *   }
2463 *   }
2464 *   };
2465 *  
2466 *  
2473 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2474 *   void
2476 *   {
2478 *  
2479 *   DataOut<dim> data_out;
2481 *   data_out.attach_dof_handler(output_dof);
2482 *   Vector<double> fe_degrees(this->triangulation->n_active_cells());
2483 *   for (const auto &cell : output_dof.active_cell_iterators())
2484 *   fe_degrees(cell->active_cell_index()) =
2485 *   (*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
2486 *   .degree;
2487 *   data_out.add_data_vector(fe_degrees, "fe_degree");
2488 *  
2489 *   data_out.add_data_vector(estimated_error_per_cell, "error");
2490 *   Vector<double> smoothness_out(this->triangulation->n_active_cells());
2491 *   for (const auto &cell : output_dof.active_cell_iterators())
2492 *   {
2493 *   auto i = cell->active_cell_index();
2494 *   if (!cell->refine_flag_set() && !cell->coarsen_flag_set())
2495 *   smoothness_out(i) = -1;
2496 *   else
2498 *   }
2499 *   data_out.add_data_vector(smoothness_out, "smoothness");
2500 *   data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2501 *   std::string("eigenfunction_no_") +
2503 *   data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2504 *   curl_u);
2505 *  
2508 *  
2509 *   std::cout << "Number of DoFs: " << (this->get_primal_DoFHandler())->n_dofs()
2510 *   << std::endl;
2511 *  
2512 *  
2513 *   data_out.build_patches();
2514 *   std::ofstream output("eigenvectors-" + ErrorIndicator::name() + "-" +
2515 *   std::to_string(this->refinement_cycle) + +".vtu");
2516 *   data_out.write_vtu(output);
2517 *   }
2518 *  
2519 *  
2520 *  
2525 *   template <int dim, class ErrorIndicator, class RegularityIndicator>
2526 *   void
2528 *   const double &smoothness_threshold_fraction)
2529 *   {
2530 * @endcode
2531 *
2532 * First initialize the RegularityIndicator...
2533 * Depending on the limits set, this may take a while
2534 *
2535 * @code
2536 *   std::cout << "Initializing RegularityIndicator..." << std::endl;
2537 *   std::cout
2538 *   << "(This may take a while if the max expansion order is set too high)"
2539 *   << std::endl;
2542 *   std::cout << "Done!" << std::endl << "Starting Refinement..." << std::endl;
2543 *  
2544 *   for (unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
2545 *   {
2546 *   this->set_refinement_cycle(cycle);
2547 *   std::cout << "Cycle: " << cycle << std::endl;
2549 *   this->estimated_error_per_cell.reinit(
2550 *   this->triangulation->n_active_cells());
2551 *  
2553 *  
2554 * @endcode
2555 *
2558 *
2559 * @code
2562 *  
2563 *  
2565 *   *this->triangulation, estimated_error_per_cell, 1. / 5., 0.000);
2566 *  
2567 * @endcode
2568 *
2569 * Now get regularity indicators
2571 * depending on the regularity threshold...
2572 *
2573
2574 *
2575 *
2576 * @code
2578 *   std::vector<double>(this->triangulation->n_active_cells(),
2579 *   std::numeric_limits<double>::max());
2580 *   if (ErrorIndicator::PrimalSolver::min_degree !=
2581 *   ErrorIndicator::PrimalSolver::max_degree)
2584 * @endcode
2585 *
2586 * save data
2587 *
2588 * @code
2589 *   this->output_solution();
2591 *   unsigned int num_refined = 0, num_coarsened = 0;
2592 *   if (ErrorIndicator::PrimalSolver::min_degree !=
2593 *   ErrorIndicator::PrimalSolver::max_degree)
2594 *   {
2595 *   for (const auto &cell :
2596 *   ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2597 *   {
2598 *   if (cell->refine_flag_set())
2599 *   ++num_refined;
2600 *   if (cell->coarsen_flag_set())
2601 *   ++num_coarsened;
2602 *   if (cell->refine_flag_set() &&
2603 *   smoothness_indicators[cell->active_cell_index()] >
2605 *   static_cast<unsigned int>(cell->active_fe_index() + 1) <
2607 *   {
2608 *   cell->clear_refine_flag();
2609 *   cell->set_active_fe_index(cell->active_fe_index() + 1);
2610 *   }
2611 *   else if (cell->coarsen_flag_set() &&
2612 *   smoothness_indicators[cell->active_cell_index()] <
2614 *   cell->active_fe_index() != 0)
2615 *   {
2616 *   cell->clear_coarsen_flag();
2617 *  
2618 *   cell->set_active_fe_index(cell->active_fe_index() - 1);
2619 *   }
2620 * @endcode
2621 *
2623 *
2624 * @code
2625 *   else if (cell->refine_flag_set() && cell->diameter() < 5.0e-6)
2626 *   {
2627 *   cell->clear_refine_flag();
2628 *   if (static_cast<unsigned int>(cell->active_fe_index() + 1) <
2630 *   cell->set_active_fe_index(cell->active_fe_index() + 1);
2631 *   }
2632 *   }
2633 *   }
2634 *  
2635 * @endcode
2636 *
2638 *
2639 * @code
2640 *   double min_diameter = std::numeric_limits<double>::max();
2641 *   for (const auto &cell :
2642 *   ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2643 *   if (cell->diameter() < min_diameter)
2644 *   min_diameter = cell->diameter();
2645 *  
2646 *   std::cout << "Min diameter: " << min_diameter << std::endl;
2647 *  
2649 *  
2650 *   (this->triangulation)->execute_coarsening_and_refinement();
2651 *   }
2652 *   }
2653 *   } // namespace Refinement
2654 *  
2655 *   int
2656 *   main(int argc, char **argv)
2657 *   {
2658 *   try
2659 *   {
2660 *   using namespace dealii;
2661 *   using namespace Maxwell;
2662 *   using namespace Refinement;
2663 *   using namespace ErrorIndicators;
2664 *   using namespace RegularityIndicators;
2665 *  
2666 *  
2668 *  
2669 *  
2670 *   AssertThrow(
2672 *   ExcMessage("This program can only be run in serial, use ./maxwell-hp"));
2673 *  
2677 *  
2679 *   "maxwell-hp.prm",
2681 *   /*Minimum Degree*/ 2,
2682 *   /*Maximum Degree*/ 5,
2683 *   /*Starting Degree*/ 2);
2684 *  
2686 *   problem_DWR("maxwell-hp.prm",
2688 *   /*Minimum Degree*/ 2,
2689 *   /*Maximum Degree*/ 5,
2690 *   /*Starting Degree*/ 2);
2691 *  
2692 * @endcode
2693 *
2694 * The threshold for the hp-decision: too small -> not enough
2696 *
2697 * @code
2698 *   double smoothness_threshold = 0.75;
2699 *  
2700 *   std::cout << "Executing refinement for the Kelly strategy!" << std::endl;
2701 *   problem_Kelly.execute_refinement(smoothness_threshold);
2702 *   std::cout << "...Done with Kelly refinement strategy!" << std::endl;
2703 *   std::cout << "Executing refinement for the DWR strategy!" << std::endl;
2704 *   problem_DWR.execute_refinement(smoothness_threshold);
2705 *   std::cout << "...Done with DWR refinement strategy!" << std::endl;
2706 *   }
2707 *  
2708 *   catch (std::exception &exc)
2709 *   {
2710 *   std::cerr << std::endl
2711 *   << std::endl
2712 *   << "----------------------------------------------------"
2713 *   << std::endl;
2714 *   std::cerr << "Exception on processing: " << std::endl
2715 *   << exc.what() << std::endl
2716 *   << "Aborting!" << std::endl
2717 *   << "----------------------------------------------------"
2718 *   << std::endl;
2719 *  
2720 *   return 1;
2721 *   }
2722 *   catch (...)
2723 *   {
2724 *   std::cerr << std::endl
2725 *   << std::endl
2726 *   << "----------------------------------------------------"
2727 *   << std::endl;
2728 *   std::cerr << "Unknown exception!" << std::endl
2729 *   << "Aborting!" << std::endl
2730 *   << "----------------------------------------------------"
2731 *   << std::endl;
2732 *   return 1;
2733 *   }
2734 *  
2735 *   std::cout << std::endl << " Job done." << std::endl;
2736 *  
2737 *   return 0;
2738 *   }
2739 * @endcode
2740
2741
2742*/
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
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)
constexpr void clear()
friend class Tensor
Definition tensor.h:865
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3020
#define DEAL_II_VERSION_GTE(major, minor, subminor)
Definition config.h:349
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &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={})
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
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())
double diameter(const Triangulation< dim, spacedim > &tria)
@ valid
Iterator points to a valid object.
constexpr types::blas_int one
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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)
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
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)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:97
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474
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)
Definition hp.h:117
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14889
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
double legendre(unsigned int l, double x)
Definition cmath.h:65
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)