fe_values.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2021 - 2024 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_fe_values
#define dealii_non_matching_fe_values
 
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/observer_pointer.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe_interface_values.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values.h>
 
#include <deal.II/grid/tria_iterator.h>
 
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/q_collection.h>
 
#include <deal.II/non_matching/fe_immersed_values.h>
#include <deal.II/non_matching/mesh_classifier.h>
#include <deal.II/non_matching/quadrature_generator.h>
 
#include <deque>
#include <optional>
 
DEAL_II_NAMESPACE_OPEN
 
namespace NonMatching
{
  struct RegionUpdateFlags
  {
    RegionUpdateFlags();
 
    UpdateFlags inside;
 
    UpdateFlags outside;
 
    UpdateFlags surface;
  };
 
 
  template <int dim>
  class FEValues
  {
  public:
    using AdditionalData = typename QuadratureGenerator<dim>::AdditionalData;
 
    template <typename Number>
    FEValues(const hp::FECollection<dim> &fe_collection,
             const Quadrature<1>         &quadrature,
             const RegionUpdateFlags      region_update_flags,
             const MeshClassifier<dim>   &mesh_classifier,
             const DoFHandler<dim>       &dof_handler,
             const ReadVector<Number>    &level_set,
             const AdditionalData        &additional_data = AdditionalData());
 
    template <typename Number>
    FEValues(const hp::MappingCollection<dim> &mapping_collection,
             const hp::FECollection<dim>      &fe_collection,
             const hp::QCollection<dim>       &q_collection,
             const hp::QCollection<1>         &q_collection_1d,
             const RegionUpdateFlags           region_update_flags,
             const MeshClassifier<dim>        &mesh_classifier,
             const DoFHandler<dim>            &dof_handler,
             const ReadVector<Number>         &level_set,
             const AdditionalData &additional_data = AdditionalData());
 
    template <bool level_dof_access>
    void
    reinit(
      const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
      const unsigned int q_index       = numbers::invalid_unsigned_int,
      const unsigned int mapping_index = numbers::invalid_unsigned_int);
 
    void
    reinit(const TriaIterator<CellAccessor<dim, dim>> &cell,
           const unsigned int q_index       = numbers::invalid_unsigned_int,
           const unsigned int mapping_index = numbers::invalid_unsigned_int,
           const unsigned int fe_index      = numbers::invalid_unsigned_int);
 
    const std::optional<::FEValues<dim>> &
    get_inside_fe_values() const;
 
    const std::optional<::FEValues<dim>> &
    get_outside_fe_values() const;
 
    const std::optional<FEImmersedSurfaceValues<dim>> &
    get_surface_fe_values() const;
 
  private:
    template <typename CellIteratorType>
    void
    reinit_internal(const CellIteratorType &cell,
                    const unsigned int      q_index,
                    const unsigned int      mapping_index,
                    const unsigned int      fe_index);
 
    void
    initialize(const hp::QCollection<dim> &q_collection);
 
    const ObserverPointer<const hp::MappingCollection<dim>> mapping_collection;
 
    const ObserverPointer<const hp::FECollection<dim>> fe_collection;
 
    const hp::QCollection<1> q_collection_1D;
 
    LocationToLevelSet current_cell_location;
 
    unsigned int active_fe_index;
 
    const RegionUpdateFlags region_update_flags;
 
    const ObserverPointer<const MeshClassifier<dim>> mesh_classifier;
 
    std::deque<std::optional<::FEValues<dim>>>
      fe_values_inside_full_quadrature;
 
    std::deque<std::optional<::FEValues<dim>>>
      fe_values_outside_full_quadrature;
 
    std::optional<::FEValues<dim>> fe_values_inside;
 
    std::optional<::FEValues<dim>> fe_values_outside;
 
    std::optional<NonMatching::FEImmersedSurfaceValues<dim>> fe_values_surface;
 
    DiscreteQuadratureGenerator<dim> quadrature_generator;
  };
 
 
 
  template <int dim>
  class FEInterfaceValues
  {
  public:
    using AdditionalData =
      typename FaceQuadratureGenerator<dim>::AdditionalData;
 
    template <typename Number>
    FEInterfaceValues(const hp::FECollection<dim> &fe_collection,
                      const Quadrature<1>         &quadrature,
                      const RegionUpdateFlags      region_update_flags,
                      const MeshClassifier<dim>   &mesh_classifier,
                      const DoFHandler<dim>       &dof_handler,
                      const ReadVector<Number>    &level_set,
                      const AdditionalData &additional_data = AdditionalData());
 
    template <typename Number>
    FEInterfaceValues(const hp::MappingCollection<dim> &mapping_collection,
                      const hp::FECollection<dim>      &fe_collection,
                      const hp::QCollection<dim - 1>   &q_collection,
                      const hp::QCollection<1>         &q_collection_1d,
                      const RegionUpdateFlags           region_update_flags,
                      const MeshClassifier<dim>        &mesh_classifier,
                      const DoFHandler<dim>            &dof_handler,
                      const ReadVector<Number>         &level_set,
                      const AdditionalData &additional_data = AdditionalData());
 
