deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-89.h
Go to the documentation of this file.
1{},
347 /* interior_face_operation= */ face_function,
348 /* boundary_face_operation= */{},
349 dst, src);
350@endcode
351
352The code to do the same with FERemoteEvaluation is shown below.
353For brevity, we assume all boundary faces are somehow connected via
354non-conforming interfaces for FERemoteEvaluation.
355
356@code
357// Initialize FERemoteEvaluation: Note, that FERemoteEvaluation internally manages
358// the memory to store precomputed values. Therefore, FERemoteEvaluation
359// should be initialized only once to avoid frequent memory
360// allocation/deallocation. At this point, remote_communicator is assumed
361// to be initialized.
362FERemoteEvaluation<dim,Number> u_p_evaluator(remote_communicator);
363
364// Precompute the interpolated values of the finite element solution at all
365// the quadrature points outside the loop, invoking the vector entries and
366// respective basis function at possibly remote MPI processes before communication.
367u_p_evaluator.gather_evaluate(src, EvaluationFlags::values);
368
369const auto boundary_function =
370 [&](const MatrixFree &data, VectorType &dst, const VectorType &src,
371 const std::pair<unsigned int, unsigned int> face_range) {
372
373 FEFaceEvaluation phi_m(data, true);
374 // To access the values in a thread safe way each thread has
375 // to create a own accessor object. A small helper function
376 // provides the accessor.
377 internal::PrecomputedEvaluationDataAccessor u_p = u_p_evaluator.get_data_accessor();
378
379 for (unsigned int f = face_range.first; f < face_range.second; ++f)
380 {
381 phi_m.reinit(f);
382 u_p.reinit(f);
383
384 for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
385 phi_m.submit_value(u_p.get_value(q), q); // access values with phi_p
386
387 phi_m.integrate_scatter(EvaluationFlags::values, dst);
388 }
389 };
390
391matrix_free.template loop<VectorType, VectorType>({}, {}, boundary_function, dst, src);
392@endcode
393The object @c remote_communicator is of type `FERemoteEvaluationCommunicator`
394and assumed to be correctly initialized prior to the above code snippet.
395`FERemoteEvaluationCommunicator` internally manages the update of ghost values
396over non-matching interfaces and keeps track of the mapping between quadrature
397point index and corresponding values/gradients. As mentioned above, the update
398of the values/gradients happens <i>before</i> the actual matrix-free loop.
400differently for the given template parameter @c value_type. If we want to access
401values at arbitrary points (e.g. in combination with @c FEPointEvaluation), then
402we need to choose @c value_type=Number. If the values are defined at quadrature
403points of a @c FEEvaluation object it is possible to get the values at the
404quadrature points of <i>batches</i> and we need to choose
405@c value_type=VectorizedArray<Number>.
406
407
408<a name="step_89-Overviewofthetestcase"></a><h3>Overview of the test case</h3>
409
410
411In this program, we implemented both the point-to-point interpolation and
412Nitsche-type mortaring mentioned in the introduction.
413
414At first we are considering the test case of a vibrating membrane, see e.g.
415@cite nguyen2011high. Standing waves of length @f$\lambda=2/M@f$ are oscillating
416with a time period of @f$T=2 / (M \sqrt{d} c)@f$ where @f$d@f$ is the dimension of the
417space in which our domain is located and @f$M@f$ is the number of modes per meter,
418i.e. the number of half-waves per meter. The corresponding analytical solution
419reads as
420
421@f{align*}{
422 p &=\cos(M \sqrt{d} \pi c t)\prod_{i=1}^{d} \sin(M \pi x_i),\\
423 u_i&=-\frac{\sin(M \sqrt{d} \pi c t)}{\sqrt{d}\rho c}
424 \cos(M \pi x_i)\prod_{j=1,j\neq i}^{d} \sin(M \pi x_j),
425@f}
426
427For simplicity, we are using homogeneous pressure Dirichlet boundary conditions
428within this tutorial. To be able to do so we have to tailor the domain size as
429well as the number of modes to conform with the homogeneous pressure Dirichlet
430boundary conditions. Within this tutorial we are using @f$M=10@f$ and a domain
431@f$\Omega=(0,1)^2@f$. The domain will be meshed so that the left and right parts of
432the domain consist of separate meshes that do not match at the interface.
433
434As will become clear from the results, the point-to-point interpolation will
435result in aliasing, which can be resolved using Nitsche-type mortaring.
436
437In a more realistic second example, we apply this implementation to a test case
438in which a wave is propagating from one fluid into another fluid. The speed of
439sound in the left part of the domain is @f$c=1@f$ and in the right part it is @f$c=3@f$.
440Since the wavelength is directly proportional to the speed of sound, three times
441larger elements can be used in the right part of the domain to resolve waves up
442to the same frequency. A test case like this has been simulated with a different
443domain and different initial conditions, e.g., in @cite bangerth2010adaptive.
444 *
445 *
446 * <a name="step_89-CommProg"></a>
447 * <h1> The commented program</h1>
448 *
449 *
450 * <a name="step_89-Includefiles"></a>
451 * <h3>Include files</h3>
452 *
453
454 *
455 * The program starts with including all the relevant header files.
456 *
457 * @code
458 *   #include <deal.II/base/conditional_ostream.h>
459 *   #include <deal.II/base/mpi.h>
460 *  
461 *   #include <deal.II/distributed/tria.h>
462 *  
463 *   #include <deal.II/fe/fe_dgq.h>
464 *   #include <deal.II/fe/fe_system.h>
465 *   #include <deal.II/fe/fe_values.h>
466 *   #include <deal.II/fe/mapping_q1.h>
467 *  
468 *   #include <deal.II/grid/grid_generator.h>
469 *   #include <deal.II/grid/grid_tools.h>
470 *   #include <deal.II/grid/grid_tools_cache.h>
471 *  
472 *   #include <deal.II/matrix_free/fe_evaluation.h>
473 *   #include <deal.II/matrix_free/matrix_free.h>
474 *   #include <deal.II/matrix_free/operators.h>
475 *  
476 *   #include <deal.II/non_matching/mapping_info.h>
477 *  
478 *   #include <deal.II/numerics/data_out.h>
479 *   #include <deal.II/numerics/vector_tools.h>
480 * @endcode
481 *
482 * The following header file provides the class FERemoteEvaluation, which allows
483 * to access values and/or gradients at remote triangulations similar to
484 * FEEvaluation.
485 *
486 * @code
487 *   #include <deal.II/matrix_free/fe_remote_evaluation.h>
488 *  
489 * @endcode
490 *
491 * We pack everything that is specific for this program into a namespace
492 * of its own.
493 *
494 * @code
495 *   namespace Step89
496 *   {
497 *   using namespace dealii;
498 *  
499 * @endcode
500 *
501 *
502 * <a name="step_89-Initialconditionsforvibratingmembrane"></a>
503 * <h3>Initial conditions for vibrating membrane</h3>
504 *
505
506 *
507 * In the following, let us first define a function that provides
508 * the initial condition for the vibrating membrane test case. It
509 * implements both the initial pressure (component 0) and velocity
510 * (components 1 to 1 + dim).
511 *
512
513 *
514 * There is also a function that computes the duration of one
515 * oscillation.
516 *
517 * @code
518 *   template <int dim>
519 *   class InitialConditionVibratingMembrane : public Function<dim>
520 *   {
521 *   public:
522 *   InitialConditionVibratingMembrane(const double modes);
523 *  
524 *   double value(const Point<dim> &p, const unsigned int comp) const final;
525 *  
526 *   double get_period_duration(const double speed_of_sound) const;
527 *  
528 *   private:
529 *   const double M;
530 *   };
531 *  
532 *   template <int dim>
533 *   InitialConditionVibratingMembrane<dim>::InitialConditionVibratingMembrane(
534 *   const double modes)
535 *   : Function<dim>(dim + 1, 0.0)
536 *   , M(modes)
537 *   {
538 *   static_assert(dim == 2, "Only implemented for dim==2");
539 *   }
540 *  
541 *   template <int dim>
542 *   double
543 *   InitialConditionVibratingMembrane<dim>::value(const Point<dim> &p,
544 *   const unsigned int comp) const
545 *   {
546 *   if (comp == 0)
547 *   return std::sin(M * numbers::PI * p[0]) *
548 *   std::sin(M * numbers::PI * p[1]);
549 *  
550 *   return 0.0;
551 *   }
552 *  
553 *   template <int dim>
554 *   double InitialConditionVibratingMembrane<dim>::get_period_duration(
555 *   const double speed_of_sound) const
556 *   {
557 *   return 2.0 / (M * std::sqrt(dim) * speed_of_sound);
558 *   }
559 *  
560 * @endcode
561 *
562 *
563 * <a name="step_89-Gausspulse"></a>
564 * <h3>Gauss pulse</h3>
565 *
566
567 *
568 * Next up is a function that provides the values of a pressure
569 * Gauss pulse. As with the previous function, it implements both
570 * the initial pressure (component 0) and velocity (components 1 to
571 * 1 + dim).
572 *
573 * @code
574 *   template <int dim>
575 *   class GaussPulse : public Function<dim>
576 *   {
577 *   public:
578 *   GaussPulse(const double shift_x, const double shift_y);
579 *  
580 *   double value(const Point<dim> &p, const unsigned int comp) const final;
581 *  
582 *   private:
583 *   const double shift_x;
584 *   const double shift_y;
585 *   };
586 *  
587 *   template <int dim>
588 *   GaussPulse<dim>::GaussPulse(const double shift_x, const double shift_y)
589 *   : Function<dim>(dim + 1, 0.0)
590 *   , shift_x(shift_x)
591 *   , shift_y(shift_y)
592 *   {
593 *   static_assert(dim == 2, "Only implemented for dim==2");
594 *   }
595 *  
596 *   template <int dim>
597 *   double GaussPulse<dim>::value(const Point<dim> &p,
598 *   const unsigned int comp) const
599 *   {
600 *   if (comp == 0)
601 *   return std::exp(-1000.0 * ((std::pow(p[0] - shift_x, 2)) +
602 *   (std::pow(p[1] - shift_y, 2))));
603 *  
604 *   return 0.0;
605 *   }
606 *  
607 * @endcode
608 *
609 *
610 * <a name="step_89-Helperfunctions"></a>
611 * <h3>Helper functions</h3>
612 *
613
614 *
615 * The following namespace contains free helper functions that are used in the
616 * tutorial.
617 *
618 * @code
619 *   namespace HelperFunctions
620 *   {
621 * @endcode
622 *
623 * Helper function to check if a boundary ID is related to a non-matching
624 * face. A @c std::set that contains all non-matching boundary IDs is
625 * handed over in addition to the face ID under question.
626 *
627 * @code
628 *   bool is_non_matching_face(
629 *   const std::set<types::boundary_id> &non_matching_face_ids,
630 *   const types::boundary_id face_id)
631 *   {
632 *   return non_matching_face_ids.find(face_id) != non_matching_face_ids.end();
633 *   }
634 *  
635 * @endcode
636 *
637 * Helper function to set the initial conditions for the vibrating membrane
638 * test case.
639 *
640 * @code
641 *   template <int dim, typename Number, typename VectorType>
642 *   void set_initial_condition(MatrixFree<dim, Number> matrix_free,
643 *   const Function<dim> &initial_solution,
644 *   VectorType &dst)
645 *   {
646 *   VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
647 *   matrix_free.get_dof_handler(),
648 *   initial_solution,
649 *   dst);
650 *   }
651 *  
652 * @endcode
653 *
654 * Helper function to compute the time step size according to the CFL
655 * condition.
656 *
657 * @code
658 *   double
659 *   compute_dt_cfl(const double hmin, const unsigned int degree, const double c)
660 *   {
661 *   return hmin / (std::pow(degree, 1.5) * c);
662 *   }
663 *  
664 * @endcode
665 *
666 * Helper function that writes vtu output.
667 *
668 * @code
669 *   template <typename VectorType, int dim>
670 *   void write_vtu(const VectorType &solution,
671 *   const DoFHandler<dim> &dof_handler,
672 *   const Mapping<dim> &mapping,
673 *   const unsigned int degree,
674 *   const std::string &name_prefix)
675 *   {
676 *   DataOut<dim> data_out;
677 *   DataOutBase::VtkFlags flags;
678 *   flags.write_higher_order_cells = true;
679 *   data_out.set_flags(flags);
680 *  
681 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
682 *   interpretation(
684 *   std::vector<std::string> names(dim + 1, "U");
685 *  
687 *   names[0] = "P";
688 *  
689 *   data_out.add_data_vector(dof_handler, solution, names, interpretation);
690 *  
691 *   data_out.build_patches(mapping, degree, DataOut<dim>::curved_inner_cells);
692 *   data_out.write_vtu_in_parallel(name_prefix + ".vtu",
693 *   dof_handler.get_mpi_communicator());
694 *   }
695 *   } // namespace HelperFunctions
696 *  
697 * @endcode
698 *
699 *
700 * <a name="step_89-Materialparameterdescription"></a>
701 * <h3>Material parameter description</h3>
702 *
703
704 *
705 * The following class stores the information if the fluid is
706 * homogeneous as well as the material properties at every cell.
707 * This class helps access the correct values without accessing a
708 * large vector of materials in the homogeneous case.
709 *
710
711 *
712 * The class is provided a map from material ids to material
713 * properties (given as a pair of values for the speed of sound and
714 * the density). If the map has only one entry, the material is
715 * homogeneous -- using the same values everywhere -- and we can
716 * remember that fact in the `homogeneous` member variable and use it
717 * to optimize some code paths below. If the material is not
718 * homogeneous, we will fill a vector with the correct materials for
719 * each batch of cells; this information can then be access via
721 *
722
723 *
724 * As is usual when working with the MatrixFree framework, we will
725 * not only access material parameters at a single site, but for
726 * whole "batches" of cells. As a consequence, the functions below
727 * returned VectorizedArray objects for a batch at a time.
728 *
729 * @code
730 *   template <typename Number>
731 *   class CellwiseMaterialData
732 *   {
733 *   public:
734 *   template <int dim>
735 *   CellwiseMaterialData(
736 *   const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
737 *   const std::map<types::material_id, std::pair<double, double>>
738 *   &material_id_map)
739 *   : homogeneous(material_id_map.size() == 1)
740 *   {
741 *   Assert(material_id_map.size() > 0,
742 *   ExcMessage("No materials given to CellwiseMaterialData."));
743 *  
744 *   if (homogeneous)
745 *   {
746 *   speed_of_sound_homogeneous = material_id_map.begin()->second.first;
747 *   density_homogeneous = material_id_map.begin()->second.second;
748 *   }
749 *   else
750 *   {
751 *   const auto n_cell_batches =
752 *   matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
753 *  
754 *   speed_of_sound.resize(n_cell_batches);
755 *   density.resize(n_cell_batches);
756 *  
757 *   for (unsigned int cell = 0; cell < n_cell_batches; ++cell)
758 *   {
759 *   speed_of_sound[cell] = 1.;
760 *   density[cell] = 1.;
761 *   for (unsigned int v = 0;
762 *   v < matrix_free.n_active_entries_per_cell_batch(cell);
763 *   ++v)
764 *   {
765 *   const auto material_id =
766 *   matrix_free.get_cell_iterator(cell, v)->material_id();
767 *  
768 *   speed_of_sound[cell][v] =
769 *   material_id_map.at(material_id).first;
770 *   density[cell][v] = material_id_map.at(material_id).second;
771 *   }
772 *   }
773 *   }
774 *   }
775 *  
776 *   bool is_homogeneous() const
777 *   {
778 *   return homogeneous;
779 *   }
780 *  
781 *   const AlignedVector<VectorizedArray<Number>> &get_speed_of_sound() const
782 *   {
783 *   Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()."));
784 *   return speed_of_sound;
785 *   }
786 *  
787 *   const AlignedVector<VectorizedArray<Number>> &get_density() const
788 *   {
789 *   Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()."));
790 *   return density;
791 *   }
792 *  
793 *   VectorizedArray<Number> get_homogeneous_speed_of_sound() const
794 *   {
795 *   Assert(homogeneous, ExcMessage("Use get_speed_of_sound()."));
796 *   return speed_of_sound_homogeneous;
797 *   }
798 *  
799 *   VectorizedArray<Number> get_homogeneous_density() const
800 *   {
801 *   Assert(homogeneous, ExcMessage("Use get_density()."));
802 *   return density_homogeneous;
803 *   }
804 *  
805 *   private:
806 *   const bool homogeneous;
807 *  
808 *   /* Materials in the inhomogeneous case. */
809 *   AlignedVector<VectorizedArray<Number>> speed_of_sound;
811 *  
812 *   /* Materials in the homogeneous case. */
813 *   VectorizedArray<Number> speed_of_sound_homogeneous;
814 *   VectorizedArray<Number> density_homogeneous;
815 *   };
816 *  
817 * @endcode
818 *
819 * To be able to access the material data in every cell in a thread
820 * safe way, the following class @c MaterialEvaluation is
821 * used. Similar to @c FEEvaluation, functions below will create
822 * their own instances of this class; thus, there can be no race
823 * conditions even if these functions run multiple times in
824 * parallel. For inhomogeneous materials, a @c reinit_cell() or @c
825 * reinit_face() function is used to set the correct material at the
826 * current cell batch. In the homogeneous case the @c _reinit()
827 * functions don't have to reset the materials.
828 *
829 * @code
830 *   template <int dim, typename Number>
831 *   class MaterialEvaluation
832 *   {
833 *   public:
834 *   MaterialEvaluation(
835 *   const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
836 *   const CellwiseMaterialData<Number> &material_data)
837 *   : phi(matrix_free)
838 *   , phi_face(matrix_free, true)
839 *   , material_data(material_data)
840 *   {
841 *   if (material_data.is_homogeneous())
842 *   {
843 *   speed_of_sound = material_data.get_homogeneous_speed_of_sound();
844 *   density = material_data.get_homogeneous_density();
845 *   }
846 *   }
847 *  
848 *   bool is_homogeneous() const
849 *   {
850 *   return material_data.is_homogeneous();
851 *   }
852 *  
853 * @endcode
854 *
855 * The following two functions update the data for the current
856 * cell or face, given a cell batch index. If the material is
857 * homogeneous, there is nothing to do. Otherwise, we reinit the
858 * FEEvaluation object and store the data for the current object.
859 *
860 * @code
861 *   void reinit_cell(const unsigned int cell)
862 *   {
863 *   if (!material_data.is_homogeneous())
864 *   {
865 *   phi.reinit(cell);
866 *   speed_of_sound =
867 *   phi.read_cell_data(material_data.get_speed_of_sound());
868 *   density = phi.read_cell_data(material_data.get_density());
869 *   }
870 *   }
871 *  
872 *   void reinit_face(const unsigned int face)
873 *   {
874 *   if (!material_data.is_homogeneous())
875 *   {
876 *   phi_face.reinit(face);
877 *   speed_of_sound =
878 *   phi_face.read_cell_data(material_data.get_speed_of_sound());
879 *   density = phi_face.read_cell_data(material_data.get_density());
880 *   }
881 *   }
882 *  
883 * @endcode
884 *
885 * The following functions then return the speed of sound and
886 * density for the current cell batch.
887 *
888 * @code
889 *   VectorizedArray<Number> get_speed_of_sound() const
890 *   {
891 *   return speed_of_sound;
892 *   }
893 *  
894 *   VectorizedArray<Number> get_density() const
895 *   {
896 *   return density;
897 *   }
898 *  
899 *   private:
900 *   /* Members needed for the inhomogeneous case. */
901 *   FEEvaluation<dim, -1, 0, 1, Number> phi;
902 *   FEFaceEvaluation<dim, -1, 0, 1, Number> phi_face;
903 *  
904 *   /* Material defined at every cell. */
905 *   const CellwiseMaterialData<Number> &material_data;
906 *  
907 *   /* Materials at current cell. */
908 *   VectorizedArray<Number> speed_of_sound;
909 *   VectorizedArray<Number> density;
910 *   };
911 *  
912 *  
913 * @endcode
914 *
915 *
916 * <a name="step_89-Boundaryconditions"></a>
917 * <h3>Boundary conditions</h3>
918 *
919
920 *
921 * To be able to use the same kernel, for all face integrals we define
922 * a class that returns the needed values at boundaries. In this tutorial
923 * homogeneous pressure Dirichlet boundary conditions are applied via
924 * the mirror principle, i.e. @f$p_h^+=-p_h^- + 2g@f$ with @f$g=0@f$.
925 *
926 * @code
927 *   template <int dim, typename Number>
928 *   class BCEvaluationP
929 *   {
930 *   public:
931 *   BCEvaluationP(const FEFaceEvaluation<dim, -1, 0, 1, Number> &pressure_m)
932 *   : pressure_m(pressure_m)
933 *   {}
934 *  
935 *   typename FEFaceEvaluation<dim, -1, 0, 1, Number>::value_type
936 *   get_value(const unsigned int q) const
937 *   {
938 *   return -pressure_m.get_value(q);
939 *   }
940 *  
941 *   private:
942 *   const FEFaceEvaluation<dim, -1, 0, 1, Number> &pressure_m;
943 *   };
944 *  
945 * @endcode
946 *
947 * We don't have to apply boundary conditions for the velocity, i.e.
948 * @f$\mathbf{u}_h^+=\mathbf{u}_h^-@f$.
949 *
950 * @code
951 *   template <int dim, typename Number>
952 *   class BCEvaluationU
953 *   {
954 *   public:
955 *   BCEvaluationU(const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m)
956 *   : velocity_m(velocity_m)
957 *   {}
958 *  
959 *   typename FEFaceEvaluation<dim, -1, 0, dim, Number>::value_type
960 *   get_value(const unsigned int q) const
961 *   {
962 *   return velocity_m.get_value(q);
963 *   }
964 *  
965 *   private:
966 *   const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m;
967 *   };
968 *  
969 * @endcode
970 *
971 *
972 * <a name="step_89-Acousticoperator"></a>
973 * <h3>Acoustic operator</h3>
974 *
975
976 *
977 * The following class then defines the acoustic operator. The class is
978 * heavily based on matrix-free methods. For a better understanding in
979 * matrix-free methods please refer to @ref step_67 "step-67".
980 *
981
982 *
983 * At the top we define a flag that decides whether we want to use
984 * mortaring. If the remote evaluators are set up with a
985 * VectorizedArray we are using point-to-point interpolation;
986 * otherwise we make use of Nitsche-type mortaring. The decision is
987 * made using `std::is_floating_point_v`, which is a variable that
988 * is true if the template argument is a floating point type (such
989 * as `float` or `double`) and false otherwise (in particular if the
990 * template argument is a vectorized array type).
991 *
992 * @code
993 *   template <int dim, typename Number, typename remote_value_type>
994 *   class AcousticOperator
995 *   {
996 *   static constexpr bool use_mortaring =
997 *   std::is_floating_point_v<remote_value_type>;
998 *  
999 *   public:
1000 * @endcode
1001 *
1002 * In case of Nitsche-type mortaring, `NonMatching::MappingInfo` has to
1003 * be provided in the constructor.
1004 *
1005 * @code
1006 *   AcousticOperator(
1007 *   const MatrixFree<dim, Number> &matrix_free,
1008 *   std::shared_ptr<CellwiseMaterialData<Number>> material_data,
1009 *   const std::set<types::boundary_id> &remote_face_ids,
1010 *   std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1011 *   pressure_r_eval,
1012 *   std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
1013 *   velocity_r_eval,
1014 *   std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> c_r_eval,
1015 *   std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> rho_r_eval,
1016 *   std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> nm_info =
1017 *   nullptr)
1018 *   : matrix_free(matrix_free)
1019 *   , material_data(material_data)
1020 *   , remote_face_ids(remote_face_ids)
1021 *   , pressure_r_eval(pressure_r_eval)
1022 *   , velocity_r_eval(velocity_r_eval)
1023 *   , c_r_eval(c_r_eval)
1024 *   , rho_r_eval(rho_r_eval)
1025 *   , nm_mapping_info(nm_info)
1026 *   {
1027 *   if (use_mortaring)
1028 *   Assert(nm_info,
1029 *   ExcMessage(
1030 *   "In case of Nitsche-type mortaring NonMatching::MappingInfo \
1031 *   has to be provided."));
1032 *   }
1033 *  
1034 * @endcode
1035 *
1036 * The following function then evaluates the acoustic operator.
1037 * It first updates the precomputed values in corresponding the
1038 * FERemoteEvaluation objects. The material parameters do not change and
1039 * thus, we do not have to update precomputed values in @c c_r_eval and @c
1040 * rho_r_eval.
1041 *
1042
1043 *
1044 * It then either performs a matrix-free loop with Nitsche-type
1045 * mortaring at non-matching faces (if `use_mortaring` is true) or
1046 * with point-to-point interpolation at non-matching faces (in the
1047 * `else` branch). The difference is only in the third argument to
1048 * the loop function, denoting the function object that is
1049 * executed at boundary faces.
1050 *
1051 * @code
1052 *   template <typename VectorType>
1053 *   void evaluate(VectorType &dst, const VectorType &src) const
1054 *   {
1055 *   pressure_r_eval->gather_evaluate(src, EvaluationFlags::values);
1056 *   velocity_r_eval->gather_evaluate(src, EvaluationFlags::values);
1057 *  
1058 *   if constexpr (use_mortaring)
1059 *   {
1060 *   matrix_free.loop(
1061 *   &AcousticOperator::local_apply_cell<VectorType>,
1062 *   &AcousticOperator::local_apply_face<VectorType>,
1063 *   &AcousticOperator::local_apply_boundary_face_mortaring<VectorType>,
1064 *   this,
1065 *   dst,
1066 *   src,
1067 *   true,
1068 *   MatrixFree<dim, Number>::DataAccessOnFaces::values,
1069 *   MatrixFree<dim, Number>::DataAccessOnFaces::values);
1070 *   }
1071 *   else
1072 *   {
1073 *   matrix_free.loop(
1074 *   &AcousticOperator::local_apply_cell<VectorType>,
1075 *   &AcousticOperator::local_apply_face<VectorType>,
1076 *   &AcousticOperator::local_apply_boundary_face_point_to_point<
1077 *   VectorType>,
1078 *   this,
1079 *   dst,
1080 *   src,
1081 *   true,
1082 *   MatrixFree<dim, Number>::DataAccessOnFaces::values,
1083 *   MatrixFree<dim, Number>::DataAccessOnFaces::values);
1084 *   }
1085 *   }
1086 *  
1087 * @endcode
1088 *
1089 * In the `private` section of the class, we then define the
1090 * functions that evaluate volume, interior face, and boundary
1091 * face integrals. The concrete terms these functions evaluate are
1092 * stated in the documentation at the top of this tutorial
1093 * program. Each of these functions has its own object of type
1094 * `MaterialEvaluation` that provides access to the material at
1095 * each cell or face.
1096 *
1097 * @code
1098 *   private:
1099 *   template <typename VectorType>
1100 *   void local_apply_cell(
1101 *   const MatrixFree<dim, Number> &matrix_free,
1102 *   VectorType &dst,
1103 *   const VectorType &src,
1104 *   const std::pair<unsigned int, unsigned int> &cell_range) const
1105 *   {
1106 *   FEEvaluation<dim, -1, 0, 1, Number> pressure(matrix_free, 0, 0, 0);
1107 *   FEEvaluation<dim, -1, 0, dim, Number> velocity(matrix_free, 0, 0, 1);
1108 *  
1109 *   MaterialEvaluation material(matrix_free, *material_data);
1110 *  
1111 *   for (unsigned int cell = cell_range.first; cell < cell_range.second;
1112 *   ++cell)
1113 *   {
1114 *   velocity.reinit(cell);
1115 *   pressure.reinit(cell);
1116 *  
1117 *   pressure.gather_evaluate(src, EvaluationFlags::gradients);
1118 *   velocity.gather_evaluate(src, EvaluationFlags::gradients);
1119 *  
1120 * @endcode
1121 *
1122 * Get the materials on the corresponding cell. Since we
1123 * introduced @c MaterialEvaluation we can write the code
1124 * independent of whether the material is homogeneous or
1125 * inhomogeneous.
1126 *
1127 * @code
1128 *   material.reinit_cell(cell);
1129 *   const auto c = material.get_speed_of_sound();
1130 *   const auto rho = material.get_density();
1131 *  
1132 *   for (const unsigned int q : pressure.quadrature_point_indices())
1133 *   {
1134 *   pressure.submit_value(rho * c * c * velocity.get_divergence(q),
1135 *   q);
1136 *   velocity.submit_value(1.0 / rho * pressure.get_gradient(q), q);
1137 *   }
1138 *  
1139 *   pressure.integrate_scatter(EvaluationFlags::values, dst);
1140 *   velocity.integrate_scatter(EvaluationFlags::values, dst);
1141 *   }
1142 *   }
1143 *  
1144 * @endcode
1145 *
1146 * This next function evaluates the fluxes at faces between cells with the
1147 * same material. If boundary faces are under consideration fluxes into
1148 * neighboring faces do not have to be considered which is enforced via
1149 * `weight_neighbor=false`. For non-matching faces the fluxes into
1150 * neighboring faces are not considered as well. This is because we iterate
1151 * over each side of the non-matching face separately (similar to a cell
1152 * centric loop).
1153 *
1154
1155 *
1156 * In this and following functions, we also introduce the factors
1157 * @f$\tau@f$ and @f$\gamma@f$ that are computed from @f$\rho@f$ and @f$c@f$ along
1158 * interfaces and that appear in the bilinear forms. Their
1159 * definitions are provided in the introduction.
1160 *
1161 * @code
1162 *   template <bool weight_neighbor,
1163 *   typename InternalFaceIntegratorPressure,
1164 *   typename InternalFaceIntegratorVelocity,
1165 *   typename ExternalFaceIntegratorPressure,
1166 *   typename ExternalFaceIntegratorVelocity>
1167 *   inline DEAL_II_ALWAYS_INLINE void evaluate_face_kernel(
1168 *   InternalFaceIntegratorPressure &pressure_m,
1169 *   InternalFaceIntegratorVelocity &velocity_m,
1170 *   ExternalFaceIntegratorPressure &pressure_p,
1171 *   ExternalFaceIntegratorVelocity &velocity_p,
1172 *   const typename InternalFaceIntegratorPressure::value_type c,
1173 *   const typename InternalFaceIntegratorPressure::value_type rho) const
1174 *   {
1175 *   const auto tau = 0.5 * rho * c;
1176 *   const auto gamma = 0.5 / (rho * c);
1177 *  
1178 *   for (const unsigned int q : pressure_m.quadrature_point_indices())
1179 *   {
1180 *   const auto n = pressure_m.normal_vector(q);
1181 *   const auto pm = pressure_m.get_value(q);
1182 *   const auto um = velocity_m.get_value(q);
1183 *  
1184 *   const auto pp = pressure_p.get_value(q);
1185 *   const auto up = velocity_p.get_value(q);
1186 *  
1187 * @endcode
1188 *
1189 * Compute homogeneous local Lax-Friedrichs fluxes and submit the
1190 * corresponding values to the integrators.
1191 *
1192 * @code
1193 *   const auto momentum_flux =
1194 *   0.5 * (pm + pp) + 0.5 * tau * (um - up) * n;
1195 *   velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
1196 *   if constexpr (weight_neighbor)
1197 *   velocity_p.submit_value(1.0 / rho * (momentum_flux - pp) * (-n), q);
1198 *  
1199 *   const auto mass_flux = 0.5 * (um + up) + 0.5 * gamma * (pm - pp) * n;
1200 *   pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
1201 *   if constexpr (weight_neighbor)
1202 *   pressure_p.submit_value(rho * c * c * (mass_flux - up) * (-n), q);
1203 *   }
1204 *   }
1205 *  
1206 * @endcode
1207 *
1208 * This function evaluates the fluxes at faces between cells with different
1209 * materials. This can only happen over non-matching interfaces. Therefore,
1210 * it is implicitly known that `weight_neighbor=false` and we can omit the
1211 * parameter.
1212 *
1213 * @code
1214 *   template <typename InternalFaceIntegratorPressure,
1215 *   typename InternalFaceIntegratorVelocity,
1216 *   typename ExternalFaceIntegratorPressure,
1217 *   typename ExternalFaceIntegratorVelocity,
1218 *   typename MaterialIntegrator>
1219 *   void evaluate_face_kernel_inhomogeneous(
1220 *   InternalFaceIntegratorPressure &pressure_m,
1221 *   InternalFaceIntegratorVelocity &velocity_m,
1222 *   const ExternalFaceIntegratorPressure &pressure_p,
1223 *   const ExternalFaceIntegratorVelocity &velocity_p,
1224 *   const typename InternalFaceIntegratorPressure::value_type c,
1225 *   const typename InternalFaceIntegratorPressure::value_type rho,
1226 *   const MaterialIntegrator &c_r,
1227 *   const MaterialIntegrator &rho_r) const
1228 *   {
1229 *   const auto tau_m = 0.5 * rho * c;
1230 *   const auto gamma_m = 0.5 / (rho * c);
1231 *  
1232 *   for (const unsigned int q : pressure_m.quadrature_point_indices())
1233 *   {
1234 *   const auto c_p = c_r.get_value(q);
1235 *   const auto rho_p = rho_r.get_value(q);
1236 *   const auto tau_p = 0.5 * rho_p * c_p;
1237 *   const auto gamma_p = 0.5 / (rho_p * c_p);
1238 *   const auto tau_sum_inv = 1.0 / (tau_m + tau_p);
1239 *   const auto gamma_sum_inv = 1.0 / (gamma_m + gamma_p);
1240 *  
1241 *   const auto n = pressure_m.normal_vector(q);
1242 *   const auto pm = pressure_m.get_value(q);
1243 *   const auto um = velocity_m.get_value(q);
1244 *  
1245 *   const auto pp = pressure_p.get_value(q);
1246 *   const auto up = velocity_p.get_value(q);
1247 *  
1248 *  
1249 * @endcode
1250 *
1251 * Compute inhomogeneous fluxes and submit the corresponding values
1252 * to the integrators.
1253 *
1254 * @code
1255 *   const auto momentum_flux =
1256 *   pm - tau_m * tau_sum_inv * (pm - pp) +
1257 *   tau_m * tau_p * tau_sum_inv * (um - up) * n;
1258 *   velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
1259 *  
1260 *  
1261 *   const auto mass_flux =
1262 *   um - gamma_m * gamma_sum_inv * (um - up) +
1263 *   gamma_m * gamma_p * gamma_sum_inv * (pm - pp) * n;
1264 *  
1265 *   pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
1266 *   }
1267 *   }
1268 *  
1269 * @endcode
1270 *
1271 * This function evaluates the inner face integrals.
1272 *
1273 * @code
1274 *   template <typename VectorType>
1275 *   void local_apply_face(
1276 *   const MatrixFree<dim, Number> &matrix_free,
1277 *   VectorType &dst,
1278 *   const VectorType &src,
1279 *   const std::pair<unsigned int, unsigned int> &face_range) const
1280 *   {
1281 *   FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
1282 *   matrix_free, true, 0, 0, 0);
1283 *   FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_p(
1284 *   matrix_free, false, 0, 0, 0);
1285 *   FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
1286 *   matrix_free, true, 0, 0, 1);
1287 *   FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_p(
1288 *   matrix_free, false, 0, 0, 1);
1289 *  
1290 *   MaterialEvaluation material(matrix_free, *material_data);
1291 *  
1292 *   for (unsigned int face = face_range.first; face < face_range.second;
1293 *   ++face)
1294 *   {
1295 *   velocity_m.reinit(face);
1296 *   velocity_p.reinit(face);
1297 *  
1298 *   pressure_m.reinit(face);
1299 *   pressure_p.reinit(face);
1300 *  
1301 *   pressure_m.gather_evaluate(src, EvaluationFlags::values);
1302 *   pressure_p.gather_evaluate(src, EvaluationFlags::values);
1303 *  
1304 *   velocity_m.gather_evaluate(src, EvaluationFlags::values);
1305 *   velocity_p.gather_evaluate(src, EvaluationFlags::values);
1306 *  
1307 *   material.reinit_face(face);
1308 *   evaluate_face_kernel<true>(pressure_m,
1309 *   velocity_m,
1310 *   pressure_p,
1311 *   velocity_p,
1312 *   material.get_speed_of_sound(),
1313 *   material.get_density());
1314 *  
1315 *   pressure_m.integrate_scatter(EvaluationFlags::values, dst);
1316 *   pressure_p.integrate_scatter(EvaluationFlags::values, dst);
1317 *   velocity_m.integrate_scatter(EvaluationFlags::values, dst);
1318 *   velocity_p.integrate_scatter(EvaluationFlags::values, dst);
1319 *   }
1320 *   }
1321 *  
1322 *  
1323 * @endcode
1324 *
1325 *
1326 * <a name="step_89-Matrixfreeboundaryfunctionforpointtopointinterpolation"></a>
1327 * <h4>Matrix-free boundary function for point-to-point interpolation</h4>
1328 *
1329
1330 *
1331 * The following function then evaluates the boundary face integrals and the
1332 * non-matching face integrals using point-to-point interpolation.
1333 *
1334 * @code
1335 *   template <typename VectorType>
1336 *   void local_apply_boundary_face_point_to_point(
1337 *   const MatrixFree<dim, Number> &matrix_free,
1338 *   VectorType &dst,
1339 *   const VectorType &src,
1340 *   const std::pair<unsigned int, unsigned int> &face_range) const
1341 *   {
1342 *   FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
1343 *   matrix_free, true, 0, 0, 0);
1344 *   FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
1345 *   matrix_free, true, 0, 0, 1);
1346 *  
1347 *   BCEvaluationP pressure_bc(pressure_m);
1348 *   BCEvaluationU velocity_bc(velocity_m);
1349 *  
1350 *   MaterialEvaluation material(matrix_free, *material_data);
1351 *  
1352 * @endcode
1353 *
1354 * Remote evaluators.
1355 *
1356 * @code
1357 *   auto pressure_r = pressure_r_eval->get_data_accessor();
1358 *   auto velocity_r = velocity_r_eval->get_data_accessor();
1359 *   auto c_r = c_r_eval->get_data_accessor();
1360 *   auto rho_r = rho_r_eval->get_data_accessor();
1361 *  
1362 *   for (unsigned int face = face_range.first; face < face_range.second;
1363 *   ++face)
1364 *   {
1365 *   velocity_m.reinit(face);
1366 *   pressure_m.reinit(face);
1367 *  
1368 *   pressure_m.gather_evaluate(src, EvaluationFlags::values);
1369 *   velocity_m.gather_evaluate(src, EvaluationFlags::values);
1370 *  
1371 *   if (HelperFunctions::is_non_matching_face(
1372 *   remote_face_ids, matrix_free.get_boundary_id(face)))
1373 *   {
1374 * @endcode
1375 *
1376 * If @c face is non-matching we have to query values via the
1377 * FERemoteEvaluaton objects. This is done by passing the
1378 * corresponding FERemoteEvaluaton objects to the function that
1379 * evaluates the kernel. As mentioned above, each side of the
1380 * non-matching interface is traversed separately and we do not
1381 * have to consider the neighbor in the kernel. Note, that the
1382 * values in the FERemoteEvaluaton objects are already updated at
1383 * this point.
1384 *
1385
1386 *
1387 * For point-to-point interpolation we simply use the
1388 * corresponding FERemoteEvaluaton objects in combination with the
1389 * standard FEFaceEvaluation objects.
1390 *
1391 * @code
1392 *   velocity_r.reinit(face);
1393 *   pressure_r.reinit(face);
1394 *  
1395 *   material.reinit_face(face);
1396 *  
1397 * @endcode
1398 *
1399 * If we are considering a homogeneous material, do not use the
1400 * inhomogeneous fluxes. While it would be possible
1401 * to use the inhomogeneous fluxes they are more expensive to
1402 * compute.
1403 *
1404 * @code
1405 *   if (material.is_homogeneous())
1406 *   {
1407 *   evaluate_face_kernel<false>(pressure_m,
1408 *   velocity_m,
1409 *   pressure_r,
1410 *   velocity_r,
1411 *   material.get_speed_of_sound(),
1412 *   material.get_density());
1413 *   }
1414 *   else
1415 *   {
1416 *   c_r.reinit(face);
1417 *   rho_r.reinit(face);
1418 *   evaluate_face_kernel_inhomogeneous(
1419 *   pressure_m,
1420 *   velocity_m,
1421 *   pressure_r,
1422 *   velocity_r,
1423 *   material.get_speed_of_sound(),
1424 *   material.get_density(),
1425 *   c_r,
1426 *   rho_r);
1427 *   }
1428 *   }
1429 *   else
1430 *   {
1431 * @endcode
1432 *
1433 * If @c face is a standard boundary face, evaluate the integral
1434 * as usual in the matrix free context. To be able to use the same
1435 * kernel as for inner faces we pass the boundary condition
1436 * objects to the function that evaluates the kernel. As detailed
1437 * above `weight_neighbor=false`.
1438 *
1439 * @code
1440 *   material.reinit_face(face);
1441 *   evaluate_face_kernel<false>(pressure_m,
1442 *   velocity_m,
1443 *   pressure_bc,
1444 *   velocity_bc,
1445 *   material.get_speed_of_sound(),
1446 *   material.get_density());
1447 *   }
1448 *  
1449 *   pressure_m.integrate_scatter(EvaluationFlags::values, dst);
1450 *   velocity_m.integrate_scatter(EvaluationFlags::values, dst);
1451 *   }
1452 *   }
1453 *  
1454 * @endcode
1455 *
1456 *
1457 * <a name="step_89-MatrixfreeboundaryfunctionforNitschetypemortaring"></a>
1458 * <h4>Matrix-free boundary function for Nitsche-type mortaring</h4>
1459 *
1460
1461 *
1462 * This function evaluates the boundary face integrals and the
1463 * non-matching face integrals using Nitsche-type mortaring.
1464 *
1465 * @code
1466 *   template <typename VectorType>
1467 *   void local_apply_boundary_face_mortaring(
1468 *   const MatrixFree<dim, Number> &matrix_free,
1469 *   VectorType &dst,
1470 *   const VectorType &src,
1471 *   const std::pair<unsigned int, unsigned int> &face_range) const
1472 *   {
1473 *   FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
1474 *   matrix_free, true, 0, 0, 0);
1475 *   FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
1476 *   matrix_free, true, 0, 0, 1);
1477 *  
1478 * @endcode
1479 *
1480 * For Nitsche-type mortaring we are evaluating the integrals over
1481 * intersections. This is why, quadrature points are arbitrarily
1482 * distributed on every face. Thus, we can not make use of face batches
1483 * and FEFaceEvaluation but have to consider each face individually and
1484 * make use of @c FEFacePointEvaluation to evaluate the integrals in the
1485 * arbitrarily distributed quadrature points.
1486 * Since the setup of FEFacePointEvaluation is more expensive than that of
1487 * FEEvaluation we do the setup only once. For this we are using the
1488 * helper function @c get_thread_safe_fe_face_point_evaluation_object().
1489 *
1490 * @code
1491 *   FEFacePointEvaluation<1, dim, dim, Number> &pressure_m_mortar =
1492 *   get_thread_safe_fe_face_point_evaluation_object<1>(
1493 *   thread_local_pressure_m_mortar, 0);
1494 *   FEFacePointEvaluation<dim, dim, dim, Number> &velocity_m_mortar =
1495 *   get_thread_safe_fe_face_point_evaluation_object<dim>(
1496 *   thread_local_velocity_m_mortar, 1);
1497 *  
1498 *   BCEvaluationP pressure_bc(pressure_m);
1499 *   BCEvaluationU velocity_bc(velocity_m);
1500 *  
1501 *   MaterialEvaluation material(matrix_free, *material_data);
1502 *  
1503 *   auto pressure_r_mortar = pressure_r_eval->get_data_accessor();
1504 *   auto velocity_r_mortar = velocity_r_eval->get_data_accessor();
1505 *   auto c_r = c_r_eval->get_data_accessor();
1506 *   auto rho_r = rho_r_eval->get_data_accessor();
1507 *  
1508 *   for (unsigned int face = face_range.first; face < face_range.second;
1509 *   ++face)
1510 *   {
1511 *   if (HelperFunctions::is_non_matching_face(
1512 *   remote_face_ids, matrix_free.get_boundary_id(face)))
1513 *   {
1514 *   material.reinit_face(face);
1515 *  
1516 * @endcode
1517 *
1518 * First fetch the DoF values with standard FEFaceEvaluation
1519 * objects.
1520 *
1521 * @code
1522 *   pressure_m.reinit(face);
1523 *   velocity_m.reinit(face);
1524 *  
1525 *   pressure_m.read_dof_values(src);
1526 *   velocity_m.read_dof_values(src);
1527 *  
1528 * @endcode
1529 *
1530 * Project the internally stored values into the face DoFs
1531 * of the current face.
1532 *
1533 * @code
1534 *   pressure_m.project_to_face(EvaluationFlags::values);
1535 *   velocity_m.project_to_face(EvaluationFlags::values);
1536 *  
1537 * @endcode
1538 *
1539 * For mortaring, we have to consider every face from the face
1540 * batches separately and have to use the FEFacePointEvaluation
1541 * objects to be able to evaluate the integrals with the
1542 * arbitrarily distributed quadrature points.
1543 *
1544 * @code
1545 *   for (unsigned int v = 0;
1546 *   v < matrix_free.n_active_entries_per_face_batch(face);
1547 *   ++v)
1548 *   {
1549 *   constexpr unsigned int n_lanes =
1550 *   VectorizedArray<Number>::size();
1551 *   velocity_m_mortar.reinit(face * n_lanes + v);
1552 *   pressure_m_mortar.reinit(face * n_lanes + v);
1553 *  
1554 * @endcode
1555 *
1556 * Evaluate using FEFacePointEvaluation. As buffer,
1557 * simply use the internal buffers from the
1558 * FEFaceEvaluation objects.
1559 *
1560 * @code
1561 *   velocity_m_mortar.evaluate_in_face(
1562 *   &velocity_m.get_scratch_data().begin()[0][v],
1563 *   EvaluationFlags::values);
1564 *  
1565 *   pressure_m_mortar.evaluate_in_face(
1566 *   &pressure_m.get_scratch_data().begin()[0][v],
1567 *   EvaluationFlags::values);
1568 *  
1569 *   velocity_r_mortar.reinit(face * n_lanes + v);
1570 *   pressure_r_mortar.reinit(face * n_lanes + v);
1571 *  
1572 * @endcode
1573 *
1574 * As above, if we are considering a homogeneous
1575 * material, do not use the inhomogeneous
1576 * fluxes. Since we are operating on face @c v we
1577 * call @c material.get_density()[v].
1578 *
1579 * @code
1580 *   if (material.is_homogeneous())
1581 *   {
1582 *   evaluate_face_kernel<false>(
1583 *   pressure_m_mortar,
1584 *   velocity_m_mortar,
1585 *   pressure_r_mortar,
1586 *   velocity_r_mortar,
1587 *   material.get_speed_of_sound()[v],
1588 *   material.get_density()[v]);
1589 *   }
1590 *   else
1591 *   {
1592 *   c_r.reinit(face * n_lanes + v);
1593 *   rho_r.reinit(face * n_lanes + v);
1594 *  
1595 *   evaluate_face_kernel_inhomogeneous(
1596 *   pressure_m_mortar,
1597 *   velocity_m_mortar,
1598 *   pressure_r_mortar,
1599 *   velocity_r_mortar,
1600 *   material.get_speed_of_sound()[v],
1601 *   material.get_density()[v],
1602 *   c_r,
1603 *   rho_r);
1604 *   }
1605 *  
1606 * @endcode
1607 *
1608 * Integrate using FEFacePointEvaluation. As buffer,
1609 * simply use the internal buffers from the
1610 * FEFaceEvaluation objects.
1611 *
1612 * @code
1613 *   velocity_m_mortar.integrate_in_face(
1614 *   &velocity_m.get_scratch_data().begin()[0][v],
1615 *   EvaluationFlags::values);
1616 *  
1617 *   pressure_m_mortar.integrate_in_face(
1618 *   &pressure_m.get_scratch_data().begin()[0][v],
1619 *   EvaluationFlags::values);
1620 *   }
1621 *  
1622 * @endcode
1623 *
1624 * Collect the contributions from the face DoFs to
1625 * the internal cell DoFs to be able to use the
1626 * member function @c distribute_local_to_global().
1627 *
1628 * @code
1629 *   pressure_m.collect_from_face(EvaluationFlags::values);
1630 *   velocity_m.collect_from_face(EvaluationFlags::values);
1631 *  
1632 *   pressure_m.distribute_local_to_global(dst);
1633 *   velocity_m.distribute_local_to_global(dst);
1634 *   }
1635 *   else
1636 *   {
1637 *   velocity_m.reinit(face);
1638 *   pressure_m.reinit(face);
1639 *  
1640 *   pressure_m.gather_evaluate(src, EvaluationFlags::values);
1641 *   velocity_m.gather_evaluate(src, EvaluationFlags::values);
1642 *  
1643 *   material.reinit_face(face);
1644 *   evaluate_face_kernel<false>(pressure_m,
1645 *   velocity_m,
1646 *   pressure_bc,
1647 *   velocity_bc,
1648 *   material.get_speed_of_sound(),
1649 *   material.get_density());
1650 *  
1651 *   pressure_m.integrate_scatter(EvaluationFlags::values, dst);
1652 *   velocity_m.integrate_scatter(EvaluationFlags::values, dst);
1653 *   }
1654 *   }
1655 *   }
1656 *  
1657 *   const MatrixFree<dim, Number> &matrix_free;
1658 *  
1659 *   const std::shared_ptr<CellwiseMaterialData<Number>> material_data;
1660 *  
1661 *   const std::set<types::boundary_id> remote_face_ids;
1662 *  
1663 * @endcode
1664 *
1665 * FERemoteEvaluation objects are stored as shared pointers. This way,
1666 * they can also be used for other operators without precomputing the values
1667 * multiple times.
1668 *
1669 * @code
1670 *   const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1671 *   pressure_r_eval;
1672 *   const std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
1673 *   velocity_r_eval;
1674 *  
1675 *   const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1676 *   c_r_eval;
1677 *   const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1678 *   rho_r_eval;
1679 *  
1680 *   const std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>>
1681 *   nm_mapping_info;
1682 *  
1683 * @endcode
1684 *
1685 * We store FEFacePointEvaluation objects as members in a thread local
1686 * way, since its creation is more expensive compared to FEEvaluation
1687 * objects.
1688 *
1689 * @code
1690 *   mutable Threads::ThreadLocalStorage<
1691 *   std::unique_ptr<FEFacePointEvaluation<1, dim, dim, Number>>>
1692 *   thread_local_pressure_m_mortar;
1693 *  
1694 *   mutable Threads::ThreadLocalStorage<
1695 *   std::unique_ptr<FEFacePointEvaluation<dim, dim, dim, Number>>>
1696 *   thread_local_velocity_m_mortar;
1697 *  
1698 * @endcode
1699 *
1700 * Helper function to create and get FEFacePointEvaluation objects in a
1701 * thread safe way. On each thread, FEFacePointEvaluation is created if it
1702 * has not been created by now. After that, simply return the object
1703 * corresponding to the thread under consideration.
1704 *
1705 * @code
1706 *   template <int n_components>
1707 *   FEFacePointEvaluation<n_components, dim, dim, Number> &
1708 *   get_thread_safe_fe_face_point_evaluation_object(
1709 *   Threads::ThreadLocalStorage<
1710 *   std::unique_ptr<FEFacePointEvaluation<n_components, dim, dim, Number>>>
1711 *   &fe_face_point_eval_thread_local,
1712 *   unsigned int fist_selected_comp) const
1713 *   {
1714 *   if (fe_face_point_eval_thread_local.get() == nullptr)
1715 *   {
1716 *   fe_face_point_eval_thread_local = std::make_unique<
1717 *   FEFacePointEvaluation<n_components, dim, dim, Number>>(
1718 *   *nm_mapping_info,
1719 *   matrix_free.get_dof_handler().get_fe(),
1720 *   true,
1721 *   fist_selected_comp);
1722 *   }
1723 *   return *fe_face_point_eval_thread_local.get();
1724 *   }
1725 *   };
1726 *  
1727 * @endcode
1728 *
1729 *
1730 * <a name="step_89-Inversemassoperator"></a>
1731 * <h3>Inverse mass operator</h3>
1732 *
1733
1734 *
1735 * For the time stepping methods below, we need the inverse mass
1736 * operator. We apply it via a loop over all (batches of) cells as
1737 * always, where the contribution of each cell is computed in a
1738 * matrix-free way:
1739 *
1740 * @code
1741 *   template <int dim, typename Number>
1742 *   class InverseMassOperator
1743 *   {
1744 *   public:
1745 *   InverseMassOperator(const MatrixFree<dim, Number> &matrix_free)
1746 *   : matrix_free(matrix_free)
1747 *   {}
1748 *  
1749 *   template <typename VectorType>
1750 *   void apply(VectorType &dst, const VectorType &src) const
1751 *   {
1752 *   dst.zero_out_ghost_values();
1753 *   matrix_free.cell_loop(&InverseMassOperator::local_apply_cell<VectorType>,
1754 *   this,
1755 *   dst,
1756 *   src);
1757 *   }
1758 *  
1759 *   private:
1760 *   template <typename VectorType>
1761 *   void local_apply_cell(
1762 *   const MatrixFree<dim, Number> &mf,
1763 *   VectorType &dst,
1764 *   const VectorType &src,
1765 *   const std::pair<unsigned int, unsigned int> &cell_range) const
1766 *   {
1767 *   FEEvaluation<dim, -1, 0, dim + 1, Number> phi(mf);
1768 *   MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1, dim + 1, Number>
1769 *   minv(phi);
1770 *  
1771 *   for (unsigned int cell = cell_range.first; cell < cell_range.second;
1772 *   ++cell)
1773 *   {
1774 *   phi.reinit(cell);
1775 *   phi.read_dof_values(src);
1776 *   minv.apply(phi.begin_dof_values(), phi.begin_dof_values());
1777 *   phi.set_dof_values(dst);
1778 *   }
1779 *   }
1780 *  
1781 *   const MatrixFree<dim, Number> &matrix_free;
1782 *   };
1783 *  
1784 * @endcode
1785 *
1786 *
1787 * <a name="step_89-RungeKuttatimestepping"></a>
1788 * <h3>Runge-Kutta time-stepping</h3>
1789 *
1790
1791 *
1792 * This class implements a Runge-Kutta scheme of order 2.
1793 *
1794 * @code
1795 *   template <int dim, typename Number, typename remote_value_type>
1796 *   class RungeKutta2
1797 *   {
1798 *   using VectorType = LinearAlgebra::distributed::Vector<Number>;
1799 *  
1800 *   public:
1801 *   RungeKutta2(
1802 *   const std::shared_ptr<InverseMassOperator<dim, Number>>
1803 *   inverse_mass_operator,
1804 *   const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
1805 *   acoustic_operator)
1806 *   : inverse_mass_operator(inverse_mass_operator)
1807 *   , acoustic_operator(acoustic_operator)
1808 *   {}
1809 *  
1810 * @endcode
1811 *
1812 * The `run()` function of this class contains the time loop. It
1813 * starts by initializing some member variables (such as
1814 * short-cuts to objects that describe the finite element, its
1815 * properties, and the mapping; an by initializing vectors). It
1816 * then computes the time step size via minimum element edge
1817 * length. Assuming non-distorted elements, we can compute the
1818 * edge length as the distance between two vertices. From this,
1819 * we can then obtain the time step size via the CFL condition.
1820 *
1821 * @code
1822 *   void run(const MatrixFree<dim, Number> &matrix_free,
1823 *   const double cr,
1824 *   const double end_time,
1825 *   const double speed_of_sound,
1826 *   const Function<dim> &initial_condition,
1827 *   const std::string &vtk_prefix)
1828 *   {
1829 *   const auto &dof_handler = matrix_free.get_dof_handler();
1830 *   const auto &mapping = *matrix_free.get_mapping_info().mapping;
1831 *   const auto degree = dof_handler.get_fe().degree;
1832 *  
1833 *   VectorType solution;
1834 *   matrix_free.initialize_dof_vector(solution);
1835 *   VectorType solution_temp;
1836 *   matrix_free.initialize_dof_vector(solution_temp);
1837 *  
1838 *   HelperFunctions::set_initial_condition(matrix_free,
1839 *   initial_condition,
1840 *   solution);
1841 *  
1842 *   double h_local_min = std::numeric_limits<double>::max();
1843 *   for (const auto &cell : dof_handler.active_cell_iterators())
1844 *   h_local_min =
1845 *   std::min(h_local_min,
1846 *   (cell->vertex(1) - cell->vertex(0)).norm_square());
1847 *   h_local_min = std::sqrt(h_local_min);
1848 *   const double h_min =
1849 *   Utilities::MPI::min(h_local_min, dof_handler.get_mpi_communicator());
1850 *  
1851 *   const double dt =
1852 *   cr * HelperFunctions::compute_dt_cfl(h_min, degree, speed_of_sound);
1853 *  
1854 * @endcode
1855 *
1856 * We can then perform the time integration loop:
1857 *
1858 * @code
1859 *   double time = 0.0;
1860 *   unsigned int timestep = 0;
1861 *   while (time < end_time)
1862 *   {
1863 *   HelperFunctions::write_vtu(solution,
1864 *   matrix_free.get_dof_handler(),
1865 *   mapping,
1866 *   degree,
1867 *   "step_89-" + vtk_prefix +
1868 *   std::to_string(timestep));
1869 *  
1870 *   std::swap(solution, solution_temp);
1871 *   time += dt;
1872 *   ++timestep;
1873 *   perform_time_step(dt, solution, solution_temp);
1874 *   }
1875 *   }
1876 *  
1877 * @endcode
1878 *
1879 * The main work of this class is done by a `private` member
1880 * function that performs one Runge-Kutta 2 time step. Recall that
1881 * this method requires two sub-steps ("stages") computing
1882 * intermediate values `k1` and `k2`, after which the intermediate
1883 * values are summed with weights to obtain the new solution at
1884 * the end of the time step. The RK2 method allows for the
1885 * elimination of the intermediate vector `k2` by utilizing the
1886 * `dst` vector as temporary storage.
1887 *
1888 * @code
1889 *   private:
1890 *   void
1891 *   perform_time_step(const double dt, VectorType &dst, const VectorType &src)
1892 *   {
1893 *   VectorType k1 = src;
1894 *  
1895 *   /* First stage. */
1896 *   evaluate_stage(k1, src);
1897 *  
1898 *   /* Second stage. */
1899 *   k1.sadd(0.5 * dt, 1.0, src);
1900 *   evaluate_stage(dst, k1);
1901 *  
1902 *   /* Summing things into the output vector. */
1903 *   dst.sadd(dt, 1.0, src);
1904 *   }
1905 *  
1906 * @endcode
1907 *
1908 * Evaluating a single Runge-Kutta stage is a straightforward step
1909 * that really only requires applying the operator once, and then
1910 * applying the inverse of the mass matrix.
1911 *
1912 * @code
1913 *   void evaluate_stage(VectorType &dst, const VectorType &src)
1914 *   {
1915 *   acoustic_operator->evaluate(dst, src);
1916 *   dst *= -1.0;
1917 *   inverse_mass_operator->apply(dst, dst);
1918 *   }
1919 *  
1920 *   const std::shared_ptr<InverseMassOperator<dim, Number>>
1921 *   inverse_mass_operator;
1922 *   const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
1923 *   acoustic_operator;
1924 *   };
1925 *  
1926 *  
1927 * @endcode
1928 *
1929 *
1930 * <a name="step_89-Constructionofnonmatchingtriangulations"></a>
1931 * <h3>Construction of non-matching triangulations</h3>
1932 *
1933
1934 *
1935 * Let us now make our way to the higher-level functions of this program.
1936 *
1937
1938 *
1939 * The first of these functions creates a two dimensional square
1940 * triangulation that spans from @f$(0,0)@f$ to @f$(1,1)@f$. It consists of
1941 * two sub-domains. The left sub-domain spans from @f$(0,0)@f$ to
1942 * @f$(0.525,1)@f$. The right sub-domain spans from @f$(0.525,0)@f$ to
1943 * @f$(1,1)@f$. The left sub-domain has elements that are three times
1944 * smaller compared to the ones for the right sub-domain.
1945 *
1946
1947 *
1948 * At non-matching interfaces, we need to provide different boundary
1949 * IDs for the cells that make up the two parts (because, while they
1950 * may be physically adjacent, they are not logically neighbors
1951 * given that the faces of cells on both sides do not match, and the
1952 * Triangulation class will therefore treat the interface between
1953 * the two parts as a "boundary"). These boundary IDs have to differ
1954 * because later on RemotePointEvaluation has to search for remote
1955 * points for each face, that are defined in the same mesh (since we
1956 * merge the mesh) but not on the same side of the non-matching
1957 * interface. As a consequence, we declare at the top symbolic names
1958 * for these boundary indicators, and ensure that we return a set
1959 * with these values to the caller for later use.
1960 *
1961
1962 *
1963 * The actual mesh is then constructed by first constructing the
1964 * left and right parts separately (setting material ids to zero and
1965 * one, respectively), and using the appropriate boundary ids for
1966 * all parts of the mesh. We then use
1967 * GridGenerator::merge_triangulations() to combine them into one
1968 * (non-matching) mesh. We have to pay attention that should the two
1969 * sub-triangulations have vertices at the same locations, that they
1970 * are not merged (connecting the two triangulations logically)
1971 * since we want the interface to be an actual boundary. We achieve
1972 * this by providing a tolerance of zero for the merge, see the
1973 * documentation of the function
1974 * GridGenerator::merge_triangulations().
1975 *
1976 * @code
1977 *   template <int dim>
1978 *   void build_non_matching_triangulation(
1979 *   Triangulation<dim> &tria,
1980 *   std::set<types::boundary_id> &non_matching_faces,
1981 *   const unsigned int refinements)
1982 *   {
1983 *   const double length = 1.0;
1984 *  
1985 *   const types::boundary_id non_matching_id_left = 98;
1986 *   const types::boundary_id non_matching_id_right = 99;
1987 *  
1988 *   non_matching_faces = {non_matching_id_left, non_matching_id_right};
1989 *  
1990 *   /* Construct left part of mesh. */
1991 *   Triangulation<dim> tria_left;
1992 *   const unsigned int subdiv_left = 11;
1993 *   GridGenerator::subdivided_hyper_rectangle(tria_left,
1994 *   {subdiv_left, 2 * subdiv_left},
1995 *   {0.0, 0.0},
1996 *   {0.525 * length, length});
1997 *  
1998 *   for (const auto &cell : tria_left.active_cell_iterators())
1999 *   cell->set_material_id(0);
2000 *   for (const auto &face : tria_left.active_face_iterators())
2001 *   if (face->at_boundary())
2002 *   {
2003 *   face->set_boundary_id(0);
2004 *   if (face->center()[0] > 0.525 * length - 1e-6)
2005 *   face->set_boundary_id(non_matching_id_left);
2006 *   }
2007 *  
2008 *   /* Construct right part of mesh. */
2009 *   Triangulation<dim> tria_right;
2010 *   const unsigned int subdiv_right = 4;
2011 *   GridGenerator::subdivided_hyper_rectangle(tria_right,
2012 *   {subdiv_right, 2 * subdiv_right},
2013 *   {0.525 * length, 0.0},
2014 *   {length, length});
2015 *  
2016 *   for (const auto &cell : tria_right.active_cell_iterators())
2017 *   cell->set_material_id(1);
2018 *   for (const auto &face : tria_right.active_face_iterators())
2019 *   if (face->at_boundary())
2020 *   {
2021 *   face->set_boundary_id(0);
2022 *   if (face->center()[0] < 0.525 * length + 1e-6)
2023 *   face->set_boundary_id(non_matching_id_right);
2024 *   }
2025 *  
2026 *   /* Merge triangulations. */
2027 *   GridGenerator::merge_triangulations(tria_left,
2028 *   tria_right,
2029 *   tria,
2030 *   /*tolerance*/ 0.,
2031 *   /*copy_manifold_ids*/ false,
2032 *   /*copy_boundary_ids*/ true);
2033 *  
2034 *   /* Refine the result. */
2035 *   tria.refine_global(refinements);
2036 *   }
2037 *  
2038 * @endcode
2039 *
2040 *
2041 * <a name="step_89-Setupandrunningofthetwoschemes"></a>
2042 * <h3>Set up and running of the two schemes</h3>
2043 *
2044
2045 *
2046 *
2047 * <a name="step_89-Setupandrunningofthepointtopointinterpolationscheme"></a>
2048 * <h4>Setup and running of the point-to-point interpolation scheme</h4>
2049 *
2050
2051 *
2052 * We are now at the two functions that run the overall schemes (the
2053 * point-to-point and the mortaring schemes). The first of these
2054 * functions fills a FERemoteEvaluationCommunicator object that is
2055 * needed for point-to-point interpolation. Additionally, the
2056 * corresponding remote evaluators are set up using this remote
2057 * communicator. Eventually, the operators are handed to the time
2058 * integrator that runs the simulation.
2059 *
2060 * @code
2061 *   template <int dim, typename Number>
2062 *   void run_with_point_to_point_interpolation(
2063 *   const MatrixFree<dim, Number> &matrix_free,
2064 *   const std::set<types::boundary_id> &non_matching_faces,
2065 *   const std::map<types::material_id, std::pair<double, double>> &materials,
2066 *   const double end_time,
2067 *   const Function<dim> &initial_condition,
2068 *   const std::string &vtk_prefix)
2069 *   {
2070 *   const auto &dof_handler = matrix_free.get_dof_handler();
2071 *   const auto &tria = dof_handler.get_triangulation();
2072 *  
2073 * @endcode
2074 *
2075 * Communication objects know about the communication pattern. That is,
2076 * they know about the cells and quadrature points that have to be
2077 * evaluated at remote faces. This information is given via
2078 * RemotePointEvaluation. Additionally, the communication objects
2079 * have to be able to match the quadrature points of the remote
2080 * points (that provide exterior information) to the quadrature points
2081 * defined at the interior cell. In case of point-to-point interpolation
2082 * a vector of pairs with face batch Ids and the number of faces in the
2083 * batch is needed. @c FERemoteCommunicationObjectEntityBatches
2084 * is a container to store this information.
2085 *
2086
2087 *
2088 * The information is filled outside of the actual class since in some cases
2089 * the information is available from some heuristic and
2090 * it is possible to skip some expensive operations. This is for example
2091 * the case for sliding rotating interfaces with equally spaced elements on
2092 * both sides of the non-matching interface @cite duerrwaechter2021an.
2093 *
2094
2095 *
2096 * For the standard case of point to point-to-point interpolation without
2097 * any heuristic we make use of the utility function
2098 * Utilities::compute_remote_communicator_faces_point_to_point_interpolation().
2099 * Please refer to the documentation of this function to see how to manually
2100 * set up the remote communicator from outside.
2101 *
2102
2103 *
2104 * Among the inputs for the remote communicator we need a list of
2105 * boundary ids for the non-matching faces, along with a function
2106 * object for each boundary id that returns a vector of true/false
2107 * flags in which exactly the vertices of cells of the
2108 * triangulation are marked that have a face at the boundary id in
2109 * question.
2110 *
2111 * @code
2112 *   std::vector<
2113 *   std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
2114 *   non_matching_faces_marked_vertices;
2115 *   for (const auto &nm_face : non_matching_faces)
2116 *   {
2117 *   auto marked_vertices = [&nm_face, &tria]() -> std::vector<bool> {
2118 *   std::vector<bool> mask(tria.n_vertices(), true);
2119 *  
2120 *   for (const auto &cell : tria.active_cell_iterators())
2121 *   for (auto const &f : cell->face_indices())
2122 *   if (cell->face(f)->at_boundary() &&
2123 *   cell->face(f)->boundary_id() == nm_face)
2124 *   for (const auto v : cell->vertex_indices())
2125 *   mask[cell->vertex_index(v)] = false;
2126 *  
2127 *   return mask;
2128 *   };
2129 *  
2130 *   non_matching_faces_marked_vertices.emplace_back(
2131 *   std::make_pair(nm_face, marked_vertices));
2132 *   }
2133 *  
2134 *   auto remote_communicator =
2135 *   Utilities::compute_remote_communicator_faces_point_to_point_interpolation(
2136 *   matrix_free, non_matching_faces_marked_vertices);
2137 *  
2138 * @endcode
2139 *
2140 * We are using point-to-point interpolation and can therefore
2141 * easily access the corresponding data at face batches. This
2142 * is why we use a @c VectorizedArray as @c remote_value_type:
2143 *
2144 * @code
2145 *   using remote_value_type = VectorizedArray<Number>;
2146 *  
2147 * @endcode
2148 *
2149 * We then set up FERemoteEvaluation objects that access the
2150 * pressure and velocity at remote faces, along with an object to
2151 * describe cell-wise material data.
2152 *
2153 * @code
2154 *   const auto pressure_r =
2155 *   std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2156 *   remote_communicator, dof_handler, /*first_selected_component*/ 0);
2157 *  
2158 *   const auto velocity_r =
2159 *   std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
2160 *   remote_communicator, dof_handler, /*first_selected_component*/ 1);
2161 *  
2162 *   const auto material_data =
2163 *   std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
2164 *  
2165 * @endcode
2166 *
2167 * If we have an inhomogeneous problem, we have to set up the
2168 * material handler that accesses the materials at remote faces.
2169 *
2170 * @code
2171 *   const auto c_r =
2172 *   std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2173 *   remote_communicator,
2174 *   matrix_free.get_dof_handler().get_triangulation(),
2175 *   /*first_selected_component*/ 0);
2176 *   const auto rho_r =
2177 *   std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2178 *   remote_communicator,
2179 *   matrix_free.get_dof_handler().get_triangulation(),
2180 *   /*first_selected_component*/ 0);
2181 *  
2182 * @endcode
2183 *
2184 * If the domain is not homogeneous, i.e., if material parameters
2185 * change from cell to cell, we initialize and fill DoF vectors
2186 * that contain the material properties. Materials do not change
2187 * during the simulation, therefore there is no need to ever
2188 * compute the values after the first @c gather_evaluate() (at the
2189 * end of the following block) again.
2190 *
2191 * @code
2192 *   if (!material_data->is_homogeneous())
2193 *   {
2194 *   Vector<Number> c(
2195 *   matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2196 *   Vector<Number> rho(
2197 *   matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2198 *  
2199 *   for (const auto &cell : matrix_free.get_dof_handler()
2200 *   .get_triangulation()
2201 *   .active_cell_iterators())
2202 *   {
2203 *   c[cell->active_cell_index()] =
2204 *   materials.at(cell->material_id()).first;
2205 *   rho[cell->active_cell_index()] =
2206 *   materials.at(cell->material_id()).second;
2207 *   }
2208 *  
2209 *   c_r->gather_evaluate(c, EvaluationFlags::values);
2210 *   rho_r->gather_evaluate(rho, EvaluationFlags::values);
2211 *   }
2212 *  
2213 *  
2214 * @endcode
2215 *
2216 * Next, we set up the inverse mass operator and the acoustic
2217 * operator. Using `remote_value_type=VectorizedArray<Number>`
2218 * makes the operator use point-to-point interpolation. These two
2219 * objects are then used to create a `RungeKutta2` object to
2220 * perform the time integration.
2221 *
2222
2223 *
2224 * We also compute the maximum speed of sound, needed for the
2225 * computation of the time-step size, and then run the time integrator.
2226 * For the examples considered here, we found a limiting Courant number of
2227 * @f$\mathrm{Cr}\approx 0.36@f$ to maintain stability. To ensure, the
2228 * error of the temporal discretization is small, we use a considerably
2229 * smaller Courant number of @f$0.2@f$.
2230 *
2231 * @code
2232 *   const auto inverse_mass_operator =
2233 *   std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
2234 *  
2235 *   const auto acoustic_operator =
2236 *   std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
2237 *   matrix_free,
2238 *   material_data,
2239 *   non_matching_faces,
2240 *   pressure_r,
2241 *   velocity_r,
2242 *   c_r,
2243 *   rho_r);
2244 *  
2245 *   RungeKutta2<dim, Number, remote_value_type> time_integrator(
2246 *   inverse_mass_operator, acoustic_operator);
2247 *  
2248 *   double speed_of_sound_max = 0.0;
2249 *   for (const auto &mat : materials)
2250 *   speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
2251 *  
2252 *  
2253 *   time_integrator.run(matrix_free,
2254 *   /*Cr*/ 0.2,
2255 *   end_time,
2256 *   speed_of_sound_max,
2257 *   initial_condition,
2258 *   vtk_prefix);
2259 *   }
2260 *  
2261 * @endcode
2262 *
2263 *
2264 * <a name="step_89-SetupandrunningoftheNitschetypemortaringscheme"></a>
2265 * <h4>Setup and running of the Nitsche-type mortaring scheme</h4>
2266 *
2267
2268 *
2269 * The alternative to the previous function is to use the mortaring
2270 * scheme -- implemented in the following function. This function
2271 * can only be run when deal.II is configured using CGAL (but the
2272 * previous function can still be used without CGAL), and so errors
2273 * out if CGAL is not available.
2274 *
2275 * @code
2276 *   template <int dim, typename Number>
2277 *   void run_with_nitsche_type_mortaring(
2278 *   const MatrixFree<dim, Number> &matrix_free,
2279 *   const std::set<types::boundary_id> &non_matching_faces,
2280 *   const std::map<types::material_id, std::pair<double, double>> &materials,
2281 *   const double end_time,
2282 *   const Function<dim> &initial_condition,
2283 *   const std::string &vtk_prefix)
2284 *   {
2285 *   #ifndef DEAL_II_WITH_CGAL
2286 *   (void)matrix_free;
2287 *   (void)non_matching_faces;
2288 *   (void)materials;
2289 *   (void)end_time;
2290 *   (void)initial_condition;
2291 *   (void)vtk_prefix;
2292 *  
2293 *   ConditionalOStream pcout(
2294 *   std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
2295 *  
2296 *   pcout << "In this function, mortars are computed using CGAL. "
2297 *   "Configure deal.II with DEAL_II_WITH_CGAL to run this function.\n";
2298 *  
2299 *   return;
2300 *   #else
2301 *  
2302 *   const auto &dof_handler = matrix_free.get_dof_handler();
2303 *   const auto &tria = dof_handler.get_triangulation();
2304 *   const auto &mapping = *matrix_free.get_mapping_info().mapping;
2305 *   const auto n_quadrature_pnts = matrix_free.get_quadrature().size();
2306 *  
2307 *   std::vector<
2308 *   std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
2309 *   non_matching_faces_marked_vertices;
2310 *  
2311 *   for (const auto &nm_face : non_matching_faces)
2312 *   {
2313 *   auto marked_vertices = [&]() {
2314 *   std::vector<bool> mask(tria.n_vertices(), true);
2315 *  
2316 *   for (const auto &cell : tria.active_cell_iterators())
2317 *   for (auto const &f : cell->face_indices())
2318 *   if (cell->face(f)->at_boundary() &&
2319 *   cell->face(f)->boundary_id() == nm_face)
2320 *   for (const auto v : cell->vertex_indices())
2321 *   mask[cell->vertex_index(v)] = false;
2322 *  
2323 *   return mask;
2324 *   };
2325 *  
2326 *   non_matching_faces_marked_vertices.emplace_back(
2327 *   std::make_pair(nm_face, marked_vertices));
2328 *   }
2329 *  
2330 * @endcode
2331 *
2332 * The only parts in this function that are functionally different
2333 * from the previous one follow here.
2334 *
2335
2336 *
2337 * First, quadrature points are arbitrarily distributed on each non-matching
2338 * face. Therefore, we have to make use of FEFacePointEvaluation.
2339 * FEFacePointEvaluation needs NonMatching::MappingInfo to work at the
2340 * correct quadrature points that are in sync with used FERemoteEvaluation
2341 * object. Using
2342 * `compute_remote_communicator_faces_nitsche_type_mortaring()` to reinit
2343 * NonMatching::MappingInfo ensures this. In the case of mortaring, we have
2344 * to use the weights provided by the quadrature rules that are used to set
2345 * up NonMatching::MappingInfo. Therefore we set the flag @c
2346 * use_global_weights.
2347 *
2348 * @code
2349 *   typename NonMatching::MappingInfo<dim, dim, Number>::AdditionalData
2350 *   additional_data;
2351 *   additional_data.use_global_weights = true;
2352 *  
2353 * @endcode
2354 *
2355 * Set up NonMatching::MappingInfo with needed update flags and
2356 * @c additional_data.
2357 *
2358 * @code
2359 *   auto nm_mapping_info =
2360 *   std::make_shared<NonMatching::MappingInfo<dim, dim, Number>>(
2361 *   mapping,
2362 *   update_values | update_JxW_values | update_normal_vectors |
2363 *   update_quadrature_points,
2364 *   additional_data);
2365 *  
2366 *   auto remote_communicator =
2367 *   Utilities::compute_remote_communicator_faces_nitsche_type_mortaring(
2368 *   matrix_free,
2369 *   non_matching_faces_marked_vertices,
2370 *   n_quadrature_pnts,
2371 *   0,
2372 *   nm_mapping_info.get());
2373 *  
2374 * @endcode
2375 *
2376 * Second, since quadrature points are arbitrarily distributed we
2377 * have to consider each face in a batch separately and can not
2378 * make use of @c VecorizedArray.
2379 *
2380 * @code
2381 *   using remote_value_type = Number;
2382 *  
2383 * @endcode
2384 *
2385 * The rest of the code is then identical to what we had in the
2386 * previous function (though it functions differently because of
2387 * the difference in `remote_value_type`).
2388 *
2389 * @code
2390 *   const auto pressure_r =
2391 *   std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2392 *   remote_communicator, dof_handler, /*first_selected_component*/ 0);
2393 *  
2394 *   const auto velocity_r =
2395 *   std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
2396 *   remote_communicator, dof_handler, /*first_selected_component*/ 1);
2397 *  
2398 *   const auto material_data =
2399 *   std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
2400 *  
2401 *   const auto c_r =
2402 *   std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2403 *   remote_communicator,
2404 *   matrix_free.get_dof_handler().get_triangulation(),
2405 *   /*first_selected_component*/ 0);
2406 *   const auto rho_r =
2407 *   std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2408 *   remote_communicator,
2409 *   matrix_free.get_dof_handler().get_triangulation(),
2410 *   /*first_selected_component*/ 0);
2411 *  
2412 *   if (!material_data->is_homogeneous())
2413 *   {
2414 *   Vector<Number> c(
2415 *   matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2416 *   Vector<Number> rho(
2417 *   matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2418 *  
2419 *   for (const auto &cell : matrix_free.get_dof_handler()
2420 *   .get_triangulation()
2421 *   .active_cell_iterators())
2422 *   {
2423 *   c[cell->active_cell_index()] =
2424 *   materials.at(cell->material_id()).first;
2425 *   rho[cell->active_cell_index()] =
2426 *   materials.at(cell->material_id()).second;
2427 *   }
2428 *  
2429 *   c_r->gather_evaluate(c, EvaluationFlags::values);
2430 *   rho_r->gather_evaluate(rho, EvaluationFlags::values);
2431 *   }
2432 *  
2433 *   const auto inverse_mass_operator =
2434 *   std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
2435 *  
2436 *   const auto acoustic_operator =
2437 *   std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
2438 *   matrix_free,
2439 *   material_data,
2440 *   non_matching_faces,
2441 *   pressure_r,
2442 *   velocity_r,
2443 *   c_r,
2444 *   rho_r,
2445 *   nm_mapping_info);
2446 *  
2447 *   RungeKutta2<dim, Number, remote_value_type> time_integrator(
2448 *   inverse_mass_operator, acoustic_operator);
2449 *  
2450 *   double speed_of_sound_max = 0.0;
2451 *   for (const auto &mat : materials)
2452 *   speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
2453 *  
2454 *   time_integrator.run(matrix_free,
2455 *   /*Cr*/ 0.2,
2456 *   end_time,
2457 *   speed_of_sound_max,
2458 *   initial_condition,
2459 *   vtk_prefix);
2460 *   #endif
2461 *   }
2462 *   } // namespace Step89
2463 *  
2464 *  
2465 * @endcode
2466 *
2467 *
2468 * <a name="step_89-main"></a>
2469 * <h3>main()</h3>
2470 *
2471
2472 *
2473 * Finally, the `main()` function executes the different versions of handling
2474 * non-matching interfaces.
2475 *
2476
2477 *
2478 * Similar to @ref step_87 "step-87", the minimum requirement of this tutorial is MPI.
2479 * The parallel::distributed::Triangulation class is used if deal.II is
2480 * configured with p4est. Otherwise parallel::shared::Triangulation is used.
2481 *
2482 * @code
2483 *   int main(int argc, char *argv[])
2484 *   {
2485 *   using namespace dealii;
2486 *   constexpr int dim = 2;
2487 *   using Number = double;
2488 *  
2489 *   Utilities::MPI::MPI_InitFinalize mpi(argc, argv);
2490 *   std::cout.precision(5);
2491 *   ConditionalOStream pcout(std::cout,
2492 *   (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
2493 *   0));
2494 *  
2495 *   const unsigned int refinements = 1;
2496 *   const unsigned int degree = 3;
2497 *  
2498 *   #ifdef DEAL_II_WITH_P4EST
2499 *   parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
2500 *   #else
2501 *   parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
2502 *   #endif
2503 *  
2504 *   pcout << "Create non-matching grid..." << std::endl;
2505 *  
2506 *   std::set<types::boundary_id> non_matching_faces;
2507 *   Step89::build_non_matching_triangulation(tria,
2508 *   non_matching_faces,
2509 *   refinements);
2510 *  
2511 *   pcout << " - Refinement level: " << refinements << std::endl;
2512 *   pcout << " - Number of cells: " << tria.n_cells() << std::endl;
2513 *  
2514 *   pcout << "Create DoFHandler..." << std::endl;
2515 *   DoFHandler<dim> dof_handler(tria);
2516 *   dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(degree) ^ (dim + 1)));
2517 *   pcout << " - Number of DoFs: " << dof_handler.n_dofs() << std::endl;
2518 *  
2519 *   AffineConstraints<Number> constraints;
2520 *   constraints.close();
2521 *  
2522 *   pcout << "Set up MatrixFree..." << std::endl;
2523 *   typename MatrixFree<dim, Number>::AdditionalData data;
2524 *   data.mapping_update_flags = update_gradients | update_values;
2525 *   data.mapping_update_flags_inner_faces = update_values;
2526 *   data.mapping_update_flags_boundary_faces =
2527 *   update_quadrature_points | update_values;
2528 *  
2529 *   MatrixFree<dim, Number> matrix_free;
2530 *   matrix_free.reinit(
2531 *   MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(degree + 1), data);
2532 *  
2533 *  
2534 * @endcode
2535 *
2536 *
2537 * <a name="step_89-RunvibratingmembranetestcaseHomogeneouspressure"></a>
2538 * <h4>Run vibrating membrane test case} Homogeneous pressure</h4>
2539 * Dirichlet boundary conditions are applied for
2540 * simplicity. Therefore, modes can not be chosen arbitrarily.
2541 *
2542 * @code
2543 *   pcout << "Run vibrating membrane test case..." << std::endl;
2544 *   const double modes = 10.0;
2545 *   std::map<types::material_id, std::pair<double, double>> homogeneous_material;
2546 *   homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0);
2547 *   const auto initial_solution_membrane =
2548 *   Step89::InitialConditionVibratingMembrane<dim>(modes);
2549 *  
2550 *   /* Run vibrating membrane test case using point-to-point interpolation: */
2551 *   pcout << " - Point-to-point interpolation: " << std::endl;
2552 *   Step89::run_with_point_to_point_interpolation(
2553 *   matrix_free,
2554 *   non_matching_faces,
2555 *   homogeneous_material,
2556 *   8.0 * initial_solution_membrane.get_period_duration(
2557 *   homogeneous_material.begin()->second.first),
2558 *   initial_solution_membrane,
2559 *   "vm-p2p");
2560 *  
2561 *   /* Run vibrating membrane test case using Nitsche-type mortaring: */
2562 *   pcout << " - Nitsche-type mortaring: " << std::endl;
2563 *   Step89::run_with_nitsche_type_mortaring(
2564 *   matrix_free,
2565 *   non_matching_faces,
2566 *   homogeneous_material,
2567 *   8.0 * initial_solution_membrane.get_period_duration(
2568 *   homogeneous_material.begin()->second.first),
2569 *   initial_solution_membrane,
2570 *   "vm-nitsche");
2571 *  
2572 * @endcode
2573 *
2574 *
2575 * <a name="step_89-Runtestcasewithinhomogeneousmaterial"></a>
2576 * <h4>Run test case with inhomogeneous material</h4>
2577 * Run simple test case with inhomogeneous material and Nitsche-type
2578 * mortaring:
2579 *
2580 * @code
2581 *   pcout << "Run test case with inhomogeneous material..." << std::endl;
2582 *   std::map<types::material_id, std::pair<double, double>>
2583 *   inhomogeneous_material;
2584 *   inhomogeneous_material[0] = std::make_pair(1.0, 1.0);
2585 *   inhomogeneous_material[1] = std::make_pair(3.0, 1.0);
2586 *   Step89::run_with_nitsche_type_mortaring(matrix_free,
2587 *   non_matching_faces,
2588 *   inhomogeneous_material,
2589 *   /*runtime*/ 0.3,
2590 *   Step89::GaussPulse<dim>(0.3, 0.5),
2591 *   "inhomogeneous");
2592 *  
2593 *  
2594 *   return 0;
2595 *   }
2596 * @endcode
2597<a name="step_89-Results"></a><h1>Results</h1>
2598
2599
2600<a name="step_89-VibratingmembranePointtopointinterpolationvsNitschetypemortaring"></a><h3>Vibrating membrane: Point-to-point interpolation vs. Nitsche-type mortaring</h3>
2601
2602
2603We compare the results of the simulations after the last time step, i.e. at @f$t=8T@f$.
2604The @f$y@f$-component of the velocity field using Nitsche-type mortaring is depicted on the left.
2605The same field using point-to-point interpolation is depicted on the right.
2606
2607<table align="center" class="doxtable">
2608 <tr>
2609 <td>
2610 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_velocity_Y.png "" width=60%
2611 </td>
2612 <td>
2613 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_point_to_point_velocity_Y.png "" width=60%
2614 </td>
2615 </tr>
2616</table>
2617
2618Besides this, the results for the pressure and the velocity in @f$y@f$ direction
2619are plotted along the horizontal line that spans from (0,0.69) to (1,0.69).
2620
2621<table align="center" class="doxtable">
2622 <tr>
2623 <td>
2624 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_pressure.svg "" width=150%
2625 </td>
2626 <td>
2627 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_velocity_Y.svg "" width=150%
2628 </td>
2629 </tr>
2630</table>
2631
2632While the results of the pressure are similar, @f$u_y@f$ clearly differs. At certain
2633positions we can see aliasing errors for the point-to-point interpolation.
2634For different mesh configurations and/or longer run times, the aliasing effects
2635of the point-to-point simulation accumulate and the simulation becomes instable.
2636To keep the tutorial short we have chosen one mesh that can be used for all
2637examples. For a configuration that yields instable results for a wide range of
2638polynomial degrees, see @cite heinz2022high.
2639
2640<a name="step_89-Wavepropagationthroughinhomogeneousfluid"></a><h3>Wave propagation through in-homogeneous fluid</h3>
2641
2642
2643The example that follows is just one example in which non-matching discretizations can be efficiently
2644used to reduce the number of DoFs. The example is nice, since results for a similar
2645test case are shown in multiple publications. As before, we slightly adapted the
2646test case to be able to use the same mesh for all simulations. The pressure field
2647at @f$t=0.3@f$ is depicted below.
2648
2649@image html https://www.dealii.org/images/steps/developer/step_89_inhomogenous_test_case_pressure.png "" width=30%
2650
2651As expected, we can easily see that the wave length in the right domain is roughly
2652three times times the wave length in the left domain. Hence, the wave can be
2653resolved with a coarser discretization.
2654
2655Using the same element size in the whole domain, we can compute a reference result.
2656The displayed reference result is obtained by choosing the same subdivision level
2657for both sub-domains, i.e. @c subdiv_right = 11. In this particular example the
2658reference result uses @f$92,928@f$ DoFs, while the non-matching result uses @f$52,608@f$ DoFs.
2659The pressure result is plotted along the horizontal line that spans from (0,0.5) to (1,0.5).
2660
2661@image html https://www.dealii.org/images/steps/developer/step_89_inhomogenous_test_case_conforming_vs_nonmatching.svg "" width=60%
2662
2663The results computed using the non-matching discretization are clearly in good agreement with
2664the reference result.
2665
2666<a name="step_89-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2667
2668
2669All the implementations are done with overlapping triangulations in mind. In particular the
2670intersections in the mortaring case are constructed such that they are computed correctly
2671for overlapping triangulations. For this the intersection requests are of dimension @f$dim-1@f$.
2672The cells which are intersected with the intersection requests are of dimension @f$dim@f$. For the
2673simple case of non-conforming interfaces it would be sufficient to compute the intersections
2674between @f$dim-1@f$ and @f$dim-1@f$ entities. Furthermore, the lambda could be adapted, such that cells are
2675only marked if they are connected to a certain boundary ID (in this case, e.g., 99) instead of
2676marking every cell that is <i>not</i> connected to the opposite boundary ID (in this case, e.g., 98).
2677Marking fewer cells can reduce the setup costs significantly.
2678
2679Note that the use of inhomogeneous materials in this procedure is questionable, since it is not clear which
2680material is present in the overlapping part of the mesh.
2681 *
2682 *
2683<a name="step_89-PlainProg"></a>
2684<h1> The plain program</h1>
2685@include "step-89.cc"
2686*/
value_type get_value(const unsigned int q_point) const
T read_cell_data(const AlignedVector< T > &array) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Abstract base class for mapping classes.
Definition mapping.h:318
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
Definition point.h:111
const value_type get_value(const unsigned int q) const
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
const Event initial
Definition event.cc:64
void write_vtu(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
@ matrix
Contents is actually a matrix.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static constexpr double PI
Definition numbers.h:254
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned int material_id
Definition types.h:167