fe_system.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_fe_system_h
#  define dealii_fe_system_h


/*----------------------------   fe_system.h     ---------------------------*/


#  include <deal.II/base/config.h>

#  include <deal.II/base/thread_management.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

template <int dim, int spacedim>
class FE_Enriched;


template <int dim, int spacedim = dim>
class FESystem : public FiniteElement<dim, spacedim>
{
public:
  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 = typename enable_if_all<
      (std::is_same<typename std::decay<FEPairs>::type,
                    std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
                              unsigned int>>::value ||
       std::is_base_of<FiniteElement<dim, spacedim>,
                       typename std::decay<FEPairs>::type>::value)...>::type>
  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> &&) = default; // NOLINT

  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 bool         face_orientation = true,
                     const bool         face_flip        = false,
                     const bool         face_rotation = false) 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 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 override;


  virtual void
  get_subface_interpolation_matrix(const FiniteElement<dim, spacedim> &source,
                                   const unsigned int                  subface,
                                   FullMatrix<double> &matrix) 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 override;

  virtual FiniteElementDomination::Domination
  compare_for_face_domination(
    const FiniteElement<dim, spacedim> &fe_other) const override;

  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;

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;

  virtual std::unique_ptr<
    typename FiniteElement<dim, spacedim>::InternalDataBase>
  get_face_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 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;

  virtual void
  fill_fe_face_values(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
    const unsigned int                                          face_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;

  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 <int dim_1>
  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 Quadrature<dim_1> &                                   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<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;

  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;
  };

  /*
   * Mutex for protecting initialization of restriction and embedding matrix.
   */
  mutable Threads::Mutex 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::unique_ptr<FiniteElement<dim, spacedim>>,
                            unsigned int>(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
/*----------------------------  fe_system.h  ---------------------------*/