integration_info.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#ifndef dealii_mesh_worker_integration_info_h
#define dealii_mesh_worker_integration_info_h

#include <deal.II/base/config.h>

#include <deal.II/base/quadrature_lib.h>

#include <deal.II/dofs/block_info.h>

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

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/local_results.h>
#include <deal.II/meshworker/vector_selector.h>

#include <memory>

DEAL_II_NAMESPACE_OPEN

namespace MeshWorker
{
  template <int dim, int spacedim = dim>
  class IntegrationInfo
  {
  private:
    std::vector<std::shared_ptr<FEValuesBase<dim, spacedim>>> fevalv;

  public:
    static const unsigned int dimension       = dim;
    static const unsigned int space_dimension = spacedim;

    IntegrationInfo();

    IntegrationInfo(const IntegrationInfo<dim, spacedim> &other);

    template <class FEVALUES>
    void
    initialize(const FiniteElement<dim, spacedim> &            el,
               const Mapping<dim, spacedim> &                  mapping,
               const Quadrature<FEVALUES::integral_dimension> &quadrature,
               const UpdateFlags                               flags,
               const BlockInfo *local_block_info = nullptr);

    void
    initialize_data(const std::shared_ptr<VectorDataBase<dim, spacedim>> &data);

    void
    clear();

    const FiniteElement<dim, spacedim> &
    finite_element() const;

    bool multigrid;

    const FEValuesBase<dim, spacedim> &
    fe_values() const;


    const FEValuesBase<dim, spacedim> &
    fe_values(const unsigned int i) const;

    std::vector<std::vector<std::vector<double>>> values;

    std::vector<std::vector<std::vector<Tensor<1, spacedim>>>> gradients;

    std::vector<std::vector<std::vector<Tensor<2, spacedim>>>> hessians;

    template <typename number>
    void
    reinit(const DoFInfo<dim, spacedim, number> &i);

    template <typename number>
    void
    fill_local_data(const DoFInfo<dim, spacedim, number> &info,
                    bool                                  split_fevalues);

    std::shared_ptr<VectorDataBase<dim, spacedim>> global_data;

    std::size_t
    memory_consumption() const;

  private:
    SmartPointer<const FiniteElement<dim, spacedim>,
                 IntegrationInfo<dim, spacedim>>
      fe_pointer;

    template <typename TYPE>
    void
    fill_local_data(std::vector<std::vector<std::vector<TYPE>>> &data,
                    VectorSelector &                             selector,
                    bool split_fevalues) const;
    unsigned int n_components;
  };

  template <int dim, int spacedim = dim>
  class IntegrationInfoBox
  {
  public:
    using CellInfo = IntegrationInfo<dim, spacedim>;

    IntegrationInfoBox();

    void
    initialize(const FiniteElement<dim, spacedim> &el,
               const Mapping<dim, spacedim> &      mapping,
               const BlockInfo *                   block_info = nullptr);

    template <typename VectorType>
    void
    initialize(const FiniteElement<dim, spacedim> &el,
               const Mapping<dim, spacedim> &      mapping,
               const AnyData &                     data,
               const VectorType &                  dummy,
               const BlockInfo *                   block_info = nullptr);
    template <typename VectorType>
    void
    initialize(const FiniteElement<dim, spacedim> &el,
               const Mapping<dim, spacedim> &      mapping,
               const AnyData &                     data,
               const MGLevelObject<VectorType> &   dummy,
               const BlockInfo *                   block_info = nullptr);
    /* @{ */

    void
    initialize_update_flags(bool neighbor_geometry = false);

    void
    add_update_flags_all(const UpdateFlags flags);

    void
    add_update_flags_cell(const UpdateFlags flags);

    void
    add_update_flags_boundary(const UpdateFlags flags);

    void
    add_update_flags_face(const UpdateFlags flags);

    void
    add_update_flags(const UpdateFlags flags,
                     const bool        cell     = true,
                     const bool        boundary = true,
                     const bool        face     = true,
                     const bool        neighbor = true);

