fe_system.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 1999 - 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_fe_system_h
#define dealii_fe_system_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_tools.h>
 
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declaration
#ifndef DOXYGEN
template <int dim, int spacedim>
class FE_Enriched;
#endif
 
template <int dim, int spacedim = dim>
class FESystem : public FiniteElement<dim, spacedim>
{
public:
  FESystem() = delete;
 
  FESystem(const FiniteElement<dim, spacedim> &fe,
           const unsigned int                  n_elements);
 
  FESystem(const FiniteElement<dim, spacedim> &fe1,
           const unsigned int                  n1,
           const FiniteElement<dim, spacedim> &fe2,
           const unsigned int                  n2);
 
  FESystem(const FiniteElement<dim, spacedim> &fe1,
           const unsigned int                  n1,
           const FiniteElement<dim, spacedim> &fe2,
           const unsigned int                  n2,
           const FiniteElement<dim, spacedim> &fe3,
           const unsigned int                  n3);
 
  FESystem(const FiniteElement<dim, spacedim> &fe1,
           const unsigned int                  n1,
           const FiniteElement<dim, spacedim> &fe2,
           const unsigned int                  n2,
           const FiniteElement<dim, spacedim> &fe3,
           const unsigned int                  n3,
           const FiniteElement<dim, spacedim> &fe4,
           const unsigned int                  n4);
 
  FESystem(const FiniteElement<dim, spacedim> &fe1,
           const unsigned int                  n1,
           const FiniteElement<dim, spacedim> &fe2,
           const unsigned int                  n2,
           const FiniteElement<dim, spacedim> &fe3,
           const unsigned int                  n3,
           const FiniteElement<dim, spacedim> &fe4,
           const unsigned int                  n4,
           const FiniteElement<dim, spacedim> &fe5,
           const unsigned int                  n5);
 
  FESystem(const std::vector<const FiniteElement<dim, spacedim> *> &fes,
           const std::vector<unsigned int> &multiplicities);
 
#if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
  template <
    class... FEPairs,
    typename = enable_if_all_t<
      (std::is_same_v<std::decay_t<FEPairs>,
                      std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
                                unsigned int>> ||
       std::is_base_of_v<FiniteElement<dim, spacedim>,
                         std::decay_t<FEPairs>>)...>>
  FESystem(FEPairs &&...fe_pairs);
 
  FESystem(
    const std::initializer_list<
      std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
      &fe_systems);
#endif
 
  FESystem(const FESystem<dim, spacedim> &) = delete;
 
  FESystem(FESystem<dim, spacedim> &&other_fe_system) noexcept
    : FiniteElement<dim, spacedim>(std::move(other_fe_system))
    , base_elements(std::move(other_fe_system.base_elements))
    , generalized_support_points_index_table(
        std::move(other_fe_system.generalized_support_points_index_table))
  {}
 
  virtual ~FESystem() override = default;
 
  virtual std::string
  get_name() const override;
 
  virtual std::unique_ptr<FiniteElement<dim, spacedim>>
  clone() const override;
 
  virtual UpdateFlags
  requires_update_flags(const UpdateFlags update_flags) const override;
 
  // make variant with ComponentMask also available:
  using FiniteElement<dim, spacedim>::get_sub_fe;
 
  virtual const FiniteElement<dim, spacedim> &
  get_sub_fe(const unsigned int first_component,
             const unsigned int n_selected_components) const override;
 
  virtual double
  shape_value(const unsigned int i, const Point<dim> &p) const override;
 
  virtual double
  shape_value_component(const unsigned int i,
                        const Point<dim>  &p,
                        const unsigned int component) const override;
 
  virtual Tensor<1, dim>
  shape_grad(const unsigned int i, const Point<dim> &p) const override;
 
  virtual Tensor<1, dim>
  shape_grad_component(const unsigned int i,
                       const Point<dim>  &p,
                       const unsigned int component) const override;
 
  virtual Tensor<2, dim>
  shape_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
  virtual Tensor<2, dim>
  shape_grad_grad_component(const unsigned int i,
                            const Point<dim>  &p,
                            const unsigned int component) const override;
 
  virtual Tensor<3, dim>
  shape_3rd_derivative(const unsigned int i,
                       const Point<dim>  &p) const override;
 
  virtual Tensor<3, dim>
  shape_3rd_derivative_component(const unsigned int i,
                                 const Point<dim>  &p,
                                 const unsigned int component) const override;
 
  virtual Tensor<4, dim>
  shape_4th_derivative(const unsigned int i,
                       const Point<dim>  &p) const override;
 
  virtual Tensor<4, dim>
  shape_4th_derivative_component(const unsigned int i,
                                 const Point<dim>  &p,
                                 const unsigned int component) const override;
 
