librheolef-dev 6.7-6 / usr / include / rheolef / interpolate.h

#ifndef _RHEOLEF_INTERPOLATE_H
#define _RHEOLEF_INTERPOLATE_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
// 
// interpolation of functions and nonlinear expressions
// also re-interpolation of fields and linear expressions
//
// author: Pierre.Saramito@imag.fr
//
// date: 15 september 2015
//
// Notes; use template expressions and SFINAE techiques
//   The interpolation operator is required, as in
//   1/uh and uh*vh that do not belong to Xh when uh, vh in Xh
//
// SUMMARY:
// 1. implementation
// 1.1. scalar-valued result field
// 1.2. vector-valued case
// 1.3. tensor-valued case
// 1.4. undeterminated-valued case: determined at run-time
// 1.5. interface of the internal interpolate() function
// 2. the interpolate() function
// 2.1. nonlinear expression and function/functor
// 2.2. re-interpolation of fields and linear field expressions

#include "rheolef/field.h"
#include "rheolef/field_expr_v2_nonlinear.h"
namespace rheolef {

namespace details {
// --------------------------------------------------------------------------
// 1. implementation: general nonlinear expr
// --------------------------------------------------------------------------
// notes:
// * function template partial specialization is not allowed --> use class-function
// * interpolation on boundary domain spaces (geo_domain) could generate
//   some communications (see interpolate_dom[234]_tst.cc for tests)
// * interpolation of functions and functor are be performed in all cases
//   without communications, see below for this implementation.
//   here we suppose a general nonlinear expression that is interpolated
//   by using a loop on mesh elements
//
template<class T, class M, class Expr, class Result, class Status = typename details::is_equal<Result,typename Expr::value_type>::type>
struct interpolate_internal_check {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&  Xh,
    const Expr&              expr) const
{
  trace_macro ("Expr="<<pretty_typename_macro(Expr));
  trace_macro ("Result="<<typename_macro(Result));
  trace_macro ("Status="<<typename_macro(Status));
  trace_macro ("Expr::value_type="<<typename_macro(typename Expr::value_type));
  fatal_macro ("invalid type resolution");
  return field_basic<T,M>();
}};
// 1.1. scalar-valued result field
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,T,std::true_type> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&  Xh,
    const Expr&              expr) const
{
  check_macro (Xh.valued_tag() == space_constant::scalar,
	"interpolate: invalid "<<Xh.valued()<<"-valued space and scalar function");
  typedef typename field_basic<T,M>::size_type size_type;
  bool is_homogeneous = expr.initialize (Xh);
  expr.template valued_check<T>();
  const geo_basic<T,M>& omega = Xh.get_geo();
  field_basic<T,M> uh (Xh, std::numeric_limits<T>::max());
  typename field_basic<T,M>::iterator random_iter_uh = uh.begin_dof();
  std::vector<size_type> dis_idof;
  distributor ownership = Xh.ownership();
  size_type first_dis_idof = ownership.first_index();
  std::vector<T> value;
  for (typename geo_basic<T,M>::const_iterator
        iter_ie = omega.begin(),
        last_ie = omega.end(); iter_ie != last_ie; ++iter_ie) {
    const geo_element& K = *iter_ie;
    Xh.dis_idof   (K, dis_idof);
    expr.evaluate (K, value);
    for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; ++loc_idof) {
      if (ownership.is_owned (dis_idof [loc_idof])) {
        size_type idof = dis_idof [loc_idof] - first_dis_idof;
        random_iter_uh [idof] = value [loc_idof];
      } else {
	uh.dis_dof_entry (dis_idof [loc_idof]) = value [loc_idof];
      }
    }
  }
  uh.dis_dof_update();
  return uh;
}};
// 1.2. vector-valued case
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,point_basic<T>,std::true_type> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&  Xh,
    const Expr&              expr) const
{
  typedef typename field_basic<T,M>::size_type size_type;
  bool is_homogeneous = expr.initialize (Xh);
  expr.template valued_check<point_basic<T> >();
  const geo_basic<T,M>& omega = Xh.get_geo();
  field_basic<T,M> uh (Xh, std::numeric_limits<T>::max());
  std::vector<size_type> dis_idof;
  distributor ownership = Xh.ownership();
  size_type first_dis_idof = ownership.first_index();
  typename field_basic<T,M>::iterator random_iter_uh = uh.begin_dof();
  check_macro (Xh.valued_tag() == space_constant::vector,
	"interpolate: incompatible "<<Xh.valued()<<"-valued space and vector function");
  std::vector<point_basic<T> > value;
  size_type dim = Xh.get_geo().dimension();
  size_type n_comp = dim;
  for (typename geo_basic<T,M>::const_iterator
        iter_ie = omega.begin(),
        last_ie = omega.end(); iter_ie != last_ie; ++iter_ie) {
    const geo_element& K = *iter_ie;
    Xh.dis_idof   (K, dis_idof);
    expr.evaluate (K, value);
    size_type dis_ndof = dis_idof.size();
    size_type loc_comp_ndof = dis_ndof/n_comp;
    for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
      for (size_type i_comp = 0; i_comp < n_comp; ++i_comp) {
        size_type loc_idof = loc_comp_idof + i_comp*loc_comp_ndof;
        if (ownership.is_owned (dis_idof [loc_idof])) {
          size_type idof = dis_idof [loc_idof] - first_dis_idof;
          random_iter_uh [idof] = value [loc_comp_idof][i_comp];
        } else {
	  uh.dis_dof_entry (dis_idof [loc_idof]) = value [loc_comp_idof][i_comp];
        }
      }
    }
  }
  uh.dis_dof_update();
  return uh;
}};
// 1.3. tensor-valued case
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,tensor_basic<T>,std::true_type> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>& Xh,
    const Expr&             expr) const
{
  typedef typename field_basic<T,M>::size_type size_type;
  bool is_homogeneous = expr.initialize (Xh);
  expr.template valued_check<tensor_basic<T> >();
  const geo_basic<T,M>& omega = Xh.get_geo();
  field_basic<T,M> uh (Xh, std::numeric_limits<T>::max());
  std::vector<size_type> dis_idof;
  distributor ownership = Xh.ownership();
  size_type first_dis_idof = ownership.first_index();
  typename field_basic<T,M>::iterator random_iter_uh = uh.begin_dof();
  check_macro (Xh.valued_tag() == space_constant::tensor ||
               Xh.valued_tag() == space_constant::unsymmetric_tensor,
	"interpolate: invalid "<<Xh.valued()<<"-valued space and vector function");
  std::vector<tensor_basic<T> > value;
  size_type dim = Xh.get_geo().dimension();
  space_constant::coordinate_type sys_coord = uh.get_geo().coordinate_system();
  size_type n_comp = space_constant::n_component (space_constant::tensor, dim, sys_coord);
  for (typename geo_basic<T,M>::const_iterator
        iter_ie = omega.begin(),
        last_ie = omega.end(); iter_ie != last_ie; ++iter_ie) {
    const geo_element& K = *iter_ie;
    expr.evaluate (K, value);
    Xh.dis_idof   (K, dis_idof);
    size_type loc_ndof = dis_idof.size();
    size_type loc_comp_ndof = loc_ndof/n_comp;
    for (size_type loc_comp_idof = 0; loc_comp_idof < loc_comp_ndof; ++loc_comp_idof) {
      for (size_type ij_comp = 0; ij_comp < n_comp; ++ij_comp) {
        size_type loc_idof = loc_comp_idof + ij_comp*loc_comp_ndof;
        std::pair<size_type,size_type> ij
          = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
        if (ownership.is_owned (dis_idof [loc_idof])) {
          size_type idof = dis_idof [loc_idof] - first_dis_idof;
          random_iter_uh [idof] = value [loc_comp_idof](ij.first, ij.second);
        } else {
	  uh.dis_dof_entry (dis_idof [loc_idof]) = value [loc_comp_idof](ij.first, ij.second);
        }
      }
    }
  }
  uh.dis_dof_update();
  return uh;
}};
// 1.4. undeterminated-valued case: determined at run-time
template<class T, class M, class Expr, class Status>
struct interpolate_internal_check<T,M,Expr,undeterminated_basic<T>,Status> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&  Xh,
    const Expr&              expr) const 
{
  switch (expr.valued_tag()) {
    case space_constant::scalar: {
	interpolate_internal_check<T,M,Expr,T,std::true_type> eval;
	return eval (Xh, expr);
    }
    case space_constant::vector: {
	interpolate_internal_check<T,M,Expr,point_basic<T>,std::true_type> eval;
	return eval (Xh, expr);
    }
    case space_constant::tensor:
    case space_constant::unsymmetric_tensor: {
	interpolate_internal_check<T,M,Expr,tensor_basic<T>,std::true_type> eval;
	return eval (Xh, expr);
    }
    default:
        warning_macro ("Expr="<<pretty_typename_macro(Expr));
        warning_macro ("Status="<<typename_macro(Status));
        error_macro ("unexpected `"
	<< space_constant::valued_name (expr.valued_tag())
        << "' valued expression");
        return field_basic<T,M>();
  }
}};
// 1.5. interface of the internal interpolate() function
template<class T, class M, class Expr, class Result>
field_basic<T,M>
interpolate_internal (
    const space_basic<T,M>& Xh,
    const Expr&             expr)
{
  interpolate_internal_check<T,M,Expr,Result> eval;
  return eval (Xh,expr);
}

} // namespace details
// --------------------------------------------------------------------------
// 2. implementation: function & functors
// --------------------------------------------------------------------------
// notes:
// * interpolation of functions and functor are be performed in all cases
//   without communications: the loop bases on Lagrange nodes directly.
//
namespace details {

// 2.1. scalar valued
template <class T, class M, class Function>
field_basic<T,M>
interpolate_tag (const space_basic<T,M>& Xh, const Function& f, const T&)
{
  check_macro (Xh.valued_tag() == space_constant::scalar,
	"interpolate: invalid "<<Xh.valued()<<"-valued " <<Xh.stamp()
	<< " space and scalar function");
  typedef typename space_basic<T,M>::size_type size_type;
  field_basic<T,M> u (Xh, std::numeric_limits<T>::max());
  for (size_type idof = 0, ndof = Xh.ndof(); idof < ndof; idof++) {
    u.dof (idof) = Xh.momentum (f, idof);
  }
  u.dis_dof_update();
  return u;
}
// 2.2. vector valued
template <class T, class M, class Function>
field_basic<T,M>
interpolate_tag (const space_basic<T,M>& Xh, const Function& f, const point_basic<T>&)
{
  check_macro (Xh.valued_tag() == space_constant::vector,
	"interpolate: invalid "<<Xh.valued()<<"-valued space and vector function");
  typedef typename space_basic<T,M>::size_type size_type;
  field_basic<T,M> u (Xh, std::numeric_limits<T>::max());
  size_type n_comp = Xh.get_geo().dimension();
  point_basic<T> value;
  for (size_type comp_idof = 0, comp_ndof = Xh.ndof()/n_comp; comp_idof < comp_ndof; comp_idof++) {
    value = Xh.vector_momentum (f, comp_idof);
    for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
      size_type idof = comp_idof + i_comp*comp_ndof;
      u.dof (idof) = value [i_comp];
    }
  }
  u.dis_dof_update();
  return u;
}
// 2.2. tensor valued
template <class T, class M, class Function>
field_basic<T,M>
interpolate_tag (const space_basic<T,M>& Xh, const Function& f, const tensor_basic<T>&)
{
  check_macro (Xh.valued_tag() == space_constant::tensor,
	"interpolate: invalid "<<Xh.valued()<<"-valued space and tensor function");
  typedef typename space_basic<T,M>::size_type size_type;
  field_basic<T,M> u (Xh, std::numeric_limits<T>::max());
  size_type d = Xh.get_geo().dimension();
  space_constant::coordinate_type sys_coord = Xh.get_geo().coordinate_system();
  size_type n_comp = space_constant::n_component (space_constant::tensor, d, sys_coord);
  tensor_basic<T> value;
  for (size_type comp_idof = 0, comp_ndof = Xh.ndof()/n_comp; comp_idof < comp_ndof; comp_idof++) {
    value = Xh.tensor_momentum (f, comp_idof);
    for (size_type ij_comp = 0; ij_comp < n_comp; ij_comp++) {
      size_type idof = comp_idof + ij_comp*comp_ndof;
      std::pair<size_type,size_type> ij
          = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
      u.dof (idof) = value (ij.first, ij.second);
    }
  }
  u.dis_dof_update();
  return u;
}

} // namespace details