    void
    initialize_gauss_quadrature(unsigned int n_cell_points,
                                unsigned int n_boundary_points,
                                unsigned int n_face_points,
                                const bool   force = true);

    std::size_t
    memory_consumption() const;

    UpdateFlags cell_flags;
    UpdateFlags boundary_flags;

    UpdateFlags face_flags;

    UpdateFlags neighbor_flags;

    Quadrature<dim> cell_quadrature;

    Quadrature<dim - 1> boundary_quadrature;

    Quadrature<dim - 1> face_quadrature;
    /* @} */

    /* @{ */

    MeshWorker::VectorSelector cell_selector;

    MeshWorker::VectorSelector boundary_selector;

    MeshWorker::VectorSelector face_selector;

    std::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim>> cell_data;
    std::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim>> boundary_data;
    std::shared_ptr<MeshWorker::VectorDataBase<dim, spacedim>> face_data;
    /* @} */

    /* @{ */
    template <class DOFINFO>
    void
    post_cell(const DoFInfoBox<dim, DOFINFO> &);

    template <class DOFINFO>
    void
    post_faces(const DoFInfoBox<dim, DOFINFO> &);

    CellInfo cell;
    CellInfo boundary;
    CellInfo face;
    CellInfo subface;
    CellInfo neighbor;

    /* @} */
  };


  //----------------------------------------------------------------------//

  template <int dim, int sdim>
  inline IntegrationInfo<dim, sdim>::IntegrationInfo()
    : fevalv(0)
    , multigrid(false)
    , global_data(std::make_shared<VectorDataBase<dim, sdim>>())
    , n_components(numbers::invalid_unsigned_int)
  {}


  template <int dim, int sdim>
  inline IntegrationInfo<dim, sdim>::IntegrationInfo(
    const IntegrationInfo<dim, sdim> &other)
    : multigrid(other.multigrid)
    , values(other.values)
    , gradients(other.gradients)
    , hessians(other.hessians)
    , global_data(other.global_data)
    , fe_pointer(other.fe_pointer)
    , n_components(other.n_components)
  {
    fevalv.resize(other.fevalv.size());
    for (unsigned int i = 0; i < other.fevalv.size(); ++i)
      {
        const FEValuesBase<dim, sdim> &p = *other.fevalv[i];
        const FEValues<dim, sdim> *    pc =
          dynamic_cast<const FEValues<dim, sdim> *>(&p);
        const FEFaceValues<dim, sdim> *pf =
          dynamic_cast<const FEFaceValues<dim, sdim> *>(&p);
        const FESubfaceValues<dim, sdim> *ps =
          dynamic_cast<const FESubfaceValues<dim, sdim> *>(&p);

        if (pc != nullptr)
          fevalv[i] =
            std::make_shared<FEValues<dim, sdim>>(pc->get_mapping(),
                                                  pc->get_fe(),
                                                  pc->get_quadrature(),
                                                  pc->get_update_flags());
        else if (pf != nullptr)
          fevalv[i] =
            std::make_shared<FEFaceValues<dim, sdim>>(pf->get_mapping(),
                                                      pf->get_fe(),
                                                      pf->get_quadrature(),
                                                      pf->get_update_flags());
        else if (ps != nullptr)
          fevalv[i] = std::make_shared<FESubfaceValues<dim, sdim>>(
            ps->get_mapping(),
            ps->get_fe(),
            ps->get_quadrature(),
            ps->get_update_flags());
        else
          Assert(false, ExcInternalError());
      }
  }



  template <int dim, int sdim>
  template <class FEVALUES>
  inline void
  IntegrationInfo<dim, sdim>::initialize(
    const FiniteElement<dim, sdim> &                el,
    const Mapping<dim, sdim> &                      mapping,
    const Quadrature<FEVALUES::integral_dimension> &quadrature,
    const UpdateFlags                               flags,
    const BlockInfo *                               block_info)
  {
    fe_pointer = &el;
    if (block_info == nullptr || block_info->local().size() == 0)
      {
        fevalv.resize(1);
        fevalv[0] = std::make_shared<FEVALUES>(mapping, el, quadrature, flags);
      }
    else
      {
        fevalv.resize(el.n_base_elements());
        for (unsigned int i = 0; i < fevalv.size(); ++i)
          fevalv[i] = std::make_shared<FEVALUES>(mapping,
                                                 el.base_element(i),
                                                 quadrature,
                                                 flags);
      }
    n_components = el.n_components();
  }