  virtual void
  get_interpolation_matrix(const FiniteElement<dim, spacedim> &source,
                           FullMatrix<double> &matrix) const override;
 
  virtual const FiniteElement<dim, spacedim> &
  base_element(const unsigned int index) const override;
 
  virtual bool
  has_support_on_face(const unsigned int shape_index,
                      const unsigned int face_index) const override;
 
  virtual const FullMatrix<double> &
  get_restriction_matrix(
    const unsigned int         child,
    const RefinementCase<dim> &refinement_case =
      RefinementCase<dim>::isotropic_refinement) const override;
 
  virtual const FullMatrix<double> &
  get_prolongation_matrix(
    const unsigned int         child,
    const RefinementCase<dim> &refinement_case =
      RefinementCase<dim>::isotropic_refinement) const override;
 
  virtual unsigned int
  face_to_cell_index(const unsigned int                 face_dof_index,
                     const unsigned int                 face,
                     const types::geometric_orientation combined_orientation =
                       numbers::default_geometric_orientation) const override;
 
  virtual Point<dim>
  unit_support_point(const unsigned int index) const override;
 
  virtual Point<dim - 1>
  unit_face_support_point(const unsigned int index,
                          const unsigned int face_no = 0) const override;
 
  virtual std::pair<Table<2, bool>, std::vector<unsigned int>>
  get_constant_modes() const override;
 
  virtual bool
  hp_constraints_are_implemented() const override;
 
  virtual void
  get_face_interpolation_matrix(const FiniteElement<dim, spacedim> &source,
                                FullMatrix<double>                 &matrix,
                                const unsigned int face_no = 0) const override;
 
 
  virtual void
  get_subface_interpolation_matrix(
    const FiniteElement<dim, spacedim> &source,
    const unsigned int                  subface,
    FullMatrix<double>                 &matrix,
    const unsigned int                  face_no = 0) const override;
 
  virtual std::vector<std::pair<unsigned int, unsigned int>>
  hp_vertex_dof_identities(
    const FiniteElement<dim, spacedim> &fe_other) const override;
 
  virtual std::vector<std::pair<unsigned int, unsigned int>>
  hp_line_dof_identities(
    const FiniteElement<dim, spacedim> &fe_other) const override;
 
  virtual std::vector<std::pair<unsigned int, unsigned int>>
  hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
                         const unsigned int face_no = 0) const override;
 
  virtual FiniteElementDomination::Domination
  compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
                         const unsigned int codim = 0) const override final;
 
  virtual void
  convert_generalized_support_point_values_to_dof_values(
    const std::vector<Vector<double>> &support_point_values,
    std::vector<double>               &dof_values) const override;
 
  virtual std::size_t
  memory_consumption() const override;
 
  struct BaseOffsets
  {
    unsigned int n_nonzero_components;
 
    unsigned int in_index;
 
    unsigned int out_index;
  };
 
