fe_tools.h

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



#include <deal.II/base/config.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/lac/parallel_vector.h>



#include <vector>
#include <string>


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;
class ConstraintMatrix;





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

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

  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 ConstraintMatrix              &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 ConstraintMatrix         &constraints1,
                         const InVector                 &u1,
                         const DoFHandler<dim,spacedim> &dof2,
                         const ConstraintMatrix         &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 ConstraintMatrix         &constraints1,
                                const InVector                 &z1,
                                const DoFHandler<dim,spacedim> &dof2,
                                const ConstraintMatrix         &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 ConstraintMatrix         &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);

  }


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

  DeclException1 (ExcHangingNodesNotAllowed,
                  int,
                  << "You are using continuous elements on a grid with "
                  << "hanging nodes but without providing hanging node "
                  << "constraints. Use the respective function with "
                  << "additional ConstraintMatrix argument(s), instead."
                  << (arg1?"":""));
  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);
}


#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);
    }
  } //Compositing
} //FETools

#endif

DEAL_II_NAMESPACE_CLOSE

#endif /* dealii_fe_tools_H */