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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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 <vector>
#include <memory>
#include <utility>
#include <type_traits>


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

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

  FESystem (const FESystem<dim,spacedim> &) = delete;

  FESystem (FESystem<dim,spacedim> &&) = default; // NOLINT

  virtual ~FESystem () = default;

  virtual std::string get_name () const;

  virtual
  std::unique_ptr<FiniteElement<dim,spacedim> >
  clone() const;

  virtual
  UpdateFlags
  requires_update_flags (const UpdateFlags update_flags) const;

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

  virtual double shape_value (const unsigned int i,
                              const Point<dim> &p) const;

  virtual double shape_value_component (const unsigned int i,
                                        const Point<dim> &p,
                                        const unsigned int component) const;

  virtual Tensor<1,dim> shape_grad (const unsigned int  i,
                                    const Point<dim>   &p) const;

  virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
                                              const Point<dim> &p,
                                              const unsigned int component) const;

  virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
                                         const Point<dim> &p) const;

  virtual
  Tensor<2,dim>
  shape_grad_grad_component (const unsigned int i,
                             const Point<dim> &p,
                             const unsigned int component) const;

  virtual Tensor<3,dim> shape_3rd_derivative (const unsigned int  i,
                                              const Point<dim>   &p) const;

  virtual Tensor<3,dim> shape_3rd_derivative_component (const unsigned int i,
                                                        const Point<dim>   &p,
                                                        const unsigned int component) const;

  virtual Tensor<4,dim> shape_4th_derivative (const unsigned int  i,
                                              const Point<dim>   &p) const;

  virtual Tensor<4,dim> shape_4th_derivative_component (const unsigned int i,
                                                        const Point<dim>   &p,
                                                        const unsigned int component) const;

  virtual void
  get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                            FullMatrix<double>           &matrix) const;

  virtual const FiniteElement<dim,spacedim> &
  base_element (const unsigned int index) const;

  virtual bool has_support_on_face (const unsigned int shape_index,
                                    const unsigned int face_index) const;

  virtual const FullMatrix<double> &
  get_restriction_matrix (const unsigned int child,
                          const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;

  virtual const FullMatrix<double> &
  get_prolongation_matrix (const unsigned int child,
                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;

  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;

  virtual
  Point<dim>
  unit_support_point (const unsigned int index) const;

  virtual
  Point<dim-1>
  unit_face_support_point (const unsigned int index) const;

  virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
  get_constant_modes () const;

  virtual bool hp_constraints_are_implemented () const;

  virtual void
  get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                                 FullMatrix<double>       &matrix) const;


  virtual void
  get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                                    const unsigned int        subface,
                                    FullMatrix<double>       &matrix) const;

  virtual
  std::vector<std::pair<unsigned int, unsigned int> >
  hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;

  virtual
  std::vector<std::pair<unsigned int, unsigned int> >
  hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;

  virtual
  std::vector<std::pair<unsigned int, unsigned int> >
  hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const;

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

  virtual
  void
  convert_generalized_support_point_values_to_dof_values (const std::vector<Vector<double> > &support_point_values,
                                                          std::vector<double>                &dof_values) const;

  virtual std::size_t memory_consumption () const;

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;

  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;

  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;

  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;

  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;

  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;

  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();

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



// 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> ({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(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 //DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif
/*----------------------------  fe_system.h  ---------------------------*/