quadrature_generator.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2020 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_non_matching_quadrature_generator_h
#define dealii_non_matching_quadrature_generator_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_restriction.h>
#include <deal.II/base/quadrature.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/grid/tria.h>
 
#include <deal.II/hp/q_collection.h>
 
#include <deal.II/non_matching/immersed_surface_quadrature.h>
 
#include <functional>
#include <optional>
 
DEAL_II_NAMESPACE_OPEN
namespace NonMatching
{
  namespace internal
  {
    namespace QuadratureGeneratorImplementation
    {
      template <int dim, int spacedim>
      class QGenerator;
    } // namespace QuadratureGeneratorImplementation
 
 
    namespace DiscreteQuadratureGeneratorImplementation
    {
      template <int dim>
      class CellWiseFunction;
    } // namespace DiscreteQuadratureGeneratorImplementation
  }   // namespace internal
 
 
  struct AdditionalQGeneratorData
  {
    AdditionalQGeneratorData(const unsigned int max_box_splits          = 4,
                             const double lower_bound_implicit_function = 1e-11,
                             const double min_distance_between_roots    = 1e-12,
                             const double limit_to_be_definite          = 1e-11,
                             const double root_finder_tolerance         = 1e-12,
                             const unsigned int max_root_finder_splits  = 2,
                             bool               split_in_half           = true);
 
    unsigned int max_box_splits;
 
    double lower_bound_implicit_function;
 
    double min_distance_between_roots;
 
    double limit_to_be_definite;
 
    double root_finder_tolerance;
 
    unsigned int max_root_finder_splits;
 
    bool split_in_half;
  };
 
 
 
  template <int dim>
  class QuadratureGenerator
  {
  public:
    using AdditionalData = AdditionalQGeneratorData;
 
    QuadratureGenerator(
      const hp::QCollection<1> &quadratures1D,
      const AdditionalData     &additional_data = AdditionalData());
 
    void
    clear_quadratures();
 
    void
    generate(const Function<dim> &level_set, const BoundingBox<dim> &box);
 
    void
    generate_append(const Function<dim>    &level_set,
                    const BoundingBox<dim> &box);
 
    const Quadrature<dim> &
    get_inside_quadrature() const;
 
    const Quadrature<dim> &
    get_outside_quadrature() const;
 
    const ImmersedSurfaceQuadrature<dim> &
    get_surface_quadrature() const;
 
    void
    set_1D_quadrature(const unsigned int q_index);
 
  private:
    internal::QuadratureGeneratorImplementation::QGenerator<dim, dim>
      q_generator;
  };
 
 
  template <int dim>
  class FaceQuadratureGenerator
  {
  public:
    using AdditionalData = AdditionalQGeneratorData;
 
    FaceQuadratureGenerator(
      const hp::QCollection<1> &quadratures1D,
      const AdditionalData     &additional_data = AdditionalData());
 
    void
    clear_quadratures();
 
    void
    generate(const Function<dim>    &level_set,
             const BoundingBox<dim> &box,
             const unsigned int      face_index);
 
    void
    generate_append(const Function<dim>    &level_set,
                    const BoundingBox<dim> &box,
                    const unsigned int      face_index);
 
    const Quadrature<dim - 1> &
    get_inside_quadrature() const;
 
    const Quadrature<dim - 1> &
    get_outside_quadrature() const;
 
    const ImmersedSurfaceQuadrature<dim - 1, dim> &
    get_surface_quadrature() const;
 
    void
    set_1D_quadrature(const unsigned int q_index);
 
  private:
    QuadratureGenerator<dim - 1> quadrature_generator;
 
    ImmersedSurfaceQuadrature<dim - 1, dim> surface_quadrature;
  };
 
 
  template <>
  class FaceQuadratureGenerator<1>
  {
  public:
    using AdditionalData = AdditionalQGeneratorData;
 
    FaceQuadratureGenerator(
      const hp::QCollection<1> &quadratures1D,
      const AdditionalData     &additional_data = AdditionalData());
 
    void
    clear_quadratures();
 
    void
    generate(const Function<1>    &level_set,
             const BoundingBox<1> &box,
             const unsigned int    face_index);
 
    void
    generate_append(const Function<1>    &level_set,
                    const BoundingBox<1> &box,
                    const unsigned int    face_index);
 
    const Quadrature<0> &
    get_inside_quadrature() const;
 
    const Quadrature<0> &
    get_outside_quadrature() const;
 
    const ImmersedSurfaceQuadrature<0, 1> &
    get_surface_quadrature() const;
 
    void
    set_1D_quadrature(const unsigned int q_index);
 
  private:
    Quadrature<0> inside_quadrature;
 
    Quadrature<0> outside_quadrature;
 
    const ImmersedSurfaceQuadrature<0, 1> surface_quadrature;
  };
 
 
  template <int dim>
  class DiscreteQuadratureGenerator : public QuadratureGenerator<dim>
  {
  public:
    using AdditionalData = AdditionalQGeneratorData;
 