// --------------------------------------------------------------------------
// 2. the interpolate() function
// --------------------------------------------------------------------------

/*Class:interpolate
NAME: @code{interpolate} - Lagrange interpolation of a function
@findex  interpolate
@clindex space
@clindex field

DESCRIPTION:
 The function @code{interpolation} implements the
 Lagrange interpolation of a function or a class-function.
SYNOPSYS:
 template <class Function>
 field interpolate (const space& Xh, const Function& f);
EXAMPLE:
@noindent
 The following code compute the Lagrange interpolation
 @code{pi_h_u} of u(x).
@example
  Float u(const point& x);
  ...
  geo omega("square");
  space Xh (omega, "P1");
  field pi_h_u = interpolate (Xh, u);
@end example
ADVANCED EXAMPLE:
  It is possible the replace the function @code{u} 
  by a variable of the @code{field} type that represents
  a picewise polynomial function: this invocation allows
  the reinterpolation of a field on another mesh or with
  another approximation.
@example
  geo omega2 ("square2");
  space X2h (omega2, "P1");
  field uh2 = interpolate (X2h, pi_h_u);
@end example
End: */
//>interpolate:

// ----------------------------------------------------
// 2. public interface
// ----------------------------------------------------
// 2.1. general nonlinear expression
template<class T, class M, class Expr>
typename std::enable_if<
       details::is_field_expr_v2_nonlinear_arg<Expr>::value
  && ! details::is_field_expr_v2_linear_arg<Expr>::value
  && ! details::is_field_function<Expr>::value
 ,field_basic<T,M>
>::type
interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type wrap_t;
  typedef typename wrap_t::value_type result_t;
  return details::interpolate_internal<T,M,wrap_t,result_t> (Xh, wrap_t(expr));
}
// 2.2. function & functor
template <class T, class M, class Expr>
inline
typename std::enable_if<
  details::is_field_function<Expr>::value
 ,field_basic<T,M>
>::type
interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
  typedef typename details::function_traits<Expr>::result_type result_t;
  return details::interpolate_tag (Xh, expr, result_t());
}

// 2.3. re-interpolation of fields and linear field expressions
//      for change of mesh, of approx, ect
template<class T, class M>
field_basic<T,M>
interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h);

template <class T, class M, class Expr>
inline
typename std::enable_if<
     details::is_field_expr_v2_linear_arg<Expr>::value
  && ! details::is_field<Expr>::value
 ,field_basic<T,M>
>::type
interpolate (const space_basic<T,M>& Xh, const Expr& expr)
{
  return interpolate (Xh, field_basic<T,M>(expr));
}

}// namespace rheolef
#endif // _RHEOLEF_INTERPOLATE_H