mapping_info.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 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_matrix_free_mapping_info_h
#define dealii_matrix_free_mapping_info_h


#include <deal.II/base/exceptions.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/base/aligned_vector.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/matrix_free/helper_functions.h>
#include <deal.II/matrix_free/face_info.h>

#include <memory>


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  namespace MatrixFreeFunctions
  {
    enum GeometryType : unsigned char
    {
      cartesian  = 0,
      affine     = 1,
      flat_faces = 2,
      general    = 3
    };



    template <int structdim, int spacedim, typename Number>
    struct MappingInfoStorage
    {
      struct QuadratureDescriptor
      {
        QuadratureDescriptor();

        void initialize(const Quadrature<1> &quadrature_1d,
                        const UpdateFlags    update_flags_inner_faces = update_default);

        std::size_t memory_consumption() const;

        unsigned int n_q_points;

        Quadrature<structdim> quadrature;

        std::array<AlignedVector<Number>, structdim> tensor_quadrature_weights;

        AlignedVector<Number> quadrature_weights;

        ::Table<2,unsigned int> face_orientations;
      };

      std::vector<QuadratureDescriptor> descriptor;

      AlignedVector<unsigned int> data_index_offsets;

      AlignedVector<VectorizedArray<Number> > JxW_values;

      AlignedVector<Tensor<1,spacedim,VectorizedArray<Number> > > normal_vectors;

      AlignedVector<Tensor<2,spacedim,VectorizedArray<Number> > > jacobians[2];

      AlignedVector<Tensor<1,spacedim *(spacedim+1)/2,
                    Tensor<1,spacedim,VectorizedArray<Number> > > > jacobian_gradients[2];

      AlignedVector<Tensor<1,spacedim,VectorizedArray<Number> > > normals_times_jacobians [2];

      AlignedVector<unsigned int> quadrature_point_offsets;

      AlignedVector<Point<spacedim,VectorizedArray<Number> > > quadrature_points;

      unsigned int
      quad_index_from_n_q_points (const unsigned int n_q_points) const;

      template <typename StreamType>
      void print_memory_consumption(StreamType     &out,
                                    const SizeInfo &task_info) const;

      std::size_t memory_consumption () const;
    };



    template <int dim, typename Number>
    struct MappingInfo
    {
      MappingInfo();

      void initialize (const ::Triangulation<dim>                &tria,
                       const std::vector<std::pair<unsigned int,unsigned int> >  &cells,
                       const FaceInfo<VectorizedArray<Number>::n_array_elements> &faces,
                       const std::vector<unsigned int>         &active_fe_index,
                       const Mapping<dim>                      &mapping,
                       const std::vector<::hp::QCollection<1> >  &quad,
                       const UpdateFlags                        update_flags_cells,
                       const UpdateFlags                        update_flags_boundary_faces,
                       const UpdateFlags                        update_flags_inner_faces,
                       const UpdateFlags                        update_flags_faces_by_cells);

      GeometryType get_cell_type (const unsigned int cell_chunk_no) const;

      void clear ();

      std::size_t memory_consumption() const;

      template <typename StreamType>
      void print_memory_consumption(StreamType     &out,
                                    const TaskInfo &task_info) const;

      std::vector<GeometryType> cell_type;

      std::vector<GeometryType> face_type;

      std::vector<MappingInfoStorage<dim,dim,Number> > cell_data;

      std::vector<MappingInfoStorage<dim-1,dim,Number> > face_data;

      std::vector<MappingInfoStorage<dim-1,dim,Number> > face_data_by_cells;

      void initialize_cells (const ::Triangulation<dim>                &tria,
                             const std::vector<std::pair<unsigned int,unsigned int> > &cells,
                             const std::vector<unsigned int>         &active_fe_index,
                             const Mapping<dim>                      &mapping,
                             const std::vector<::hp::QCollection<1> >  &quad,
                             const UpdateFlags                        update_flags_cells);

      void initialize_faces (const ::Triangulation<dim>        &tria,
                             const std::vector<std::pair<unsigned int,unsigned int> > &cells,
                             const std::vector<FaceToCellTopology<VectorizedArray<Number>::n_array_elements> > &faces,
                             const Mapping<dim>                      &mapping,
                             const std::vector<::hp::QCollection<1> >  &quad,
                             const UpdateFlags                        update_flags_boundary_faces,
                             const UpdateFlags                        update_flags_inner_faces);

      void initialize_faces_by_cells
      (const ::Triangulation<dim>        &tria,
       const std::vector<std::pair<unsigned int,unsigned int> > &cells,
       const Mapping<dim>                      &mapping,
       const std::vector<::hp::QCollection<1> > &quad,
       const UpdateFlags                        update_flags_faces_by_cells);

      static UpdateFlags
      compute_update_flags (const UpdateFlags                        update_flags,
                            const std::vector<::hp::QCollection<1> > &quad =
                              std::vector<::hp::QCollection<1> >());
    };



    template <int, typename, bool> struct MappingInfoCellsOrFaces;

    template <int dim, typename Number>
    struct MappingInfoCellsOrFaces<dim,Number,false>
    {
      static const MappingInfoStorage<dim,dim,Number> *
      get(const MappingInfo<dim,Number> &mapping_info,
          const unsigned int quad_no)
      {
        AssertIndexRange(quad_no, mapping_info.cell_data.size());
        return &mapping_info.cell_data[quad_no];
      }
    };

    template <int dim, typename Number>
    struct MappingInfoCellsOrFaces<dim,Number,true>
    {
      static const MappingInfoStorage<dim-1,dim,Number> *
      get(const MappingInfo<dim,Number> &mapping_info,
          const unsigned int quad_no)
      {
        AssertIndexRange(quad_no, mapping_info.face_data.size());
        return &mapping_info.face_data[quad_no];
      }
    };



    template <typename Number>
    struct FPArrayComparator
    {
      FPArrayComparator (const Number scaling);

      bool operator() (const std::vector<Number> &v1,
                       const std::vector<Number> &v2) const;

      bool operator ()(const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t1,
                       const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t2) const;

      template <int dim>
      bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
                       const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;

      template <int dim>
      bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
                       const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;

      Number tolerance;
    };



    /* ------------------- inline functions ----------------------------- */

    template <int structdim, int spacedim, typename Number>
    inline
    unsigned int
    MappingInfoStorage<structdim,spacedim,Number>
    ::quad_index_from_n_q_points (const unsigned int n_q_points) const
    {
      for (unsigned int i=0; i<descriptor.size(); ++i)
        if (n_q_points == descriptor[i].n_q_points)
          return i;
      return 0;
    }



    template <int dim, typename Number>
    inline
    GeometryType
    MappingInfo<dim,Number>::get_cell_type (const unsigned int cell_no) const
    {
      AssertIndexRange (cell_no, cell_type.size());
      return cell_type[cell_no];
    }

  } // end of namespace MatrixFreeFunctions
} // end of namespace internal

DEAL_II_NAMESPACE_CLOSE

#endif