fe_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2019 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_tools_H
#define dealii_fe_tools_H



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

#include <deal.II/base/exceptions.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/fe/component_mask.h>

#include <deal.II/lac/la_parallel_vector.h>

#include <string>
#include <vector>


DEAL_II_NAMESPACE_OPEN

template <typename number>
class FullMatrix;
template <int dim>
class Quadrature;
template <int dim, int spacedim>
class FiniteElement;
template <int dim, int spacedim>
class DoFHandler;
template <int dim>
class FiniteElementData;
template <typename number>
class AffineConstraints;





namespace FETools
{
  template <int dim, int spacedim = dim>
  class FEFactoryBase : public Subscriptor
  {
  public:
    virtual std::unique_ptr<FiniteElement<dim, spacedim>>
    get(const unsigned int degree) const = 0;

    virtual std::unique_ptr<FiniteElement<dim, spacedim>>
    get(const Quadrature<1> &quad) const = 0;

    virtual ~FEFactoryBase() override = default;
  };

  template <class FE>
  class FEFactory : public FEFactoryBase<FE::dimension, FE::space_dimension>
  {
  public:
    virtual std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>>
    get(const unsigned int degree) const override;

    virtual std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>>
    get(const Quadrature<1> &quad) const override;
  };

  template <int dim, int spacedim>
  void
  compute_component_wise(const FiniteElement<dim, spacedim> &    fe,
                         std::vector<unsigned int> &             renumbering,
                         std::vector<std::vector<unsigned int>> &start_indices);

  template <int dim, int spacedim>
  void
  compute_block_renumbering(const FiniteElement<dim, spacedim> &  fe,
                            std::vector<types::global_dof_index> &renumbering,
                            std::vector<types::global_dof_index> &block_data,
                            bool return_start_indices = true);

  template <int dim, typename number, int spacedim>
  void
  get_interpolation_matrix(const FiniteElement<dim, spacedim> &fe1,
                           const FiniteElement<dim, spacedim> &fe2,
                           FullMatrix<number> &interpolation_matrix);

  template <int dim, typename number, int spacedim>
  void
  get_back_interpolation_matrix(const FiniteElement<dim, spacedim> &fe1,
                                const FiniteElement<dim, spacedim> &fe2,
                                FullMatrix<number> &interpolation_matrix);

  template <int dim, typename number, int spacedim>
  void
  get_interpolation_difference_matrix(const FiniteElement<dim, spacedim> &fe1,
                                      const FiniteElement<dim, spacedim> &fe2,
                                      FullMatrix<number> &difference_matrix);

  template <int dim, typename number, int spacedim>
  void
  get_projection_matrix(const FiniteElement<dim, spacedim> &fe1,
                        const FiniteElement<dim, spacedim> &fe2,
                        FullMatrix<number> &                matrix);

  template <int dim, int spacedim>
  FullMatrix<double>
  compute_node_matrix(const FiniteElement<dim, spacedim> &fe);

  template <int dim, typename number, int spacedim>
  void
  compute_embedding_matrices(
    const FiniteElement<dim, spacedim> &          fe,
    std::vector<std::vector<FullMatrix<number>>> &matrices,
    const bool                                    isotropic_only = false,
    const double                                  threshold      = 1.e-12);

  template <int dim, typename number, int spacedim>
  void
  compute_face_embedding_matrices(
    const FiniteElement<dim, spacedim> &fe,
    FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
    const unsigned int face_coarse,
    const unsigned int face_fine,
    const double       threshold = 1.e-12);

  template <int dim, typename number, int spacedim>
  void
  compute_projection_matrices(
    const FiniteElement<dim, spacedim> &          fe,
    std::vector<std::vector<FullMatrix<number>>> &matrices,
    const bool                                    isotropic_only = false);

  template <int dim, int spacedim>
  void
  compute_projection_from_quadrature_points_matrix(
    const FiniteElement<dim, spacedim> &fe,
    const Quadrature<dim> &             lhs_quadrature,
    const Quadrature<dim> &             rhs_quadrature,
    FullMatrix<double> &                X);

