mapping_info.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2020 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/config.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,
              typename VectorizedArrayType>
    struct MappingInfoStorage
    {
      static_assert(
        std::is_same<Number, typename VectorizedArrayType::value_type>::value,
        "Type of Number and of VectorizedArrayType do not match.");

      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<1> quadrature_1d;

        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<VectorizedArrayType> JxW_values;

      AlignedVector<Tensor<1, spacedim, VectorizedArrayType>> normal_vectors;

      AlignedVector<Tensor<2, spacedim, VectorizedArrayType>> jacobians[2];

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

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

      AlignedVector<unsigned int> quadrature_point_offsets;

      AlignedVector<Point<spacedim, VectorizedArrayType>> quadrature_points;

      void
      clear_data_fields();

      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 TaskInfo &task_info) const;

      std::size_t
      memory_consumption() const;
    };



    template <int dim, typename Number, typename VectorizedArrayType>
    struct MappingInfo
    {
      void
      initialize(
        const ::Triangulation<dim> &                        tria,
        const std::vector<std::pair<unsigned int, unsigned int>> &cells,
        const FaceInfo<VectorizedArrayType::size()> &             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);

      void
      update_mapping(
        const ::Triangulation<dim> &                        tria,
        const std::vector<std::pair<unsigned int, unsigned int>> &cells,
        const FaceInfo<VectorizedArrayType::size()> &             faces,
        const std::vector<unsigned int> &active_fe_index,
        const Mapping<dim> &             mapping);

      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;

      UpdateFlags update_flags_cells;

      UpdateFlags update_flags_boundary_faces;

      UpdateFlags update_flags_inner_faces;

      UpdateFlags update_flags_faces_by_cells;

      std::vector<GeometryType> cell_type;

      std::vector<GeometryType> face_type;

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

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

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

      SmartPointer<const Mapping<dim>> mapping;

      void
      compute_mapping_q(
        const ::Triangulation<dim> &                        tria,
        const std::vector<std::pair<unsigned int, unsigned int>> &cells,
        const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
          &faces);

      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);

      void
      initialize_faces(
        const ::Triangulation<dim> &                        tria,
        const std::vector<std::pair<unsigned int, unsigned int>> &cells,
        const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
          &                 faces,
        const Mapping<dim> &mapping);

      void
      initialize_faces_by_cells(
        const ::Triangulation<dim> &                        tria,
        const std::vector<std::pair<unsigned int, unsigned int>> &cells,
        const Mapping<dim> &                                      mapping);

      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, typename>
    struct MappingInfoCellsOrFaces;

    template <int dim, typename Number, typename VectorizedArrayType>
    struct MappingInfoCellsOrFaces<dim, Number, false, VectorizedArrayType>
    {
      static const MappingInfoStorage<dim, dim, Number, VectorizedArrayType> *
      get(const MappingInfo<dim, Number, VectorizedArrayType> &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, typename VectorizedArrayType>
    struct MappingInfoCellsOrFaces<dim, Number, true, VectorizedArrayType>
    {
      static const MappingInfoStorage<dim - 1, dim, Number, VectorizedArrayType>
        *
        get(const MappingInfo<dim, Number, VectorizedArrayType> &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,
              typename VectorizedArrayType = VectorizedArray<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, VectorizedArrayType::size(), Number> &t1,
        const Tensor<1, VectorizedArrayType::size(), Number> &t2) const;

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

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

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

      Number tolerance;
    };



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

    template <int structdim,
              int spacedim,
              typename Number,
              typename VectorizedArrayType>
    inline unsigned int
    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
      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, typename VectorizedArrayType>
    inline GeometryType
    MappingInfo<dim, Number, VectorizedArrayType>::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