    template <typename Number>
    DiscreteQuadratureGenerator(
      const hp::QCollection<1> &quadratures1D,
      const DoFHandler<dim>    &dof_handler,
      const ReadVector<Number> &level_set,
      const AdditionalData     &additional_data = AdditionalData());
 
    void
    generate(const typename Triangulation<dim>::active_cell_iterator &cell);
 
  private:
    void
    generate_fe_q_iso_q1(const BoundingBox<dim> &unit_box);
 
    std::unique_ptr<internal::DiscreteQuadratureGeneratorImplementation::
                      CellWiseFunction<dim>>
      reference_space_level_set;
  };
 
  template <int dim>
  class DiscreteFaceQuadratureGenerator : public FaceQuadratureGenerator<dim>
  {
  public:
    using AdditionalData = AdditionalQGeneratorData;
 
    template <typename Number>
    DiscreteFaceQuadratureGenerator(
      const hp::QCollection<1> &quadratures1D,
      const DoFHandler<dim>    &dof_handler,
      const ReadVector<Number> &level_set,
      const AdditionalData     &additional_data = AdditionalData());
 
    void
    generate(const typename Triangulation<dim>::active_cell_iterator &cell,
             const unsigned int face_index);
 
  private:
    void
    generate_fe_q_iso_q1(const BoundingBox<dim> &unit_box,
                         unsigned int            face_index);
 
    std::unique_ptr<internal::DiscreteQuadratureGeneratorImplementation::
                      CellWiseFunction<dim>>
      reference_space_level_set;
  };
 
 
  namespace internal
  {
    namespace QuadratureGeneratorImplementation
    {
      class RootFinder
      {
      public:
        struct AdditionalData
        {
          AdditionalData(const double       tolerance           = 1e-12,
                         const unsigned int max_recursion_depth = 2,
                         const unsigned int max_iterations      = 500);
 
          double tolerance;
 
          unsigned int max_recursion_depth;
 
          unsigned int max_iterations;
        };
 
 
        RootFinder(const AdditionalData &data = AdditionalData());
 
        void
        find_roots(const std::vector<std::reference_wrapper<const Function<1>>>
                                        &functions,
                   const BoundingBox<1> &interval,
                   std::vector<double>  &roots);
 
      private:
        void
        find_roots(const Function<1>    &function,
                   const BoundingBox<1> &interval,
                   const unsigned int    recursion_depth,
                   std::vector<double>  &roots);
 
        const AdditionalData additional_data;
      };
 
 
      template <int dim>
      class ExtendableQuadrature : public Quadrature<dim>
      {
      public:
        ExtendableQuadrature() = default;
 
        ExtendableQuadrature(const Quadrature<dim> &quadrature);
 
        void
        clear();
 
        void
        push_back(const Point<dim> &point, const double weight);
      };
 
 
      enum class Definiteness
      {
        negative,
        positive,
        indefinite
      };
 
 
      template <int dim>
      class QPartitioning
      {
      public:
        ExtendableQuadrature<dim> &
        quadrature_by_definiteness(const Definiteness definiteness);
 
        void
        clear();
 
        ExtendableQuadrature<dim> negative;
 
        ExtendableQuadrature<dim> positive;
 
        ExtendableQuadrature<dim> indefinite;
 
        ImmersedSurfaceQuadrature<dim> surface;
      };
 
 
      template <int dim, int spacedim>
      class UpThroughDimensionCreator
      {
      public:
        UpThroughDimensionCreator(
          const hp::QCollection<1>       &q_collection1D,
          const AdditionalQGeneratorData &additional_data);
 
        void
        generate(const std::vector<std::reference_wrapper<const Function<dim>>>
                                           &level_sets,
                 const BoundingBox<dim>    &box,
                 const Quadrature<dim - 1> &low_dim_quadrature,
                 const unsigned int         height_function_direction,
                 QPartitioning<dim>        &q_partitioning);
 
        void
        set_1D_quadrature(const unsigned int q_index);
 
      private:
        void
        create_surface_point(
          const Point<dim - 1> &point,
          const double          weight,
          const std::vector<std::reference_wrapper<const Function<dim>>>
                                         &level_sets,
          const BoundingBox<dim>         &box,
          const unsigned int              height_function_direction,
          ImmersedSurfaceQuadrature<dim> &surface_quadrature);
 
        const ObserverPointer<const hp::QCollection<1>> q_collection1D;
 
        const AdditionalQGeneratorData additional_data;
 
        unsigned int q_index;
 
        std::vector<Functions::PointRestriction<dim - 1>> point_restrictions;
 
        RootFinder root_finder;
 
        std::vector<double> roots;
      };
 
 
      struct HeightDirectionData
      {
        HeightDirectionData();
 
 
        unsigned int direction;
 
        double min_abs_dfdx;
      };
 
 
      template <int dim, int spacedim>
      class QGeneratorBase
      {
      public:
        QGeneratorBase(const hp::QCollection<1>       &q_collection1D,
                       const AdditionalQGeneratorData &additional_data);
 