  template <int dim, int spacedim>
  void
  compute_interpolation_to_quadrature_points_matrix(
    const FiniteElement<dim, spacedim> &fe,
    const Quadrature<dim> &             quadrature,
    FullMatrix<double> &                I_q);

  template <int dim>
  void
  compute_projection_from_quadrature_points(
    const FullMatrix<double> &         projection_matrix,
    const std::vector<Tensor<1, dim>> &vector_of_tensors_at_qp,
    std::vector<Tensor<1, dim>> &      vector_of_tensors_at_nodes);



  template <int dim>
  void
  compute_projection_from_quadrature_points(
    const FullMatrix<double> &                  projection_matrix,
    const std::vector<SymmetricTensor<2, dim>> &vector_of_tensors_at_qp,
    std::vector<SymmetricTensor<2, dim>> &      vector_of_tensors_at_nodes);



  template <int dim, int spacedim>
  void
  compute_projection_from_face_quadrature_points_matrix(
    const FiniteElement<dim, spacedim> &fe,
    const Quadrature<dim - 1> &         lhs_quadrature,
    const Quadrature<dim - 1> &         rhs_quadrature,
    const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
    const unsigned int                                              face,
    FullMatrix<double> &                                            X);



  template <int dim, int spacedim, typename number>
  void
  convert_generalized_support_point_values_to_dof_values(
    const FiniteElement<dim, spacedim> &finite_element,
    const std::vector<Vector<number>> & support_point_values,
    std::vector<number> &               dof_values);




  template <int dim,
            int spacedim,
            template <int, int> class DoFHandlerType1,
            template <int, int> class DoFHandlerType2,
            class InVector,
            class OutVector>
  void
  interpolate(const DoFHandlerType1<dim, spacedim> &dof1,
              const InVector &                      u1,
              const DoFHandlerType2<dim, spacedim> &dof2,
              OutVector &                           u2);

  template <int dim,
            int spacedim,
            template <int, int> class DoFHandlerType1,
            template <int, int> class DoFHandlerType2,
            class InVector,
            class OutVector>
  void
  interpolate(
    const DoFHandlerType1<dim, spacedim> &                   dof1,
    const InVector &                                         u1,
    const DoFHandlerType2<dim, spacedim> &                   dof2,
    const AffineConstraints<typename OutVector::value_type> &constraints,
    OutVector &                                              u2);

  template <int dim,
            template <int, int> class DoFHandlerType,
            class InVector,
            class OutVector,
            int spacedim>
  void
  back_interpolate(const DoFHandlerType<dim, spacedim> &dof1,
                   const InVector &                     u1,
                   const FiniteElement<dim, spacedim> & fe2,
                   OutVector &                          u1_interpolated);

  template <int dim, class InVector, class OutVector, int spacedim>
  void
  back_interpolate(
    const DoFHandler<dim, spacedim> &                        dof1,
    const AffineConstraints<typename OutVector::value_type> &constraints1,
    const InVector &                                         u1,
    const DoFHandler<dim, spacedim> &                        dof2,
    const AffineConstraints<typename OutVector::value_type> &constraints2,
    OutVector &                                              u1_interpolated);

  template <int dim, class InVector, class OutVector, int spacedim>
  void
  interpolation_difference(const DoFHandler<dim, spacedim> &   dof1,
                           const InVector &                    z1,
                           const FiniteElement<dim, spacedim> &fe2,
                           OutVector &                         z1_difference);

  template <int dim, class InVector, class OutVector, int spacedim>
  void
  interpolation_difference(
    const DoFHandler<dim, spacedim> &                        dof1,
    const AffineConstraints<typename OutVector::value_type> &constraints1,
    const InVector &                                         z1,
    const DoFHandler<dim, spacedim> &                        dof2,
    const AffineConstraints<typename OutVector::value_type> &constraints2,
    OutVector &                                              z1_difference);



  template <int dim, class InVector, class OutVector, int spacedim>
  void
  project_dg(const DoFHandler<dim, spacedim> &dof1,
             const InVector &                 u1,
             const DoFHandler<dim, spacedim> &dof2,
             OutVector &                      u2);

  template <int dim, class InVector, class OutVector, int spacedim>
  void
  extrapolate(const DoFHandler<dim, spacedim> &dof1,
              const InVector &                 z1,
              const DoFHandler<dim, spacedim> &dof2,
              OutVector &                      z2);