    template <typename CellIteratorType, typename CellNeighborIteratorType>
    void
    reinit(const CellIteratorType         &cell,
           const unsigned int              face_no,
           const unsigned int              sub_face_no,
           const CellNeighborIteratorType &cell_neighbor,
           const unsigned int              face_no_neighbor,
           const unsigned int              sub_face_no_neighbor,
           const unsigned int q_index       = numbers::invalid_unsigned_int,
           const unsigned int mapping_index = numbers::invalid_unsigned_int,
           const unsigned int fe_index      = numbers::invalid_unsigned_int);
 
 
    template <typename CellIteratorType>
    void
    reinit(const CellIteratorType &cell,
           const unsigned int      face_no,
           const unsigned int      q_index  = numbers::invalid_unsigned_int,
           const unsigned int mapping_index = numbers::invalid_unsigned_int,
           const unsigned int fe_index      = numbers::invalid_unsigned_int);
 
    const std::optional<::FEInterfaceValues<dim>> &
    get_inside_fe_values() const;
 
    const std::optional<::FEInterfaceValues<dim>> &
    get_outside_fe_values() const;
 
  private:
    void
    initialize(const hp::QCollection<dim - 1> &q_collection);
 
    template <typename CellAccessorType>
    void
    do_reinit(const TriaIterator<CellAccessorType>          &cell,
              const unsigned int                             face_no,
              const unsigned int                             q_index,
              const unsigned int                             active_fe_index,
              const std::function<void(::FEInterfaceValues<dim> &,
                                       const unsigned int)> &call_reinit);
 
    const ObserverPointer<const hp::MappingCollection<dim>> mapping_collection;
 
    const ObserverPointer<const hp::FECollection<dim>> fe_collection;
 
    const hp::QCollection<1> q_collection_1D;
 
    LocationToLevelSet current_face_location;
 
    const RegionUpdateFlags region_update_flags;
 
    const ObserverPointer<const MeshClassifier<dim>> mesh_classifier;
 
    std::optional<::FEInterfaceValues<dim>>
      fe_values_inside_full_quadrature;
 
    std::optional<::FEInterfaceValues<dim>>
      fe_values_outside_full_quadrature;
 
    std::optional<::FEInterfaceValues<dim>> fe_values_inside;
 
    std::optional<::FEInterfaceValues<dim>> fe_values_outside;
 
    DiscreteFaceQuadratureGenerator<dim> face_quadrature_generator;
  };
 
 
#ifndef DOXYGEN
 
  /*---------------------- Inline functions ---------------------*/
 
  template <int dim>
  template <typename CellIteratorType>
  inline void
  FEInterfaceValues<dim>::reinit(const CellIteratorType &cell,
                                 const unsigned int      face_no,
                                 const unsigned int      q_index,
                                 const unsigned int      mapping_index,
                                 const unsigned int      active_fe_index)
  {
    // Lambda describing how we should call reinit on a single
    // ::FEInterfaceValues object.
    const auto reinit_operation =
      [&](::FEInterfaceValues<dim> &fe_interface_values,
          const unsigned int              q_index) {
        fe_interface_values.reinit(
          cell, face_no, q_index, mapping_index, active_fe_index);
      };
 
    do_reinit(cell, face_no, q_index, active_fe_index, reinit_operation);
  }
 
 
 
  template <int dim>
  template <typename CellIteratorType, typename CellNeighborIteratorType>
  inline void
  FEInterfaceValues<dim>::reinit(const CellIteratorType         &cell,
                                 const unsigned int              face_no,
                                 const unsigned int              sub_face_no,
                                 const CellNeighborIteratorType &cell_neighbor,
                                 const unsigned int face_no_neighbor,
                                 const unsigned int sub_face_no_neighbor,
                                 const unsigned int q_index,
                                 const unsigned int mapping_index,
                                 const unsigned int active_fe_index)
  {
    Assert(sub_face_no == numbers::invalid_unsigned_int, ExcNotImplemented());
    Assert(sub_face_no_neighbor == numbers::invalid_unsigned_int,
           ExcNotImplemented());
 
    // Lambda describing how we should call reinit on a single
    // ::FEInterfaceValues object.
    const auto reinit_operation =
      [&](::FEInterfaceValues<dim> &fe_interface_values,
          const unsigned int              q_index) {
        fe_interface_values.reinit(cell,
                                   face_no,
                                   sub_face_no,
                                   cell_neighbor,
                                   face_no_neighbor,
                                   sub_face_no_neighbor,
                                   q_index,
                                   mapping_index,
                                   active_fe_index);
      };
 
    do_reinit(cell, face_no, q_index, active_fe_index, reinit_operation);
  }
 
#endif
 
} // namespace NonMatching
DEAL_II_NAMESPACE_CLOSE
 
#endif