init_expr_quadrature.h

Go to the documentation of this file.

#ifndef _RHEOLEF_INIT_EXPR_QUADRATURE_H
#define _RHEOLEF_INIT_EXPR_QUADRATURE_H
//
// set iopt.order when not yet set, from expr polynomial degree
//
// author: Pierre.Saramito@imag.fr
//
// date: 18 march 2018 
//
#include "rheolef/form_expr_variational.h"
 
namespace rheolef { namespace details {
 
template<class T, class M>
integrate_option
expr_quadrature_init_iopt (
    const geo_basic<T,M>&   omega_K,
    const space_basic<T,M>& X,
    size_t                  n_derivative,
    const integrate_option& iopt)
{
  using size_type = quadrature_option::size_type;
  integrate_option new_iopt = iopt;
  if (iopt.get_order() == std::numeric_limits<quadrature_option::size_type>::max()) {
    size_type k = X.get_constitution().degree_max();
    size_type quad_order = 2*k + 1;
    new_iopt.set_order (quad_order);
  }
  return new_iopt;
}
template<class T, class M>
integrate_option
expr_quadrature_init_iopt (
    const geo_basic<T,M>&   omega_K,
    const space_basic<T,M>& X,
    const space_basic<T,M>& Y,
    size_t                  n_derivative,
    const integrate_option& iopt)
{
  using size_type = quadrature_option::size_type;
  integrate_option new_iopt = iopt;
  if (new_iopt.get_order() == std::numeric_limits<size_type>::max()) {
    size_type k1 = X.get_constitution().degree_max();
    size_type k2 = Y.get_constitution().degree_max();
    size_type quad_order = k1 + k2 + 1;
    if (omega_K.get_background_geo().sizes().ownership_by_variant[reference_element::q].dis_size() != 0 ||
        omega_K.get_background_geo().sizes().ownership_by_variant[reference_element::P].dis_size() != 0 ||
        omega_K.get_background_geo().sizes().ownership_by_variant[reference_element::H].dis_size() != 0) {
          // integrate exactly ??
          quad_order += 2;
    }
    if (omega_K.coordinate_system() != space_constant::cartesian) quad_order++; // multiplies by a 'r' weight
    if (quad_order >= n_derivative) quad_order -= n_derivative;
    new_iopt.set_order (quad_order);
  }
  return new_iopt;
}
template<class T>
quadrature<T>
expr_quadrature_init_quad (const integrate_option& iopt)
{
  quadrature<T> quad;
  quad.set_order (iopt.get_order());
  quad.set_family(iopt.get_family());
  return quad;
}
} // namespace details
} // namespace rheolef
#endif // _RHEOLEF_INIT_EXPR_QUADRATURE_H