fe.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2017 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_h
#define dealii__fe_h

#include <deal.II/base/config.h>
#include <deal.II/fe/fe_base.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/block_mask.h>
#include <deal.II/fe/mapping.h>

DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim> class FEValuesBase;
template <int dim, int spacedim> class FEValues;
template <int dim, int spacedim> class FEFaceValues;
template <int dim, int spacedim> class FESubfaceValues;
template <int dim, int spacedim> class FESystem;

template <int dim, int spacedim=dim>
class FiniteElement : public Subscriptor,
  public FiniteElementData<dim>
{
public:
  static const unsigned int space_dimension = spacedim;

  class InternalDataBase
  {
  private:
    InternalDataBase (const InternalDataBase &);

  public:
    InternalDataBase ();

    virtual ~InternalDataBase ();

    UpdateFlags          update_each;

    virtual std::size_t memory_consumption () const;
  };

public:
  FiniteElement (const FiniteElementData<dim>     &fe_data,
                 const std::vector<bool>          &restriction_is_additive_flags,
                 const std::vector<ComponentMask> &nonzero_components);

  virtual ~FiniteElement ();

  virtual FiniteElement<dim,spacedim> *clone() const = 0;

  virtual std::string get_name () const = 0;

  const FiniteElement<dim,spacedim> &operator[] (const unsigned int fe_index) 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 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;

  bool prolongation_is_implemented () const;

  bool isotropic_prolongation_is_implemented () const;

  bool restriction_is_implemented () const;

  bool isotropic_restriction_is_implemented () const;


  bool restriction_is_additive (const unsigned int index) const;

  const FullMatrix<double> &constraints (const ::internal::SubfaceCase<dim> &subface_case=::internal::SubfaceCase<dim>::case_isotropic) const;

  bool constraints_are_implemented (const ::internal::SubfaceCase<dim> &subface_case=::internal::SubfaceCase<dim>::case_isotropic) const;


  virtual bool hp_constraints_are_implemented () const;


  virtual void
  get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
                            FullMatrix<double>       &matrix) 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;


  bool operator == (const FiniteElement<dim,spacedim> &) const;

  std::pair<unsigned int, unsigned int>
  system_to_component_index (const unsigned int index) const;

  unsigned int component_to_system_index(const unsigned int component,
                                         const unsigned int index) const;

  std::pair<unsigned int, unsigned int>
  face_system_to_component_index (const unsigned int index) const;

  unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index,
                                                           const bool face_orientation,
                                                           const bool face_flip,
                                                           const bool face_rotation) 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;

  unsigned int adjust_line_dof_index_for_line_orientation (const unsigned int index,
                                                           const bool line_orientation) const;

  const ComponentMask &
  get_nonzero_components (const unsigned int i) const;

  unsigned int
  n_nonzero_components (const unsigned int i) const;

  bool is_primitive () const;

  bool
  is_primitive (const unsigned int i) const;

  unsigned int n_base_elements () const;

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

  unsigned int
  element_multiplicity (const unsigned int index) const;

  std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
  system_to_base_index (const unsigned int index) const;

  std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
  face_system_to_base_index (const unsigned int index) const;

  types::global_dof_index first_block_of_base (const unsigned int b) const;

  std::pair<unsigned int, unsigned int>
  component_to_base_index (const unsigned int component) const;


  std::pair<unsigned int,unsigned int>
  block_to_base_index (const unsigned int block) const;

  std::pair<unsigned int,types::global_dof_index>
  system_to_block_index (const unsigned int component) const;

  unsigned int
  component_to_block_index (const unsigned int component) const;


  ComponentMask
  component_mask (const FEValuesExtractors::Scalar &scalar) const;

  ComponentMask
  component_mask (const FEValuesExtractors::Vector &vector) const;

  ComponentMask
  component_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;

  ComponentMask
  component_mask (const BlockMask &block_mask) const;

  BlockMask
  block_mask (const FEValuesExtractors::Scalar &scalar) const;

  BlockMask
  block_mask (const FEValuesExtractors::Vector &vector) const;

  BlockMask
  block_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;

  BlockMask
  block_mask (const ComponentMask &component_mask) const;

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


  const std::vector<Point<dim> > &
  get_unit_support_points () const;

  bool has_support_points () const;

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

  const std::vector<Point<dim-1> > &
  get_unit_face_support_points () const;

  bool has_face_support_points () const;

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

  const std::vector<Point<dim> > &
  get_generalized_support_points () const;

  bool has_generalized_support_points () const;

  const std::vector<Point<dim-1> > &
  get_generalized_face_support_points () const;

  bool
  has_generalized_face_support_points () const;

  GeometryPrimitive
  get_associated_geometry_primitive (const unsigned int cell_dof_index) const;


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

  virtual
  void
  interpolate(std::vector<double>       &local_dofs,
              const std::vector<double> &values) const DEAL_II_DEPRECATED;

  virtual
  void
  interpolate(std::vector<double>                &local_dofs,
              const std::vector<Vector<double> > &values,
              const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;

  virtual
  void
  interpolate(std::vector<double>                                         &local_dofs,
              const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;


  virtual std::size_t memory_consumption () const;

  DeclException1 (ExcShapeFunctionNotPrimitive,
                  int,
                  << "The shape function with index " << arg1
                  << " is not primitive, i.e. it is vector-valued and "
                  << "has more than one non-zero vector component. This "
                  << "function cannot be called for these shape functions. "
                  << "Maybe you want to use the same function with the "
                  << "_component suffix?");
  DeclException0 (ExcFENotPrimitive);
  DeclExceptionMsg (ExcUnitShapeValuesDoNotExist,
                    "You are trying to access the values or derivatives of shape functions "
                    "on the reference cell of an element that does not define its shape "
                    "functions through mapping from the reference cell. Consequently, "
                    "you cannot ask for shape function values or derivatives on the "
                    "reference cell.");

  DeclExceptionMsg (ExcFEHasNoSupportPoints,
                    "You are trying to access the support points of a finite "
                    "element that either has no support points at all, or for "
                    "which the corresponding tables have not been implemented.");

  DeclExceptionMsg (ExcEmbeddingVoid,
                    "You are trying to access the matrices that describe how "
                    "to embed a finite element function on one cell into the "
                    "finite element space on one of its children (i.e., the "
                    "'embedding' or 'prolongation' matrices). However, the "
                    "current finite element can either not define this sort of "
                    "operation, or it has not yet been implemented.");

  DeclExceptionMsg (ExcProjectionVoid,
                    "You are trying to access the matrices that describe how "
                    "to restrict a finite element function from the children "
                    "of one cell to the finite element space defined on their "
                    "parent (i.e., the 'restriction' or 'projection' matrices). "
                    "However, the current finite element can either not define "
                    "this sort of operation, or it has not yet been "
                    "implemented.");

  DeclException2 (ExcWrongInterfaceMatrixSize,
                  int, int,
                  << "The interface matrix has a size of " << arg1
                  << "x" << arg2
                  << ", which is not reasonable for the current element "
                  "in the present dimension.");
  DeclException0 (ExcInterpolationNotImplemented);

protected:

  void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false,
                                                    const bool isotropic_prolongation_only=false);

  std::vector<std::vector<FullMatrix<double> > > restriction;

  std::vector<std::vector<FullMatrix<double> > > prolongation;

  FullMatrix<double> interface_constraints;

  std::vector<Point<dim> > unit_support_points;

  std::vector<Point<dim-1> > unit_face_support_points;

  std::vector<Point<dim> > generalized_support_points;

  std::vector<Point<dim-1> > generalized_face_support_points;

  Table<2,int> adjust_quad_dof_index_for_face_orientation_table;

  std::vector<int> adjust_line_dof_index_for_line_orientation_table;

  std::vector<std::pair<unsigned int, unsigned int> > system_to_component_table;

  std::vector<std::pair<unsigned int, unsigned int> > face_system_to_component_table;

  std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
  system_to_base_table;

  std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
  face_system_to_base_table;

  BlockIndices base_to_block_indices;

  std::vector<std::pair<std::pair<unsigned int, unsigned int>, unsigned int> >
  component_to_base_table;

  const std::vector<bool> restriction_is_additive_flags;

  const std::vector<ComponentMask> nonzero_components;

  const std::vector<unsigned int> n_nonzero_components_table;

  const bool cached_primitivity;

  TableIndices<2>
  interface_constraints_size () const;

  static
  std::vector<unsigned int>
  compute_n_nonzero_components (const std::vector<ComponentMask> &nonzero_components);

  virtual
  UpdateFlags
  requires_update_flags (const UpdateFlags update_flags) const = 0;

  virtual
  InternalDataBase *
  get_data (const UpdateFlags                                                    update_flags,
            const Mapping<dim,spacedim>                                         &mapping,
            const Quadrature<dim>                                               &quadrature,
            ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const = 0;

  virtual
  InternalDataBase *
  get_face_data (const UpdateFlags                                                    update_flags,
                 const Mapping<dim,spacedim>                                         &mapping,
                 const Quadrature<dim-1>                                             &quadrature,
                 ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const;

  virtual
  InternalDataBase *
  get_subface_data (const UpdateFlags                                                    update_flags,
                    const Mapping<dim,spacedim>                                         &mapping,
                    const Quadrature<dim-1>                                             &quadrature,
                    ::internal::FEValues::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::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
                  const InternalDataBase                                              &fe_internal,
                  ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const = 0;

  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::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
                       const InternalDataBase                                              &fe_internal,
                       ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const = 0;

  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::FEValues::MappingRelatedData<dim, spacedim> &mapping_data,
                          const InternalDataBase                                              &fe_internal,
                          ::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const = 0;

  friend class InternalDataBase;
  friend class FEValuesBase<dim,spacedim>;
  friend class FEValues<dim,spacedim>;
  friend class FEFaceValues<dim,spacedim>;
  friend class FESubfaceValues<dim,spacedim>;
  friend class FESystem<dim,spacedim>;

  // explicitly check for sensible template arguments, but not on windows
  // because MSVC creates bogus warnings during normal compilation
#ifdef DEAL_II_WITH_CXX11
#ifndef DEAL_II_MSVC
  static_assert (dim<=spacedim,
                 "The dimension <dim> of a FiniteElement must be less than or "
                 "equal to the space dimension <spacedim> in which it lives.");
#endif
#endif

};


//----------------------------------------------------------------------//


template <int dim, int spacedim>
inline
const FiniteElement<dim,spacedim> &
FiniteElement<dim,spacedim>::operator[] (const unsigned int fe_index) const
{
  (void)fe_index;
  Assert (fe_index == 0,
          ExcMessage ("A fe_index of zero is the only index allowed here"));
  return *this;
}



template <int dim, int spacedim>
inline
std::pair<unsigned int,unsigned int>
FiniteElement<dim,spacedim>::system_to_component_index (const unsigned int index) const
{
  Assert (index < system_to_component_table.size(),
          ExcIndexRange(index, 0, system_to_component_table.size()));
  Assert (is_primitive (index),
          ( typename FiniteElement<dim,spacedim>::ExcShapeFunctionNotPrimitive(index)) );
  return system_to_component_table[index];
}



template <int dim, int spacedim>
inline
unsigned int
FiniteElement<dim,spacedim>::n_base_elements () const
{
  return base_to_block_indices.size();
}



template <int dim, int spacedim>
inline
unsigned int
FiniteElement<dim,spacedim>::element_multiplicity (const unsigned int index) const
{
  return static_cast<unsigned int>(base_to_block_indices.block_size(index));
}



template <int dim, int spacedim>
inline
unsigned int
FiniteElement<dim,spacedim>::component_to_system_index (const unsigned int component,
                                                        const unsigned int index) const
{
  AssertIndexRange(component, this->n_components());
  const std::vector<std::pair<unsigned int, unsigned int> >::const_iterator
  it = std::find(system_to_component_table.begin(), system_to_component_table.end(),
                 std::pair<unsigned int, unsigned int>(component, index));

  Assert(it != system_to_component_table.end(),
         ExcMessage ("You are asking for the number of the shape function "
                     "within a system element that corresponds to vector "
                     "component " + Utilities::int_to_string(component) + " and within this to "
                     "index " + Utilities::int_to_string(index) + ". But no such "
                     "shape function exists."));
  return std::distance(system_to_component_table.begin(), it);
}



template <int dim, int spacedim>
inline
std::pair<unsigned int,unsigned int>
FiniteElement<dim,spacedim>::face_system_to_component_index (const unsigned int index) const
{
  Assert(index < face_system_to_component_table.size(),
         ExcIndexRange(index, 0, face_system_to_component_table.size()));

  // in debug mode, check whether the
  // function is primitive, since
  // otherwise the result may have no
  // meaning
  //
  // since the primitivity tables are
  // all geared towards cell dof
  // indices, rather than face dof
  // indices, we have to work a
  // little bit...
  //
  // in 1d, the face index is equal
  // to the cell index
  Assert (is_primitive(this->face_to_cell_index(index, 0)),
          (typename FiniteElement<dim,spacedim>::ExcShapeFunctionNotPrimitive(index)) );

  return face_system_to_component_table[index];
}



template <int dim, int spacedim>
inline
std::pair<std::pair<unsigned int,unsigned int>,unsigned int>
FiniteElement<dim,spacedim>::system_to_base_index (const unsigned int index) const
{
  Assert (index < system_to_base_table.size(),
          ExcIndexRange(index, 0, system_to_base_table.size()));
  return system_to_base_table[index];
}




template <int dim, int spacedim>
inline
std::pair<std::pair<unsigned int,unsigned int>,unsigned int>
FiniteElement<dim,spacedim>::face_system_to_base_index (const unsigned int index) const
{
  Assert(index < face_system_to_base_table.size(),
         ExcIndexRange(index, 0, face_system_to_base_table.size()));
  return face_system_to_base_table[index];
}



template <int dim, int spacedim>
inline
types::global_dof_index
FiniteElement<dim,spacedim>::first_block_of_base (const unsigned int index) const
{
  return base_to_block_indices.block_start(index);
}



template <int dim, int spacedim>
inline
std::pair<unsigned int,unsigned int>
FiniteElement<dim,spacedim>::component_to_base_index (const unsigned int index) const
{
  Assert(index < component_to_base_table.size(),
         ExcIndexRange(index, 0, component_to_base_table.size()));

  return component_to_base_table[index].first;
}



template <int dim, int spacedim>
inline
std::pair<unsigned int,unsigned int>
FiniteElement<dim,spacedim>::block_to_base_index (const unsigned int index) const
{
  return base_to_block_indices.global_to_local(index);
}



template <int dim, int spacedim>
inline
std::pair<unsigned int,types::global_dof_index>
FiniteElement<dim,spacedim>::system_to_block_index (const unsigned int index) const
{
  Assert (index < this->dofs_per_cell,
          ExcIndexRange(index, 0, this->dofs_per_cell));
  // The block is computed simply as
  // first block of this base plus
  // the index within the base blocks
  return std::pair<unsigned int, types::global_dof_index>(
           first_block_of_base(system_to_base_table[index].first.first)
           + system_to_base_table[index].first.second,
           system_to_base_table[index].second);
}



template <int dim, int spacedim>
inline
bool
FiniteElement<dim,spacedim>::restriction_is_additive (const unsigned int index) const
{
  Assert(index < this->dofs_per_cell,
         ExcIndexRange(index, 0, this->dofs_per_cell));
  return restriction_is_additive_flags[index];
}



template <int dim, int spacedim>
inline
const ComponentMask &
FiniteElement<dim,spacedim>::get_nonzero_components (const unsigned int i) const
{
  Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));
  return nonzero_components[i];
}



template <int dim, int spacedim>
inline
unsigned int
FiniteElement<dim,spacedim>::n_nonzero_components (const unsigned int i) const
{
  Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));
  return n_nonzero_components_table[i];
}