  template <int dim, int spacedim>
  inline const FiniteElement<dim, spacedim> &
  IntegrationInfo<dim, spacedim>::finite_element() const
  {
    Assert(fe_pointer != nullptr, ExcNotInitialized());
    return *fe_pointer;
  }

  template <int dim, int spacedim>
  inline const FEValuesBase<dim, spacedim> &
  IntegrationInfo<dim, spacedim>::fe_values() const
  {
    AssertDimension(fevalv.size(), 1);
    return *fevalv[0];
  }


  template <int dim, int spacedim>
  inline const FEValuesBase<dim, spacedim> &
  IntegrationInfo<dim, spacedim>::fe_values(unsigned int i) const
  {
    Assert(i < fevalv.size(), ExcIndexRange(i, 0, fevalv.size()));
    return *fevalv[i];
  }


  template <int dim, int spacedim>
  template <typename number>
  inline void
  IntegrationInfo<dim, spacedim>::reinit(
    const DoFInfo<dim, spacedim, number> &info)
  {
    for (unsigned int i = 0; i < fevalv.size(); ++i)
      {
        FEValuesBase<dim, spacedim> &febase = *fevalv[i];
        if (info.sub_number != numbers::invalid_unsigned_int)
          {
            // This is a subface
            FESubfaceValues<dim, spacedim> &fe =
              dynamic_cast<FESubfaceValues<dim, spacedim> &>(febase);
            fe.reinit(info.cell, info.face_number, info.sub_number);
          }
        else if (info.face_number != numbers::invalid_unsigned_int)
          {
            // This is a face
            FEFaceValues<dim, spacedim> &fe =
              dynamic_cast<FEFaceValues<dim, spacedim> &>(febase);
            fe.reinit(info.cell, info.face_number);
          }
        else
          {
            // This is a cell
            FEValues<dim, spacedim> &fe =
              dynamic_cast<FEValues<dim, spacedim> &>(febase);
            fe.reinit(info.cell);
          }
      }

    const bool split_fevalues = info.block_info != nullptr;
    if (!global_data->empty())
      fill_local_data(info, split_fevalues);
  }



  //----------------------------------------------------------------------//

  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::initialize_gauss_quadrature(unsigned int cp,
                                                             unsigned int bp,
                                                             unsigned int fp,
                                                             bool         force)
  {
    if (force || cell_quadrature.size() == 0)
      cell_quadrature = QGauss<dim>(cp);
    if (force || boundary_quadrature.size() == 0)
      boundary_quadrature = QGauss<dim - 1>(bp);
    if (force || face_quadrature.size() == 0)
      face_quadrature = QGauss<dim - 1>(fp);
  }


  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_all(const UpdateFlags flags)
  {
    add_update_flags(flags, true, true, true, true);
  }


  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_cell(const UpdateFlags flags)
  {
    add_update_flags(flags, true, false, false, false);
  }


  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_boundary(
    const UpdateFlags flags)
  {
    add_update_flags(flags, false, true, false, false);
  }


  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::add_update_flags_face(const UpdateFlags flags)
  {
    add_update_flags(flags, false, false, true, true);
  }