        void
        clear_quadratures();
 
        const QPartitioning<dim> &
        get_quadratures() const;
 
      protected:
        const AdditionalQGeneratorData additional_data;
 
        unsigned int q_index;
 
        const ObserverPointer<const hp::QCollection<1>> q_collection1D;
 
        QPartitioning<dim> q_partitioning;
      };
 
 
      template <int dim, int spacedim>
      class QGenerator : public QGeneratorBase<dim, spacedim>
      {
      public:
        QGenerator(const hp::QCollection<1>       &q_collection1D,
                   const AdditionalQGeneratorData &additional_data);
 
        void
        generate(const std::vector<std::reference_wrapper<const Function<dim>>>
                                        &level_sets,
                 const BoundingBox<dim> &box,
                 const unsigned int      n_box_splits);
 
        void
        set_1D_quadrature(const unsigned int q_index);
 
      private:
        void
        create_low_dim_quadratures(
          const unsigned int height_function_direction,
          const std::vector<std::reference_wrapper<const Function<dim>>>
                                 &level_sets,
          const BoundingBox<dim> &box,
          const unsigned int      n_box_splits);
 
        void
        create_high_dim_quadratures(
          const unsigned int height_function_direction,
          const std::vector<std::reference_wrapper<const Function<dim>>>
                                 &level_sets,
          const BoundingBox<dim> &box);
 
        void
        split_box_and_recurse(
          const std::vector<std::reference_wrapper<const Function<dim>>>
                                                   &level_sets,
          const BoundingBox<dim>                   &box,
          const std::optional<HeightDirectionData> &direction_data,
          const unsigned int                        n_box_splits);
 
        void
        use_midpoint_method(
          const std::vector<std::reference_wrapper<const Function<dim>>>
                                 &level_sets,
          const BoundingBox<dim> &box);
 
        QGenerator<dim - 1, spacedim> low_dim_algorithm;
 
        UpThroughDimensionCreator<dim, spacedim> up_through_dimension_creator;
 
        hp::QCollection<dim> tensor_products;
      };
 
 
      template <int spacedim>
      class QGenerator<1, spacedim> : public QGeneratorBase<1, spacedim>
      {
      public:
        QGenerator(const hp::QCollection<1>       &quadratures1D,
                   const AdditionalQGeneratorData &additional_data);
 
        void
        generate(const std::vector<std::reference_wrapper<const Function<1>>>
                                      &level_sets,
                 const BoundingBox<1> &box,
                 const unsigned int    n_box_splits);
 
        void
        set_1D_quadrature(const unsigned int q_index);
 
      private:
        void
        create_surface_points(
          const std::vector<std::reference_wrapper<const Function<1>>>
            &level_sets);
 
        RootFinder root_finder;
 
        std::vector<double> roots;
 
        const unsigned int direction = 0;
 
        const Point<0> zero_dim_point;
        const double   unit_weight = 1;
      };
 
 
      template <int dim>
      void
      tensor_point_with_1D_quadrature(const Point<dim - 1> &point,
                                      const double          weight,
                                      const Quadrature<1>  &quadrature1D,
                                      const double          start,
                                      const double          end,
                                      const unsigned int    component_in_dim,
                                      ExtendableQuadrature<dim> &quadrature);
 
 
      template <int dim>
      Definiteness
      pointwise_definiteness(
        const std::vector<std::reference_wrapper<const Function<dim>>>
                         &functions,
        const Point<dim> &point);
 
 
      template <int dim>
      struct FunctionBounds
      {
      public:
        std::pair<double, double> value;
 
        std::array<std::pair<double, double>, dim> gradient;
      };
 
 
      template <int dim>
      std::pair<double, double>
      find_extreme_values(
        const std::vector<FunctionBounds<dim>> &all_function_bounds);
 
 
      template <int dim>
      std::optional<HeightDirectionData>
      find_best_height_direction(
        const std::vector<FunctionBounds<dim>> &all_function_bounds);
 
    } // namespace QuadratureGeneratorImplementation
 
 
    namespace DiscreteQuadratureGeneratorImplementation
    {
      template <int dim>
      class CellWiseFunction : public Function<dim>
      {
      public:
        virtual ~CellWiseFunction() = default;
 
        virtual void
        set_active_cell(
          const typename Triangulation<dim>::active_cell_iterator &cell) = 0;
 
        virtual void
        set_subcell(const std::vector<unsigned int> &mask,
                    const BoundingBox<dim>          &subcell_box) = 0;
 
        virtual bool
        is_fe_q_iso_q1() const = 0;
 
        virtual unsigned int
        n_subdivisions() const = 0;
      };
 
    } // namespace DiscreteQuadratureGeneratorImplementation
  }   // namespace internal
 
} // namespace NonMatching
DEAL_II_NAMESPACE_CLOSE
 
#endif /* dealii_non_matching_quadrature_generator_h */