template <int dim, int spacedim>
inline
bool
FiniteElement<dim,spacedim>::is_primitive () const
{
  return cached_primitivity;
}



template <int dim, int spacedim>
inline
bool
FiniteElement<dim,spacedim>::is_primitive (const unsigned int i) const
{
  Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));

  // return primitivity of a shape
  // function by checking whether it
  // has more than one non-zero
  // component or not. we could cache
  // this value in an array of bools,
  // but accessing a bit-vector (as
  // std::vector<bool> is) is
  // probably more expensive than
  // just comparing against 1
  //
  // for good measure, short circuit the test
  // if the entire FE is primitive
  return (is_primitive() ||
          (n_nonzero_components_table[i] == 1));
}



template <int dim, int spacedim>
inline
GeometryPrimitive
FiniteElement<dim,spacedim>::get_associated_geometry_primitive (const unsigned int cell_dof_index) const
{
  Assert (cell_dof_index < this->dofs_per_cell,
          ExcIndexRange (cell_dof_index, 0, this->dofs_per_cell));

  // just go through the usual cases, taking into account how DoFs
  // are enumerated on the reference cell
  if (cell_dof_index < this->first_line_index)
    return GeometryPrimitive::vertex;
  else if (cell_dof_index < this->first_quad_index)
    return GeometryPrimitive::line;
  else if (cell_dof_index < this->first_hex_index)
    return GeometryPrimitive::quad;
  else
    return GeometryPrimitive::hex;
}



DEAL_II_NAMESPACE_CLOSE

#endif