mapping_data_on_the_fly.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2014 - 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_data_on_the_fly_h
#define dealii_matrix_free_mapping_data_on_the_fly_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/base/aligned_vector.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/mapping_info.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/mapping_q1.h>


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  namespace MatrixFreeFunctions
  {
    template <int dim, typename Number=double>
    class MappingDataOnTheFly
    {
    public:
      MappingDataOnTheFly (const Mapping<dim> &mapping,
                           const Quadrature<1> &quadrature,
                           const UpdateFlags update_flags);

      MappingDataOnTheFly (const Quadrature<1> &quadrature,
                           const UpdateFlags update_flags);

      void reinit(typename ::Triangulation<dim>::cell_iterator cell);

      bool is_initialized() const;

      typename ::Triangulation<dim>::cell_iterator get_cell () const;

      const ::FEValues<dim> &get_fe_values () const;

      const MappingInfoStorage<dim,dim,Number> &
      get_data_storage() const;

      const Quadrature<1> &
      get_quadrature () const;

    private:
      typename ::Triangulation<dim>::cell_iterator present_cell;

      FE_Nothing<dim> fe_dummy;

      ::FEValues<dim> fe_values;

      const Quadrature<1> quadrature_1d;

      MappingInfoStorage<dim,dim,Number> mapping_info_storage;
    };


    /*----------------------- Inline functions ----------------------------------*/

    template <int dim, typename Number>
    inline
    MappingDataOnTheFly<dim,Number>::MappingDataOnTheFly (const Mapping<dim> &mapping,
                                                          const Quadrature<1> &quadrature,
                                                          const UpdateFlags update_flags)
      :
      fe_values(mapping, fe_dummy, Quadrature<dim>(quadrature),
                MappingInfo<dim,Number>::compute_update_flags(update_flags)),
      quadrature_1d(quadrature)
    {
      mapping_info_storage.descriptor.resize(1);
      mapping_info_storage.descriptor[0].initialize(quadrature);
      mapping_info_storage.data_index_offsets.resize(1);
      mapping_info_storage.JxW_values.resize(fe_values.n_quadrature_points);
      mapping_info_storage.jacobians[0].resize(fe_values.n_quadrature_points);
      if (update_flags & update_quadrature_points)
        {
          mapping_info_storage.quadrature_point_offsets.resize(1, 0);
          mapping_info_storage.quadrature_points.resize(fe_values.n_quadrature_points);
        }
      if (fe_values.get_update_flags() & update_normal_vectors)
        {
          mapping_info_storage.normal_vectors.resize(fe_values.n_quadrature_points);
          mapping_info_storage.normals_times_jacobians[0].resize(fe_values.n_quadrature_points);
        }
      Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
             ExcNotImplemented());
    }



    template <int dim, typename Number>
    inline
    MappingDataOnTheFly<dim,Number>::MappingDataOnTheFly (const Quadrature<1> &quadrature,
                                                          const UpdateFlags update_flags)
      :
      MappingDataOnTheFly(::StaticMappingQ1<dim,dim>::mapping,
                          quadrature, update_flags)
    {}



    template <int dim, typename Number>
    inline
    void
    MappingDataOnTheFly<dim,Number>::reinit(typename ::Triangulation<dim>::cell_iterator cell)
    {
      if (present_cell == cell)
        return;
      present_cell = cell;
      fe_values.reinit(present_cell);
      for (unsigned int q=0; q<fe_values.get_quadrature().size(); ++q)
        {
          if (fe_values.get_update_flags() & update_JxW_values)
            mapping_info_storage.JxW_values[q] = fe_values.JxW(q);
          if (fe_values.get_update_flags() & update_jacobians)
            {
              Tensor<2,dim> jac = fe_values.jacobian(q);
              jac = invert(transpose(jac));
              for (unsigned int d=0; d<dim; ++d)
                for (unsigned int e=0; e<dim; ++e)
                  mapping_info_storage.jacobians[0][q][d][e] = jac[d][e];
            }
          if (fe_values.get_update_flags() & update_quadrature_points)
            for (unsigned int d=0; d<dim; ++d)
              mapping_info_storage.quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
          if (fe_values.get_update_flags() & update_normal_vectors)
            {
              for (unsigned int d=0; d<dim; ++d)
                mapping_info_storage.normal_vectors[q][d] = fe_values.normal_vector(q)[d];
              mapping_info_storage.normals_times_jacobians[0][q] =
                mapping_info_storage.normal_vectors[q] *
                mapping_info_storage.jacobians[0][q];
            }
        }
    }



    template <int dim, typename Number>
    inline
    bool
    MappingDataOnTheFly<dim,Number>::is_initialized() const
    {
      return present_cell != typename ::Triangulation<dim>::cell_iterator();
    }



    template <int dim, typename Number>
    inline
    typename ::Triangulation<dim>::cell_iterator
    MappingDataOnTheFly<dim,Number>::get_cell() const
    {
      return fe_values.get_cell();
    }



    template <int dim, typename Number>
    inline
    const ::FEValues<dim> &
    MappingDataOnTheFly<dim,Number>::get_fe_values() const
    {
      return fe_values;
    }



    template <int dim, typename Number>
    inline
    const MappingInfoStorage<dim,dim,Number> &
    MappingDataOnTheFly<dim,Number>::get_data_storage() const
    {
      return mapping_info_storage;
    }



    template <int dim, typename Number>
    inline
    const Quadrature<1> &
    MappingDataOnTheFly<dim,Number>::get_quadrature() const
    {
      return quadrature_1d;
    }


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


DEAL_II_NAMESPACE_CLOSE

#endif