  template <int dim, class InVector, class OutVector, int spacedim>
  void
  extrapolate(
    const DoFHandler<dim, spacedim> &                        dof1,
    const InVector &                                         z1,
    const DoFHandler<dim, spacedim> &                        dof2,
    const AffineConstraints<typename OutVector::value_type> &constraints,
    OutVector &                                              z2);


  template <int dim>
  void
  hierarchic_to_lexicographic_numbering(unsigned int               degree,
                                        std::vector<unsigned int> &h2l);

  template <int dim>
  void
  hierarchic_to_lexicographic_numbering(const FiniteElementData<dim> &fe_data,
                                        std::vector<unsigned int> &   h2l);

  template <int dim>
  std::vector<unsigned int>
  hierarchic_to_lexicographic_numbering(const FiniteElementData<dim> &fe_data);

  template <int dim>
  void
  lexicographic_to_hierarchic_numbering(const FiniteElementData<dim> &fe_data,
                                        std::vector<unsigned int> &   l2h);

  template <int dim>
  std::vector<unsigned int>
  lexicographic_to_hierarchic_numbering(const FiniteElementData<dim> &fe_data);


  namespace Compositing
  {
    template <int dim, int spacedim>
    FiniteElementData<dim>
    multiply_dof_numbers(
      const std::vector<const FiniteElement<dim, spacedim> *> &fes,
      const std::vector<unsigned int> &                        multiplicities,
      const bool do_tensor_product = true);

    template <int dim, int spacedim>
    FiniteElementData<dim>
    multiply_dof_numbers(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &fe_systems);

    template <int dim, int spacedim>
    FiniteElementData<dim>
    multiply_dof_numbers(const FiniteElement<dim, spacedim> *fe1,
                         const unsigned int                  N1,
                         const FiniteElement<dim, spacedim> *fe2 = nullptr,
                         const unsigned int                  N2  = 0,
                         const FiniteElement<dim, spacedim> *fe3 = nullptr,
                         const unsigned int                  N3  = 0,
                         const FiniteElement<dim, spacedim> *fe4 = nullptr,
                         const unsigned int                  N4  = 0,
                         const FiniteElement<dim, spacedim> *fe5 = nullptr,
                         const unsigned int                  N5  = 0);

    template <int dim, int spacedim>
    std::vector<bool>
    compute_restriction_is_additive_flags(
      const std::vector<const FiniteElement<dim, spacedim> *> &fes,
      const std::vector<unsigned int> &                        multiplicities);

    template <int dim, int spacedim>
    std::vector<bool>
    compute_restriction_is_additive_flags(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &fe_systems);

    template <int dim, int spacedim>
    std::vector<bool>
    compute_restriction_is_additive_flags(
      const FiniteElement<dim, spacedim> *fe1,
      const unsigned int                  N1,
      const FiniteElement<dim, spacedim> *fe2 = nullptr,
      const unsigned int                  N2  = 0,
      const FiniteElement<dim, spacedim> *fe3 = nullptr,
      const unsigned int                  N3  = 0,
      const FiniteElement<dim, spacedim> *fe4 = nullptr,
      const unsigned int                  N4  = 0,
      const FiniteElement<dim, spacedim> *fe5 = nullptr,
      const unsigned int                  N5  = 0);


    template <int dim, int spacedim>
    std::vector<ComponentMask>
    compute_nonzero_components(
      const std::vector<const FiniteElement<dim, spacedim> *> &fes,
      const std::vector<unsigned int> &                        multiplicities,
      const bool do_tensor_product = true);

    template <int dim, int spacedim>
    std::vector<ComponentMask>
    compute_nonzero_components(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &fe_systems);

    template <int dim, int spacedim>
    std::vector<ComponentMask>
    compute_nonzero_components(
      const FiniteElement<dim, spacedim> *fe1,
      const unsigned int                  N1,
      const FiniteElement<dim, spacedim> *fe2               = nullptr,
      const unsigned int                  N2                = 0,
      const FiniteElement<dim, spacedim> *fe3               = nullptr,
      const unsigned int                  N3                = 0,
      const FiniteElement<dim, spacedim> *fe4               = nullptr,
      const unsigned int                  N4                = 0,
      const FiniteElement<dim, spacedim> *fe5               = nullptr,
      const unsigned int                  N5                = 0,
      const bool                          do_tensor_product = true);

