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