protected:
  virtual std::unique_ptr<
    typename FiniteElement<dim, spacedim>::InternalDataBase>
  get_data(
    const UpdateFlags             update_flags,
    const Mapping<dim, spacedim> &mapping,
    const Quadrature<dim>        &quadrature,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  using FiniteElement<dim, spacedim>::get_face_data;
 
  virtual std::unique_ptr<
    typename FiniteElement<dim, spacedim>::InternalDataBase>
  get_face_data(
    const UpdateFlags               update_flags,
    const Mapping<dim, spacedim>   &mapping,
    const hp::QCollection<dim - 1> &quadrature,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  virtual std::unique_ptr<
    typename FiniteElement<dim, spacedim>::InternalDataBase>
  get_subface_data(
    const UpdateFlags             update_flags,
    const Mapping<dim, spacedim> &mapping,
    const Quadrature<dim - 1>    &quadrature,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  virtual void
  fill_fe_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const CellSimilarity::Similarity                            cell_similarity,
    const Quadrature<dim>                                      &quadrature,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                                  &mapping_data,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  using FiniteElement<dim, spacedim>::fill_fe_face_values;
 
  virtual void
  fill_fe_face_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const unsigned int                                          face_no,
    const hp::QCollection<dim - 1>                             &quadrature,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                                  &mapping_data,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  virtual void
  fill_fe_subface_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const unsigned int                                          face_no,
    const unsigned int                                          sub_no,
    const Quadrature<dim - 1>                                  &quadrature,
    const Mapping<dim, spacedim>                               &mapping,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
                                                                  &mapping_data,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
    ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
                                                                       spacedim>
      &output_data) const override;
 
  template <class Q_or_QC>
  void
  compute_fill(
    const Mapping<dim, spacedim>                               &mapping,
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const unsigned int                                          face_no,
    const unsigned int                                          sub_no,
    const Q_or_QC                                              &quadrature,
    const CellSimilarity::Similarity                            cell_similarity,
    const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
    const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_data,
    const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
      &mapping_data,
    internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
      &output_data) const;
 
private:
  static const unsigned int invalid_face_number = numbers::invalid_unsigned_int;
 
  std::vector<Table<2, unsigned int>> primitive_offset_tables;
 
  std::vector<std::vector<BaseOffsets>> nonprimitive_offset_tables;
 
  std::vector<std::pair<std::unique_ptr<const FiniteElement<dim, spacedim>>,
                        unsigned int>>
    base_elements;
 
  std::vector<std::vector<std::size_t>> generalized_support_points_index_table;
 
  void
  initialize(const std::vector<const FiniteElement<dim, spacedim> *> &fes,
             const std::vector<unsigned int> &multiplicities);
 
  void
  build_interface_constraints();
 
  template <int structdim>
  std::vector<std::pair<unsigned int, unsigned int>>
  hp_object_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
                           const unsigned int face_no = 0) const;
 
  class InternalData : public FiniteElement<dim, spacedim>::InternalDataBase
  {
  public:
    InternalData(const unsigned int n_base_elements);
 
    ~InternalData() override;
 
    void
    set_fe_data(
      const unsigned int base_no,
      std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>);
 
    typename FiniteElement<dim, spacedim>::InternalDataBase &
    get_fe_data(const unsigned int base_no) const;
 
    internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim> &
    get_fe_output_object(const unsigned int base_no) const;
 
  private:
    typename std::vector<
      std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
      base_fe_datas;
 
    mutable std::vector<
      internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>>
      base_fe_output_objects;
  };
 
  mutable std::mutex restriction_matrix_mutex;
  mutable std::mutex prolongation_matrix_mutex;
 
  friend class FE_Enriched<dim, spacedim>;
};
 
//------------------------variadic template constructor------------------------
 
#ifndef DOXYGEN
namespace internal
{
  namespace FESystemImplementation
  {
    template <int dim, int spacedim>
    unsigned int
    count_nonzeros(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &fe_systems)
    {
      return std::count_if(
        fe_systems.begin(),
        fe_systems.end(),
        [](const std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
                           unsigned int> &fe_system) {
          return fe_system.second > 0;
        });
    }
 
 
 
    template <int dim, int spacedim>
    std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>
    promote_to_fe_pair(const FiniteElement<dim, spacedim> &fe)
    {
      return std::make_pair(std::move(fe.clone()), 1u);
    }
 
 
 
    template <int dim, int spacedim>
    auto
    promote_to_fe_pair(std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
                                 unsigned int> &&p)
      -> decltype(std::forward<
                  std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
                            unsigned int>>(p))
    {
      return std::forward<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>(
        p);
    }
  } // namespace FESystemImplementation
} // namespace internal
 
 
 
#  if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
// We are just forwarding/delegating to the constructor taking a
// std::initializer_list. If we decide to remove the deprecated constructors, we
// might just use the variadic constructor with a suitable static_assert instead
// of the std::enable_if.
template <int dim, int spacedim>
template <class... FEPairs, typename>
FESystem<dim, spacedim>::FESystem(FEPairs &&...fe_pairs)
  : FESystem<dim, spacedim>(
      {internal::FESystemImplementation::promote_to_fe_pair<dim, spacedim>(
        std::forward<FEPairs>(fe_pairs))...})
{}
 
 
 
template <int dim, int spacedim>
FESystem<dim, spacedim>::FESystem(
  const std::initializer_list<
    std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
    &fe_systems)
  : FiniteElement<dim, spacedim>(
      FETools::Compositing::multiply_dof_numbers<dim, spacedim>(fe_systems),
      FETools::Compositing::compute_restriction_is_additive_flags<dim,
                                                                  spacedim>(
        fe_systems),
      FETools::Compositing::compute_nonzero_components<dim, spacedim>(
        fe_systems))
  , base_elements(internal::FESystemImplementation::count_nonzeros(fe_systems))
{
  std::vector<const FiniteElement<dim, spacedim> *> fes;
  std::vector<unsigned int>                         multiplicities;
 
  const auto extract =
    [&fes, &multiplicities](
      const std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
                      unsigned int> &fe_system) {
      fes.push_back(fe_system.first.get());
      multiplicities.push_back(fe_system.second);
    };
 
  for (const auto &p : fe_systems)
    extract(p);
 
  initialize(fes, multiplicities);
}
#  endif
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif