mapping_q_generic.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_mapping_q_generic_h
#define dealii_mapping_q_generic_h


#include <deal.II/base/config.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/table.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/matrix_free/shape_info.h>

#include <array>
#include <cmath>

DEAL_II_NAMESPACE_OPEN

template <int,int> class MappingQ;




template <int dim, int spacedim=dim>
class MappingQGeneric : public Mapping<dim,spacedim>
{
public:
  MappingQGeneric (const unsigned int polynomial_degree);

  MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping);

  // for documentation, see the Mapping base class
  virtual
  std::unique_ptr<Mapping<dim,spacedim>> clone () const override;

  unsigned int get_degree () const;

  virtual
  bool preserves_vertex_locations () const override;

  // for documentation, see the Mapping base class
  virtual
  Point<spacedim>
  transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                               const Point<dim>                                 &p) const override;

  // for documentation, see the Mapping base class
  virtual
  Point<dim>
  transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                               const Point<spacedim>                            &p) const override;

  // for documentation, see the Mapping base class
  virtual
  void
  transform (const ArrayView<const Tensor<1,dim> >                  &input,
             const MappingType                                       type,
             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
             const ArrayView<Tensor<1,spacedim> >                   &output) const override;

  // for documentation, see the Mapping base class
  virtual
  void
  transform (const ArrayView<const DerivativeForm<1, dim, spacedim> > &input,
             const MappingType                                         type,
             const typename Mapping<dim,spacedim>::InternalDataBase   &internal,
             const ArrayView<Tensor<2,spacedim> >                     &output) const override;

  // for documentation, see the Mapping base class
  virtual
  void
  transform (const ArrayView<const Tensor<2, dim> >                 &input,
             const MappingType                                       type,
             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
             const ArrayView<Tensor<2,spacedim> >                   &output) const override;

  // for documentation, see the Mapping base class
  virtual
  void
  transform (const ArrayView<const DerivativeForm<2, dim, spacedim> > &input,
             const MappingType                                         type,
             const typename Mapping<dim,spacedim>::InternalDataBase   &internal,
             const ArrayView<Tensor<3,spacedim> >                     &output) const override;

  // for documentation, see the Mapping base class
  virtual
  void
  transform (const ArrayView<const Tensor<3, dim> >                 &input,
             const MappingType                                       type,
             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
             const ArrayView<Tensor<3,spacedim> >                   &output) const override;

public:
  class InternalData : public Mapping<dim,spacedim>::InternalDataBase
  {
  public:
    InternalData(const unsigned int polynomial_degree);

    void
    initialize (const UpdateFlags      update_flags,
                const Quadrature<dim> &quadrature,
                const unsigned int     n_original_q_points);

    void
    initialize_face (const UpdateFlags      update_flags,
                     const Quadrature<dim> &quadrature,
                     const unsigned int     n_original_q_points);

    void compute_shape_function_values (const std::vector<Point<dim> > &unit_points);


    const double &shape (const unsigned int qpoint,
                         const unsigned int shape_nr) const;

    double &shape (const unsigned int qpoint,
                   const unsigned int shape_nr);

    const Tensor<1,dim> &derivative (const unsigned int qpoint,
                                     const unsigned int shape_nr) const;

    Tensor<1,dim> &derivative (const unsigned int qpoint,
                               const unsigned int shape_nr);

    const Tensor<2,dim> &second_derivative (const unsigned int qpoint,
                                            const unsigned int shape_nr) const;

    Tensor<2,dim> &second_derivative (const unsigned int qpoint,
                                      const unsigned int shape_nr);

    const Tensor<3,dim> &third_derivative (const unsigned int qpoint,
                                           const unsigned int shape_nr) const;

    Tensor<3,dim> &third_derivative (const unsigned int qpoint,
                                     const unsigned int shape_nr);

    const Tensor<4,dim> &fourth_derivative (const unsigned int qpoint,
                                            const unsigned int shape_nr) const;

    Tensor<4,dim> &fourth_derivative (const unsigned int qpoint,
                                      const unsigned int shape_nr);

    virtual
    std::size_t memory_consumption () const override;

    std::vector<double> shape_values;

    std::vector<Tensor<1,dim> > shape_derivatives;

    std::vector<Tensor<2,dim> > shape_second_derivatives;

    std::vector<Tensor<3,dim> > shape_third_derivatives;

    std::vector<Tensor<4,dim> > shape_fourth_derivatives;

    std::array<std::vector<Tensor<1,dim> >, GeometryInfo<dim>::faces_per_cell *(dim-1)> unit_tangentials;

    const unsigned int polynomial_degree;

    const unsigned int n_shape_functions;

    /*
     * The default line support points. Is used in when the shape function
     * values are computed.
     *
     * The number of quadrature points depends on the degree of this
     * class, and it matches the number of degrees of freedom of an
     * FE_Q<1>(this->degree).
     */
    QGaussLobatto<1> line_support_points;

    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<double>> shape_info;

    mutable AlignedVector<VectorizedArray<double> > scratch;

    mutable AlignedVector<VectorizedArray<double> > values_dofs;

    mutable AlignedVector<VectorizedArray<double> > values_quad;

    mutable AlignedVector<VectorizedArray<double> > gradients_quad;

    mutable AlignedVector<VectorizedArray<double> > hessians_quad;

    bool tensor_product_quadrature;

    mutable std::vector<DerivativeForm<1,dim, spacedim > >  covariant;

    mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant;

    mutable std::vector<std::vector<Tensor<1,spacedim> > > aux;

    mutable std::vector<Point<spacedim> > mapping_support_points;

    mutable typename Triangulation<dim,spacedim>::cell_iterator cell_of_current_support_points;

    mutable std::vector<double> volume_elements;
  };


  // documentation can be found in Mapping::requires_update_flags()
  virtual
  UpdateFlags
  requires_update_flags (const UpdateFlags update_flags) const override;

  // documentation can be found in Mapping::get_data()
  virtual
  std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
  get_data (const UpdateFlags,
            const Quadrature<dim> &quadrature) const override;

  // documentation can be found in Mapping::get_face_data()
  virtual
  std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
  get_face_data (const UpdateFlags flags,
                 const Quadrature<dim-1>& quadrature) const override;

  // documentation can be found in Mapping::get_subface_data()
  virtual
  std::unique_ptr<typename Mapping<dim,spacedim>::InternalDataBase>
  get_subface_data (const UpdateFlags flags,
                    const Quadrature<dim-1>& quadrature) const override;

  // documentation can be found in Mapping::fill_fe_values()
  virtual
  CellSimilarity::Similarity
  fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator     &cell,
                  const CellSimilarity::Similarity                               cell_similarity,
                  const Quadrature<dim>                                         &quadrature,
                  const typename Mapping<dim,spacedim>::InternalDataBase        &internal_data,
                  ::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const override;

  // documentation can be found in Mapping::fill_fe_face_values()
  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 typename Mapping<dim,spacedim>::InternalDataBase        &internal_data,
                       ::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const override;

  // documentation can be found in Mapping::fill_fe_subface_values()
  virtual
  void
  fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator     &cell,
                          const unsigned int                                             face_no,
                          const unsigned int                                             subface_no,
                          const Quadrature<dim-1>                                       &quadrature,
                          const typename Mapping<dim,spacedim>::InternalDataBase        &internal_data,
                          ::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const override;

protected:

  const unsigned int polynomial_degree;

  /*
   * The default line support points. These are used when computing
   * the location in real space of the support points on lines and
   * quads, which are asked to the Manifold<dim,spacedim> class.
   *
   * The number of quadrature points depends on the degree of this
   * class, and it matches the number of degrees of freedom of an
   * FE_Q<1>(this->degree).
   */
  QGaussLobatto<1> line_support_points;

  const std::unique_ptr<FE_Q<dim> > fe_q;

  std::vector<Table<2,double> > support_point_weights_perimeter_to_interior;

  Table<2,double> support_point_weights_cell;

  virtual
  std::vector<Point<spacedim> >
  compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;

  Point<dim>
  transform_real_to_unit_cell_internal (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                                        const Point<spacedim> &p,
                                        const Point<dim> &initial_p_unit) const;

  virtual
  void
  add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                           std::vector<Point<spacedim> > &a) const;

  virtual
  void
  add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                          std::vector<Point<spacedim> > &a) const;

  template <int, int> friend class MappingQ;
};



/*----------------------------------------------------------------------*/

#ifndef DOXYGEN

template <int dim, int spacedim>
inline
const double &
MappingQGeneric<dim,spacedim>::InternalData::shape (const unsigned int qpoint,
                                                    const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_values.size()));
  return shape_values [qpoint*n_shape_functions + shape_nr];
}



template <int dim, int spacedim>
inline
double &
MappingQGeneric<dim,spacedim>::InternalData::shape (const unsigned int qpoint,
                                                    const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_values.size()));
  return shape_values [qpoint*n_shape_functions + shape_nr];
}


template <int dim, int spacedim>
inline
const Tensor<1,dim> &
MappingQGeneric<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
                                                         const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_derivatives.size()));
  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
}



template <int dim, int spacedim>
inline
Tensor<1,dim> &
MappingQGeneric<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
                                                         const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_derivatives.size()));
  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
}


template <int dim, int spacedim>
inline
const Tensor<2,dim> &
MappingQGeneric<dim,spacedim>::InternalData::second_derivative (const unsigned int qpoint,
    const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_second_derivatives.size()));
  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
}


template <int dim, int spacedim>
inline
Tensor<2,dim> &
MappingQGeneric<dim,spacedim>::InternalData::second_derivative (const unsigned int qpoint,
    const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_second_derivatives.size()));
  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
}

template <int dim, int spacedim>
inline
const Tensor<3,dim> &
MappingQGeneric<dim,spacedim>::InternalData::third_derivative (const unsigned int qpoint,
    const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_third_derivatives.size()));
  return shape_third_derivatives [qpoint*n_shape_functions + shape_nr];
}


template <int dim, int spacedim>
inline
Tensor<3,dim> &
MappingQGeneric<dim,spacedim>::InternalData::third_derivative (const unsigned int qpoint,
    const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_third_derivatives.size()));
  return shape_third_derivatives [qpoint*n_shape_functions + shape_nr];
}


template <int dim, int spacedim>
inline
const Tensor<4,dim> &
MappingQGeneric<dim,spacedim>::InternalData::fourth_derivative (const unsigned int qpoint,
    const unsigned int shape_nr) const
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_fourth_derivatives.size()));
  return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr];
}


template <int dim, int spacedim>
inline
Tensor<4,dim> &
MappingQGeneric<dim,spacedim>::InternalData::fourth_derivative (const unsigned int qpoint,
    const unsigned int shape_nr)
{
  Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(),
         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
                       shape_fourth_derivatives.size()));
  return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr];
}



template <int dim, int spacedim>
inline
bool
MappingQGeneric<dim,spacedim>::preserves_vertex_locations () const
{
  return true;
}

#endif // DOXYGEN

/* -------------- declaration of explicit specializations ------------- */


DEAL_II_NAMESPACE_CLOSE

#endif