    template <int dim, int spacedim>
    void
    build_cell_tables(
      std::vector<std::pair<std::pair<unsigned int, unsigned int>,
                            unsigned int>> &system_to_base_table,
      std::vector<std::pair<unsigned int, unsigned int>>
        &                                   system_to_component_table,
      std::vector<std::pair<std::pair<unsigned int, unsigned int>,
                            unsigned int>> &component_to_base_table,
      const FiniteElement<dim, spacedim> &  finite_element,
      const bool                            do_tensor_product = true);

    template <int dim, int spacedim>
    void
    build_face_tables(
      std::vector<std::pair<std::pair<unsigned int, unsigned int>,
                            unsigned int>> &face_system_to_base_table,
      std::vector<std::pair<unsigned int, unsigned int>>
        &                                 face_system_to_component_table,
      const FiniteElement<dim, spacedim> &finite_element,
      const bool                          do_tensor_product = true);

  } // namespace Compositing


  template <int dim, int spacedim = dim>
  std::unique_ptr<FiniteElement<dim, spacedim>>
  get_fe_by_name(const std::string &name);


  template <int dim>
  DEAL_II_DEPRECATED FiniteElement<dim, dim> *
                     get_fe_from_name(const std::string &name);


  template <int dim, int spacedim>
  void
  add_fe_name(const std::string &                 name,
              const FEFactoryBase<dim, spacedim> *factory);

  DeclException1(ExcInvalidFEName,
                 std::string,
                 << "Can't re-generate a finite element from the string '"
                 << arg1 << "'.");

  DeclException2(ExcInvalidFEDimension,
                 char,
                 int,
                 << "The dimension " << arg1
                 << " in the finite element string must match "
                 << "the space dimension " << arg2 << ".");

  DeclException0(ExcInvalidFE);

  DeclException0(ExcFENotPrimitive);
  DeclException0(ExcTriangulationMismatch);

  DeclExceptionMsg(ExcHangingNodesNotAllowed,
                   "You are using continuous elements on a grid with "
                   "hanging nodes but without providing hanging node "
                   "constraints. Use the respective function with "
                   "additional AffineConstraints argument(s), instead.");
  DeclException0(ExcGridNotRefinedAtLeastOnce);
  DeclException4(ExcMatrixDimensionMismatch,
                 int,
                 int,
                 int,
                 int,
                 << "This is a " << arg1 << "x" << arg2 << " matrix, "
                 << "but should be a " << arg3 << "x" << arg4 << " matrix.");

  DeclException1(ExcLeastSquaresError,
                 double,
                 << "Least squares fit leaves a gap of " << arg1);

  DeclException2(ExcNotGreaterThan,
                 int,
                 int,
                 << arg1 << " must be greater than " << arg2);
} // namespace FETools


#ifndef DOXYGEN

namespace FETools
{
  template <class FE>
  std::unique_ptr<FiniteElement<FE::dimension, FE::space_dimension>>
  FEFactory<FE>::get(const unsigned int degree) const
  {
    return std_cxx14::make_unique<FE>(degree);
  }

  namespace Compositing
  {
    template <int dim, int spacedim>
    std::vector<bool>
    compute_restriction_is_additive_flags(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &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);

      return compute_restriction_is_additive_flags(fes, multiplicities);
    }



    template <int dim, int spacedim>
    FiniteElementData<dim>
    multiply_dof_numbers(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &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);

      return multiply_dof_numbers(fes, multiplicities, true);
    }



    template <int dim, int spacedim>
    std::vector<ComponentMask>
    compute_nonzero_components(
      const std::initializer_list<
        std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>, unsigned int>>
        &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);

      return compute_nonzero_components(fes, multiplicities, true);
    }
  } // namespace Compositing
} // namespace FETools

#endif

DEAL_II_NAMESPACE_CLOSE

#endif /* dealii_fe_tools_H */