  template <int dim, int sdim>
  inline void
  IntegrationInfoBox<dim, sdim>::initialize(const FiniteElement<dim, sdim> &el,
                                            const Mapping<dim, sdim> &mapping,
                                            const BlockInfo *block_info)
  {
    initialize_update_flags();
    initialize_gauss_quadrature((cell_flags & update_values) ?
                                  (el.tensor_degree() + 1) :
                                  el.tensor_degree(),
                                (boundary_flags & update_values) ?
                                  (el.tensor_degree() + 1) :
                                  el.tensor_degree(),
                                (face_flags & update_values) ?
                                  (el.tensor_degree() + 1) :
                                  el.tensor_degree(),
                                false);

    cell.template initialize<FEValues<dim, sdim>>(
      el, mapping, cell_quadrature, cell_flags, block_info);
    boundary.template initialize<FEFaceValues<dim, sdim>>(
      el, mapping, boundary_quadrature, boundary_flags, block_info);
    face.template initialize<FEFaceValues<dim, sdim>>(
      el, mapping, face_quadrature, face_flags, block_info);
    subface.template initialize<FESubfaceValues<dim, sdim>>(
      el, mapping, face_quadrature, face_flags, block_info);
    neighbor.template initialize<FEFaceValues<dim, sdim>>(
      el, mapping, face_quadrature, neighbor_flags, block_info);
  }


  template <int dim, int sdim>
  template <typename VectorType>
  void
  IntegrationInfoBox<dim, sdim>::initialize(const FiniteElement<dim, sdim> &el,
                                            const Mapping<dim, sdim> &mapping,
                                            const AnyData &           data,
                                            const VectorType &,
                                            const BlockInfo *block_info)
  {
    initialize(el, mapping, block_info);
    std::shared_ptr<VectorData<VectorType, dim, sdim>> p;
    VectorDataBase<dim, sdim> *                        pp;

    p = std::make_shared<VectorData<VectorType, dim, sdim>>(cell_selector);
    // Public member function of parent class was not found without
    // explicit cast
    pp = &*p;
    pp->initialize(data);
    cell_data = p;
    cell.initialize_data(p);

    p  = std::make_shared<VectorData<VectorType, dim, sdim>>(boundary_selector);
    pp = &*p;
    pp->initialize(data);
    boundary_data = p;
    boundary.initialize_data(p);

    p  = std::make_shared<VectorData<VectorType, dim, sdim>>(face_selector);
    pp = &*p;
    pp->initialize(data);
    face_data = p;
    face.initialize_data(p);
    subface.initialize_data(p);
    neighbor.initialize_data(p);
  }

  template <int dim, int sdim>
  template <typename VectorType>
  void
  IntegrationInfoBox<dim, sdim>::initialize(const FiniteElement<dim, sdim> &el,
                                            const Mapping<dim, sdim> &mapping,
                                            const AnyData &           data,
                                            const MGLevelObject<VectorType> &,
                                            const BlockInfo *block_info)
  {
    initialize(el, mapping, block_info);
    std::shared_ptr<MGVectorData<VectorType, dim, sdim>> p;
    VectorDataBase<dim, sdim> *                          pp;

    p = std::make_shared<MGVectorData<VectorType, dim, sdim>>(cell_selector);
    // Public member function of parent class was not found without
    // explicit cast
    pp = &*p;
    pp->initialize(data);
    cell_data = p;
    cell.initialize_data(p);

    p =
      std::make_shared<MGVectorData<VectorType, dim, sdim>>(boundary_selector);
    pp = &*p;
    pp->initialize(data);
    boundary_data = p;
    boundary.initialize_data(p);

    p  = std::make_shared<MGVectorData<VectorType, dim, sdim>>(face_selector);
    pp = &*p;
    pp->initialize(data);
    face_data = p;
    face.initialize_data(p);
    subface.initialize_data(p);
    neighbor.initialize_data(p);
  }

  template <int dim, int sdim>
  template <class DOFINFO>
  void
  IntegrationInfoBox<dim, sdim>::post_cell(const DoFInfoBox<dim, DOFINFO> &)
  {}


  template <int dim, int sdim>
  template <class DOFINFO>
  void
  IntegrationInfoBox<dim, sdim>::post_faces(const DoFInfoBox<dim, DOFINFO> &)
  {}


} // namespace MeshWorker

DEAL_II_NAMESPACE_CLOSE

#endif