mapping_info.h

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


#include <deal.II/base/aligned_vector.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/vectorization.h>

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

#include <deal.II/hp/q_collection.h>

#include <deal.II/matrix_free/face_info.h>
#include <deal.II/matrix_free/helper_functions.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
    {
      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