mesh_loop_info_objects.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
//
// FCST: Fuel Cell Simulation Toolbox
//
// Copyright (C) 2006-2009 by Guido Kanschat
// Copyright (C) 2009-2014 by Energy Systems Design Laboratory, University of Alberta
//
// This software is distributed under the MIT License
// For more information, see the README file in /doc/LICENSE
//
// - Class: mesh_loop_info_objects.h
// - Description: Objects used for looping over mesh
// - Developers: Guido Kanschat,     Texas A&M University
//               Valentin N. Zingan, University of Alberta
//
// ----------------------------------------------------------------------------

#ifndef _FUEL_CELL_APPLICATION_CORE_MESH_LOOP_INFO_OBJECTS_H_
#define _FUEL_CELL_APPLICATION_CORE_MESH_LOOP_INFO_OBJECTS_H_

//-- deal.II
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/vector_slice.h>
#include <deal.II/lac/block_indices.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_tools.h>

//-- OpenFCST
#include <application_core/fe_vectors.h>

//-- boost libraries:
#include <boost/shared_ptr.hpp>

//----
namespace dealii
{
       template<int,int> class DoFHandler;
}

using namespace dealii;

namespace FuelCell
{
namespace ApplicationCore
{

// +++ IMPLEMENTATION +++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename TYPE>
inline
void
fill_data(const FEValuesBase<dim>&          fe_values,
          const FEVector&                   fe_vector,
          const std::vector<unsigned int>&  local_dof_indices,
          unsigned int                      first_index,
          unsigned int                      n_indices,
          std::vector< std::vector<TYPE> >& result)
{
       Assert(false, ExcNotImplemented());
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim>
inline
void
fill_data(const FEValuesBase<dim>&            fe_values,
          const FEVector&                     fe_vector,
          const std::vector<unsigned int>&    local_dof_indices,
          unsigned int                        first_index,
          unsigned int                        n_indices,
          std::vector< std::vector<double> >& result)
{
    
    VectorSlice<std::vector< std::vector<double> > > aux(result);
    fe_values.get_function_values(fe_vector,
                                  make_slice(local_dof_indices, first_index, n_indices),
                                  aux,
                                  true);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim>
inline
void
fill_data(const FEValuesBase<dim>&                     fe_values,
          const FEVector&                              fe_vector,
          const std::vector<unsigned int>&             local_dof_indices,
          unsigned int                                 first_index,
          unsigned int                                 n_indices,
          std::vector< std::vector< Tensor<1,dim> > >& result)
{
    VectorSlice< std::vector< std::vector< Tensor<1,dim> > > > aux(result);
    fe_values.get_function_gradients(fe_vector,
                                     make_slice(local_dof_indices, first_index, n_indices),
                                     aux,
                                     true);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim>
inline
void
fill_data(const FEValuesBase<dim>&                     fe_values,
          const FEVector&                              fe_vector,
          const std::vector<unsigned int>&             local_dof_indices,
          unsigned int                                 first_index,
          unsigned int                                 n_indices,
          std::vector< std::vector< Tensor<2,dim> > >& result)
{
       fe_values.get_function_hessians(fe_vector,
                                       make_slice(local_dof_indices, first_index, n_indices),
                                       result,
                                       true);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

struct BlockInfo : public Subscriptor
{
    BlockIndices global;
    
    BlockIndices local;
    
    std::vector<BlockIndices> levels;
    
    std::vector<unsigned int> base_element;
    
    std::vector<unsigned int> local_renumbering;
    
    template<int dim, int spacedim>
    void initialize(const DoFHandler<dim, spacedim>& dof_handler)
    {
        const FiniteElement<dim, spacedim>& fe = dof_handler.get_fe();
        std::vector<unsigned int> sizes(fe.n_blocks());
        DoFTools::count_dofs_per_block(dof_handler, sizes);
        global.reinit(sizes);
    }
    
    template<int dim, int spacedim>
    void initialize_local(const DoFHandler<dim, spacedim>& dof_handler)
    {
        const FiniteElement<dim, spacedim>& fe = dof_handler.get_fe();
        std::vector<unsigned int> sizes(fe.n_blocks());
        
        base_element.resize(fe.n_blocks());
        
        for(unsigned int i = 0; i < base_element.size(); ++i)
            base_element[i] = fe.block_to_base_index(i).first;
        
        local_renumbering.resize(fe.n_dofs_per_cell());
        FETools::compute_block_renumbering(fe,
                                           local_renumbering,
                                           sizes,
                                           false);
        local.reinit(sizes);
    }   
};


template<int dim, int spacedim = dim>
class DoFInfo
{
public:



  DoFInfo(const BlockInfo& block_info);

  DoFInfo(const FEVectors&,
          const BlockInfo& block_info);




  template<typename DHCellIterator>
  void reinit(const DHCellIterator& c);

  template<typename DHCellIterator, typename DHFaceIterator>
  void reinit(const DHCellIterator& c,
              const DHFaceIterator& f,
              const unsigned int    fn);

  template<typename DHCellIterator, typename DHFaceIterator>
  void reinit(const DHCellIterator& c,
              const DHFaceIterator& f,
              const unsigned int    fn,
              const unsigned int    sn);




  typename DoFHandler<dim>::cell_iterator dof_cell;

  typename DoFHandler<dim>::active_cell_iterator dof_active_cell;

  typename DoFHandler<dim>::face_iterator dof_face;

  typename Triangulation<dim>::cell_iterator cell;

  typename Triangulation<dim>::face_iterator face;

  unsigned int face_number;

  unsigned int sub_number;




  std::vector<unsigned int> indices;

  SmartPointer<const BlockInfo> block_info;


private:

  void get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator c);

  std::vector<unsigned int> indices_org;
};

// +++ IMPLEMENTATION +++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, int spacedim>
inline
DoFInfo<dim, spacedim>::DoFInfo(const BlockInfo& block_info)
:
block_info(&block_info,
           typeid(*this).name())
{ }

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, int spacedim>
inline
DoFInfo<dim, spacedim>::DoFInfo(const FEVectors&,
                                const BlockInfo& block_info)
:
block_info(&block_info,
           typeid(*this).name())
{ }

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, int spacedim>
template<typename DHCellIterator>
inline
void
DoFInfo<dim, spacedim>::reinit(const DHCellIterator& c)
{
       get_indices(c);

       dof_cell        = c;
       dof_active_cell = c;
       cell            = static_cast<typename Triangulation<dim>::cell_iterator> (c);

       face_number = static_cast<unsigned int>(-1);
       sub_number  = static_cast<unsigned int>(-1);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, int spacedim>
template<typename DHCellIterator, typename DHFaceIterator>
inline
void
DoFInfo<dim, spacedim>::reinit(const DHCellIterator& c,
                               const DHFaceIterator& f,
                               const unsigned int    fn)
{
       if(   cell.state() != IteratorState::valid || cell != static_cast<typename Triangulation<dim>::cell_iterator> (c)   )
       {
         get_indices(c);
       }

       dof_cell        = c;
       dof_active_cell = c;
       cell            = static_cast<typename Triangulation<dim>::cell_iterator> (c);

       dof_face = f;
       face     = static_cast<typename Triangulation<dim>::face_iterator> (f);

       face_number = fn;
       sub_number  = static_cast<unsigned int>(-1);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, int spacedim>
template<typename DHCellIterator, typename DHFaceIterator>
inline
void
DoFInfo<dim, spacedim>::reinit(const DHCellIterator& c,
                               const DHFaceIterator& f,
                               const unsigned int    fn,
                               const unsigned int    sn)
{
       if(   cell.state() != IteratorState::valid || cell != static_cast<typename Triangulation<dim>::cell_iterator> (c)   )
       {
         get_indices(c);
       }

       dof_cell        = c;
       dof_active_cell = c;
       cell            = static_cast<typename Triangulation<dim>::cell_iterator> (c);

       dof_face = f;
       face     = static_cast<typename Triangulation<dim>::face_iterator> (f);

       face_number = fn;
       sub_number  = sn;
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, int spacedim>
inline
void
DoFInfo<dim, spacedim>::get_indices(const typename DoFHandler<dim, spacedim>::cell_iterator c)
{
       indices.resize(c->get_fe().dofs_per_cell);

       if(block_info->local_renumbering.size() == 0)
       {
         c->get_dof_indices(indices);
       }
       else
       {
         indices_org.resize(c->get_fe().dofs_per_cell);
         c->get_dof_indices(indices_org);

         for(unsigned int i = 0; i < indices.size(); ++i)
           indices[this->block_info->local_renumbering[i]] = indices_org[i];
       }
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
class IntegrationInfo : public DoFInfo<dim>
{
public:



  IntegrationInfo(const BlockInfo& block_info);

  IntegrationInfo(const FEVectors& data,
                  const BlockInfo& block_info);

  template<typename FEVALUES>
  void initialize(const FEVALUES*                                 FE_VALUES,
                  const FiniteElement<dim>&                       fe,
                  const Mapping<dim>&                             mapping,
                  const Quadrature<FEVALUES::integral_dimension>& quadrature,
                  const UpdateFlags                               flags);

  void initialize_data(const FEVectors& data);




  template<typename DHCellIterator>
  void reinit(const DHCellIterator& c);

  template<typename DHCellIterator, typename DHFaceIterator>
  void reinit(const DHCellIterator& c,
              const DHFaceIterator& f,
              const unsigned int    fn);

  template<typename DHCellIterator, typename DHFaceIterator>
  void reinit(const DHCellIterator& c,
              const DHFaceIterator& f,
              const unsigned int    fn,
              const unsigned int    sn);

  template<typename TYPE>
  void fill_local_data(std::vector< std::vector< std::vector<TYPE> > >& data,
                       bool split_fevalues) const;




  const FEVALUESBASE& fe() const;

  const FEVALUESBASE& fe(unsigned int i) const;

  const FEVALUESBASE& get_fe_val_unsplit() const;

  void clear();




  bool multigrid;

  SmartPointer<const FEVectors> global_data;

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

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

  std::vector< std::vector< std::vector< Tensor<1,dim> > > >& derivatives;

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


private:



  std::vector< boost::shared_ptr<FEVALUESBASE> > fevalv;

  boost::shared_ptr<FEVALUESBASE> fe_val_unsplit;

};

// +++ IMPLEMENTATION +++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
inline
IntegrationInfo<dim, FEVALUESBASE>::IntegrationInfo(const BlockInfo& block_info)
:
DoFInfo<dim>(block_info),
multigrid(false),
global_data(0, typeid(*this).name()),
derivatives(gradients),
fevalv(0)
{ }

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
inline
IntegrationInfo<dim, FEVALUESBASE>::IntegrationInfo(const FEVectors& data,
                                                    const BlockInfo& block_info)
:
DoFInfo<dim>(block_info),
multigrid(false),
global_data(&data),
derivatives(gradients),
fevalv(0)
{ }

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
template<typename FEVALUES>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::initialize(const FEVALUES*,
                                               const FiniteElement<dim>&                       fe,
                                               const Mapping<dim>&                             mapping,
                                               const Quadrature<FEVALUES::integral_dimension>& quadrature,
                                               const UpdateFlags                               flags)
{
       if( this->block_info->local_renumbering.size() != 0 ) // MODE1
       {
              fevalv.resize(fe.n_base_elements());
              for(unsigned int i = 0; i < fevalv.size(); ++i)
              {
                fevalv[i] = boost::shared_ptr<FEVALUESBASE>( new FEVALUES(mapping,
                                                                          fe.base_element(i),
                                                                          quadrature,
                                                                          flags)
                                                           );
              }

              // --- fe_val_unsplit initialization ---

              fe_val_unsplit = boost::shared_ptr<FEVALUESBASE>( new FEVALUES(mapping,
                                                                             fe,
                                                                             quadrature,
                                                                             flags)
                                                              );
       }
       else // MODE2
       {
              fevalv.resize(1);
              fevalv[0] = boost::shared_ptr<FEVALUESBASE>( new FEVALUES(mapping,
                                                                        fe,
                                                                        quadrature,
                                                                        flags)
                                                         );
       }

       values.resize(global_data->n_vectors(),
                     std::vector< std::vector<double> >(fe.n_components(),
                                                        std::vector<double>(quadrature.size())
                                                       )
                    );
       gradients.resize(global_data->n_vectors(),
                        std::vector< std::vector< Tensor<1,dim> > >(fe.n_components(),
                                                                    std::vector< Tensor<1,dim> >(quadrature.size())
                                                                   )
                       );
       hessians.resize(global_data->n_vectors(),
                       std::vector< std::vector< Tensor<2,dim> > >(fe.n_components(),
                                                                   std::vector< Tensor<2,dim> >(quadrature.size())
                                                                  )
                      );
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::initialize_data(const FEVectors& data)
{
       global_data = &data;
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
template<typename DHCellIterator>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::reinit(const DHCellIterator& c)
{
       DoFInfo<dim>::reinit(c);

       for(unsigned int i = 0; i < fevalv.size(); ++i)
       {
         FEVALUESBASE&  fe_values_base = *fevalv[i];
         FEValues<dim>& fe_values      = dynamic_cast<FEValues<dim>&>(fe_values_base);
         fe_values.reinit(this->cell);
       }

       FEVALUESBASE&  fe_values_base_unsplit = *fe_val_unsplit;
       FEValues<dim>& fe_values_unsplit      = dynamic_cast<FEValues<dim>&>(fe_values_base_unsplit);
       fe_values_unsplit.reinit(this->dof_active_cell);

       const bool split_fevalues = this->block_info->local_renumbering.size() != 0;
       fill_local_data(values,
                       split_fevalues);
       fill_local_data(gradients,
                       split_fevalues);
     //fill_local_data(hessians,
     //                split_fevalues);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
template<typename DHCellIterator, typename DHFaceIterator>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::reinit(const DHCellIterator& c,
                                           const DHFaceIterator& f,
                                           const unsigned int    fn)
{
       DoFInfo<dim>::reinit(c,
                            f,
                            fn);

       for(unsigned int i = 0; i < fevalv.size(); ++i)
       {
         FEVALUESBASE&      fe_values_base = *fevalv[i];
         FEFaceValues<dim>& fe_face_values = dynamic_cast<FEFaceValues<dim>&>(fe_values_base);
         fe_face_values.reinit(this->cell,
                               fn);
       }

       FEVALUESBASE&      fe_values_base_unsplit = *fe_val_unsplit;
       FEFaceValues<dim>& fe_face_values_unsplit = dynamic_cast<FEFaceValues<dim>&>(fe_values_base_unsplit);
       fe_face_values_unsplit.reinit(this->dof_active_cell,
                                     fn);

       const bool split_fevalues = this->block_info->local_renumbering.size() != 0;
       fill_local_data(values,
                       split_fevalues);
       fill_local_data(gradients,
                       split_fevalues);
     //fill_local_data(hessians,
     //                split_fevalues);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
template<typename DHCellIterator, typename DHFaceIterator>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::reinit(const DHCellIterator& c,
                                           const DHFaceIterator& f,
                                           const unsigned int    fn,
                                           const unsigned int    sn)
{
       DoFInfo<dim>::reinit(c, f, fn, sn);

       for(unsigned int i = 0; i < fevalv.size(); ++i)
       {
         FEVALUESBASE&         fe_values_base    = *fevalv[i];
         FESubfaceValues<dim>& fe_subface_values = dynamic_cast<FESubfaceValues<dim>&>(fe_values_base);
         fe_subface_values.reinit(this->cell,
                                  fn,
                                  sn);
       }

       FEVALUESBASE&         fe_values_base_unsplit    = *fe_val_unsplit;
       FESubfaceValues<dim>& fe_subface_values_unsplit = dynamic_cast<FESubfaceValues<dim>&>(fe_values_base_unsplit);
       fe_subface_values_unsplit.reinit(this->dof_active_cell,
                                        fn,
                                        sn);

       const bool split_fevalues = this->block_info->local_renumbering.size() != 0;
       fill_local_data(values,
                       split_fevalues);
       fill_local_data(gradients,
                       split_fevalues);
     //fill_local_data(hessians,
     //                split_fevalues);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
template<typename TYPE>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::fill_local_data(std::vector< std::vector< std::vector<TYPE> > >& data,
                                                    bool                                             split_fevalues) const
{
       if(split_fevalues)
       {
              std::vector< std::vector<TYPE> > local_data;

              unsigned int comp = 0;
              for(unsigned int b = 0; b < this->block_info->base_element.size(); ++b)
              {
                     const unsigned int       no_fe          = this->block_info->base_element[b];
                     const FEValuesBase<dim>& fe_values_base = this->fe(no_fe);

                     const unsigned int n_comp = fe_values_base.get_fe().n_components();
                     local_data.resize(n_comp,
                                       std::vector<TYPE>(fe_values_base.n_quadrature_points));

                     for(unsigned int i = 0; i < this->global_data->n_vectors(); ++i)
                     {
                       const FEVector& src = this->global_data->vector(i);
                       fill_data(fe_values_base,
                                 src,
                                 this->indices,
                                 this->block_info->local.block_start(b),
                                 fe_values_base.dofs_per_cell,
                                 local_data);

                       for(unsigned int c = 0; c < local_data.size(); ++c)
                         for(unsigned int k = 0; k < local_data[c].size(); ++k)
                           data[i][comp+c][k] = local_data[c][k];
                     }

                     comp += n_comp;
              }
       }
       else
       {
              for(unsigned int i = 0; i < this->global_data->n_vectors(); ++i)
              {
                     const FEVector&          src            = this->global_data->vector(i);
                     const FEValuesBase<dim>& fe_values_base = this->fe();

                     fill_data(fe_values_base,
                               src,
                               this->indices,
                               0,
                               fe_values_base.get_fe().dofs_per_cell,
                               data[i]);
              }
       }
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
inline
const FEVALUESBASE&
IntegrationInfo<dim, FEVALUESBASE>::get_fe_val_unsplit() const
{
       return *fe_val_unsplit;
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

template<int dim, typename FEVALUESBASE>
inline
void
IntegrationInfo<dim, FEVALUESBASE>::clear()
{
       fevalv.clear();
       fevalv.resize(0);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

} // ApplicationCore

} // FuelCell

#endif