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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#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/local_results.h>
#include <deal.II/meshworker/dof_info.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:

    typedef IntegrationInfo<dim, spacedim> CellInfo;

    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> &)
  {}


}

DEAL_II_NAMESPACE_CLOSE

#endif