fe_values_cache.h

Go to the documentation of this file.

//-----------------------------------------------------------
//
//    Copyright (C) 2015 by the deal2lkit authors
//
//    This file is part of the deal2lkit library.
//
//    The deal2lkit 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 deal2lkit distribution.
//
//-----------------------------------------------------------

#ifndef _d2k_fe_values_cache_h_
#define _d2k_fe_values_cache_h_

#include <deal2lkit/config.h>
#include <deal.II/fe/fe_values.h>
#include <deal2lkit/any_data.h>
#include <deal2lkit/dof_utilities.h>
#include <deal2lkit/utilities.h>

D2K_NAMESPACE_OPEN

template <int dim, int spacedim=dim>
class FEValuesCache
{
public:
  FEValuesCache     (const Mapping<dim, spacedim>         &mapping,
                     const FiniteElement<dim, spacedim>   &fe,
                     const Quadrature<dim>                &quadrature,
                     const UpdateFlags                    &update_flags,
                     const Quadrature<dim-1>              &face_quadrature,
                     const UpdateFlags                    &face_update_flags):
    fe_values         (mapping, fe, quadrature, update_flags),
    fe_face_values    (mapping, fe, face_quadrature, face_update_flags),
    fe_subface_values (mapping, fe, face_quadrature, face_update_flags),
    local_dof_indices(fe.dofs_per_cell)

  {};

  FEValuesCache     (const FEValuesCache<dim,spacedim> &scratch):
    cache   (scratch.cache),
    fe_values ( scratch.fe_values.get_mapping(),
                scratch.fe_values.get_fe(),
                scratch.fe_values.get_quadrature(),
                scratch.fe_values.get_update_flags()),
    fe_face_values ( scratch.fe_values.get_mapping(),
                     scratch.fe_values.get_fe(),
                     scratch.fe_face_values.get_quadrature(),
                     scratch.fe_face_values.get_update_flags()),
    fe_subface_values ( scratch.fe_values.get_mapping(),
                        scratch.fe_values.get_fe(),
                        scratch.fe_subface_values.get_quadrature(),
                        scratch.fe_subface_values.get_update_flags()),
    local_dof_indices(scratch.fe_values.get_fe().dofs_per_cell)
  {};


  void reinit(const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell)
  {
    fe_values.reinit(cell);
    cell->get_dof_indices(local_dof_indices);
    cache.template add_ref<FEValuesBase<dim,spacedim> >(fe_values, "FEValuesBase");
  };


  void reinit(const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
              const unsigned int face_no)
  {
    fe_face_values.reinit(cell, face_no);
    cell->get_dof_indices(local_dof_indices);
    cache.template add_ref<FEValuesBase<dim,spacedim> >(fe_face_values, "FEValuesBase");
  };


  void reinit(const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
              const unsigned int face_no, const unsigned int subface_no)
  {
    fe_subface_values.reinit(cell, face_no, subface_no);
    cell->get_dof_indices(local_dof_indices);
    cache.template add_ref<FEValuesBase<dim,spacedim> >(fe_subface_values, "FEValuesBase");
  };

  AnyData &get_cache()
  {
    return cache;
  };

  const FEValuesBase<dim,spacedim> &get_current_fe_values()
  {
    Assert(cache.have("FEValuesBase"),
           ExcMessage("You have to initialize the cache using one of the "
                      "reinit function first!"));
    FEValuesBase<dim,spacedim> &fev =
      cache.template get<FEValuesBase<dim,spacedim> >("FEValuesBase");
    return fev;
  }


  template<typename Number>
  std::vector<Number> &
  get_current_independent_local_dofs(const std::string &prefix, Number dummy)
  {

    std::string dofs_name = prefix+"_independent_local_dofs_"+type(dummy);

    Assert(cache.have(dofs_name),
           ExcMessage("You did not call cache_local_solution_vector with the right types!"));

    return  cache.template get<std::vector<Number> >(dofs_name);
  }


  const std::vector<Point<spacedim> > &
  get_quadrature_points()
  {
    return get_current_fe_values().get_quadrature_points();
  }


  const std::vector<double > &
  get_JxW_values()
  {
    return get_current_fe_values().get_JxW_values();
  }


  template <class STREAM>
  void print_info(STREAM &os)
  {
    cache.print_info(os);
  }


  template<typename VEC, typename Number>
  void cache_local_solution_vector(const std::string &prefix, const VEC &input_vector, const Number dummy)
  {
    const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell;

    std::string name = prefix+"_independent_local_dofs_"+type(dummy);

    if (!cache.have(name))
      cache.add_copy(std::vector<Number>(n_dofs), name);

    std::vector<Number> &independent_local_dofs
      = cache.template get<std::vector<Number> >(name);
    DOFUtilities::extract_local_dofs(input_vector, local_dof_indices, independent_local_dofs);
  }




  template<typename Number>
  const std::vector <std::vector<Number> > &
  get_values(const std::string &prefix, const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;
    const unsigned int           n_components  = fev.get_fe().n_components();

    std::string name = prefix+"_all_values_q"+Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    if (!cache.have(name))
      {
        cache.add_copy(std::vector <std::vector<Number> >(n_q_points, std::vector<Number>(n_components)), name);
      }

    std::vector <std::vector<Number> > &ret = cache.template get<std::vector <std::vector<Number> > >(name);
    DOFUtilities::get_values(fev, independent_local_dofs, ret);
    return ret;
  }



  template<typename Number>
  const std::vector <Number> &
  get_values(const std::string &prefix,
             const std::string &additional_prefix,
             const FEValuesExtractors::Scalar &variable,
             const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_scalar_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector<Number> RetType;


    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_values(fev, independent_local_dofs, variable, ret);
    return ret;
  }



  template<typename Number>
  const std::vector <Tensor <1, spacedim, Number> > &
  get_values(const std::string &prefix,
             const std::string &additional_prefix,
             const FEValuesExtractors::Vector &vector_variable,
             const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_vector_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);


    // Now build the return type
    typedef typename std::vector<Tensor<1, spacedim, Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_values(fev, independent_local_dofs, vector_variable, ret);
    return ret;
  }



  template<typename Number>
  const std::vector <Number> &
  get_divergences(const std::string &prefix,
                  const std::string &additional_prefix,
                  const FEValuesExtractors::Vector &vector_variable,
                  const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_div_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);


    // Now build the return type
    typedef typename std::vector<Number> RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_divergences(fev, independent_local_dofs, vector_variable, ret);
    return ret;
  }



  template<typename Number>
  const std::vector <Tensor<1, spacedim, Number> > &
  get_gradients(const std::string &prefix,
                const std::string &additional_prefix,
                const FEValuesExtractors::Scalar &variable,
                const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_grad_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<1, spacedim, Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_gradients(fev, independent_local_dofs, variable, ret);
    return ret;
  }



  template<typename Number>
  const std::vector <Tensor<2, spacedim, Number> > &
  get_gradients(const std::string &prefix,
                const std::string &additional_prefix,
                const FEValuesExtractors::Vector &variable,
                const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_grad2_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<2, spacedim, Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_gradients(fev, independent_local_dofs, variable, ret);
    return ret;
  }

  template<typename Number>
  const std::vector <Tensor<1, (spacedim > 2 ? spacedim : 1), Number> > &
  get_curls(const std::string &prefix,
            const std::string &additional_prefix,
            const FEValuesExtractors::Vector &variable,
            const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_curl_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<1, (spacedim > 2 ? spacedim : 1), Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_curls(fev, independent_local_dofs, variable, ret);
    return ret;
  }

  template<typename Number>
  const std::vector <Number> &
  get_laplacians(const std::string &prefix,
                 const std::string &additional_prefix,
                 const FEValuesExtractors::Scalar &variable,
                 const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_laplacian_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Number> RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_laplacians(fev, independent_local_dofs, variable, ret);
    return ret;
  }


  template<typename Number>
  const std::vector <Tensor<1,spacedim,Number> > &
  get_laplacians(const std::string &prefix,
                 const std::string &additional_prefix,
                 const FEValuesExtractors::Vector &variable,
                 const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_laplacian2_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<1,spacedim,Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_laplacians(fev, independent_local_dofs, variable, ret);
    return ret;
  }

  template<typename Number>
  const std::vector <Tensor<2, spacedim, Number> > &
  get_deformation_gradients(const std::string &prefix,
                            const std::string &additional_prefix,
                            const FEValuesExtractors::Vector &variable,
                            const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_F_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<2, spacedim, Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_deformation_gradients(fev, independent_local_dofs, variable, ret);
    return ret;
  }



  template<typename Number>
  const std::vector <Tensor<2, spacedim, Number> > &
  get_symmetric_gradients(const std::string &prefix,
                          const std::string &additional_prefix,
                          const FEValuesExtractors::Vector &variable,
                          const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_sym_grad_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<2, spacedim, Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_symmetric_gradients(fev, independent_local_dofs, variable, ret);
    return ret;
  }

  template<typename Number>
  const std::vector <Tensor<2,spacedim,Number> > &
  get_hessians(const std::string &prefix,
               const std::string &additional_prefix,
               const FEValuesExtractors::Scalar &variable,
               const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_hessians_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<2,spacedim,Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_hessians(fev, independent_local_dofs, variable, ret);
    return ret;
  }

  template<typename Number>
  const std::vector <Tensor<3,spacedim,Number> > &
  get_hessians(const std::string &prefix,
               const std::string &additional_prefix,
               const FEValuesExtractors::Vector &variable,
               const Number dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_current_independent_local_dofs(prefix, dummy);

    const FEValuesBase<dim,spacedim> &fev = get_current_fe_values();

    const unsigned int           n_q_points    = fev.n_quadrature_points;

    std::string name = prefix+"_"+additional_prefix+"_hessians2_values_q"+
                       Utilities::int_to_string(n_q_points)+"_"+type(dummy);

    // Now build the return type
    typedef typename std::vector <Tensor<3,spacedim,Number> > RetType;

    if (!cache.have(name))
      cache.add_copy(RetType(n_q_points), name);

    RetType &ret = cache.template get<RetType>(name);
    DOFUtilities::get_hessians(fev, independent_local_dofs, variable, ret);
    return ret;
  }

private:

  AnyData                                           cache;
  FEValues<dim, spacedim>                           fe_values;
  FEFaceValues<dim, spacedim>                       fe_face_values;
  FESubfaceValues<dim, spacedim>                    fe_subface_values;
  std::vector<types::global_dof_index>  local_dof_indices;
};

D2K_NAMESPACE_CLOSE

#endif