mapping_info.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2016 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 <memory>


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  namespace MatrixFreeFunctions
  {
    template <int dim, typename Number>
    struct MappingInfo
    {
      static const std::size_t  n_cell_type_bits = 2;

      static const unsigned int n_cell_types = 1U<<n_cell_type_bits;

      static const unsigned int n_vector_elements = VectorizedArray<Number>::n_array_elements;

      MappingInfo();

      void initialize (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);

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

      CellType get_cell_type (const unsigned int cell_chunk_no) const;

      unsigned int get_cell_data_index (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 SizeInfo &size_info) const;

      std::vector<unsigned int> cell_type;

      AlignedVector<std::pair<Tensor<1,dim,VectorizedArray<Number> >,
                    VectorizedArray<Number> > > cartesian_data;

      AlignedVector<std::pair<Tensor<2,dim,VectorizedArray<Number> >,
                    VectorizedArray<Number> > > affine_data;

      struct MappingInfoDependent
      {
        std::vector<unsigned int> rowstart_jacobians;

        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians;

        AlignedVector<VectorizedArray<Number> > JxW_values;

        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians_grad_diag;

        AlignedVector<Tensor<1,(dim>1?dim*(dim-1)/2:1),
                      Tensor<1,dim,VectorizedArray<Number> > > > jacobians_grad_upper;

        std::vector<unsigned int> rowstart_q_points;

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

        ::hp::QCollection<dim>    quadrature;

        ::hp::QCollection<dim-1>  face_quadrature;

        std::vector<unsigned int> n_q_points;

        std::vector<unsigned int> n_q_points_face;

        std::vector<AlignedVector<VectorizedArray<Number> > > quadrature_weights;

        std::vector<unsigned int> quad_index_conversion;

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

        std::size_t memory_consumption () const;
      };

      std::vector<MappingInfoDependent> mapping_data_gen;

      bool JxW_values_initialized;

      bool second_derivatives_initialized;

      bool quadrature_points_initialized;

      struct CellData
      {
        CellData (const double jac_size);
        void resize (const unsigned int size);

        AlignedVector<Tensor<1,dim,VectorizedArray<Number> > >  quadrature_points;
        AlignedVector<Tensor<2,dim,VectorizedArray<Number> > >  general_jac;
        AlignedVector<Tensor<3,dim,VectorizedArray<Number> > >  general_jac_grad;
        Tensor<2,dim,VectorizedArray<Number> > const_jac;
        const double                           jac_size;
      };

      void evaluate_on_cell (const ::Triangulation<dim> &tria,
                             const std::pair<unsigned int,unsigned int> *cells,
                             const unsigned int  cell,
                             const unsigned int  my_q,
                             CellType (&cell_t_prev)[n_vector_elements],
                             CellType (&cell_t)[n_vector_elements],
                             ::FEValues<dim,dim> &fe_values,
                             CellData          &cell_data) const;
    };



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

    template <int dim, typename Number>
    inline
    unsigned int
    MappingInfo<dim,Number>::MappingInfoDependent::
    quad_index_from_n_q_points (const unsigned int n_q_points) const
    {
      for (unsigned int i=0; i<quad_index_conversion.size(); ++i)
        if (n_q_points == quad_index_conversion[i])
          return i;
      return 0;
    }



    template <int dim, typename Number>
    inline
    CellType
    MappingInfo<dim,Number>::get_cell_type (const unsigned int cell_no) const
    {
      AssertIndexRange (cell_no, cell_type.size());
      CellType enum_cell_type = (CellType)(cell_type[cell_no] % n_cell_types);
      Assert(enum_cell_type != undefined, ExcInternalError());
      return enum_cell_type;
    }



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

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

DEAL_II_NAMESPACE_CLOSE

#endif