librheolef-dev 6.5-1+b1 / 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
//
#include "rheolef/field.h"
#include "rheolef/field_nonlinear_expr.h"
namespace rheolef {

// --------------------------------------------------------------------------
// interpolate a general nonlinear expression
// --------------------------------------------------------------------------
// function template partial specialization is not allowed --> class-function
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 field_nonlinear_expr<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>();
}};
// scalar-valued result field
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,T,mpl::true_> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&           Xh,
    const field_nonlinear_expr<Expr>& expr) const
{
  typedef typename field_basic<T,M>::size_type size_type;
  bool is_homogeneous = expr.initialize (Xh);
  expr.template valued_check<T>();
  trace_macro ("is_homogeneous: " << is_homogeneous);
  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();
  check_macro (Xh.valued_tag() == space_constant::scalar,
	"interpolate: invalid "<<Xh.valued()<<"-valued space and scalar function");
  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;
    expr.evaluate (K, value);
    Xh.dis_idof   (K, dis_idof);
    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])) continue;
      size_type idof = dis_idof [loc_idof] - first_dis_idof;
      random_iter_uh [idof] = value [loc_idof];
    }
  }
  return uh;
}};
// vector-valued case
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,point_basic<T>,mpl::true_> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&           Xh,
    const field_nonlinear_expr<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: invalid "<<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;
    expr.evaluate (K, value);
    Xh.dis_idof   (K, dis_idof);
    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])) continue;
        size_type idof = dis_idof [loc_idof] - first_dis_idof;
        random_iter_uh [idof] = value [loc_comp_idof][i_comp];
      }
    }
  }
  return uh;
}};
// tensor-valued case
template<class T, class M, class Expr>
struct interpolate_internal_check<T,M,Expr,tensor_basic<T>,mpl::true_> {
  field_basic<T,M>
  operator() (
    const space_basic<T,M>&           Xh,
    const field_nonlinear_expr<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;
        if (! ownership.is_owned (dis_idof [loc_idof])) continue;
        size_type idof = dis_idof [loc_idof] - first_dis_idof;
        std::pair<size_type,size_type> ij
          = space_constant::tensor_subscript (space_constant::tensor, sys_coord, ij_comp);
        random_iter_uh [idof] = value [loc_comp_idof](ij.first, ij.second);
      }
    }
  }
  return uh;
}};
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 field_nonlinear_expr<Expr>& expr) const 
{
  switch (expr.valued_tag()) {
    case space_constant::scalar: {
	interpolate_internal_check<T,M,Expr,T,mpl::true_> eval;
	return eval (Xh, expr);
    }
    case space_constant::vector: {
	interpolate_internal_check<T,M,Expr,point_basic<T>,mpl::true_> eval;
	return eval (Xh, expr);
    }
    case space_constant::tensor:
    case space_constant::unsymmetric_tensor: {
	interpolate_internal_check<T,M,Expr,tensor_basic<T>,mpl::true_> 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>();
  }
}};
template<class T, class M, class Expr, class Result>
field_basic<T,M>
interpolate_internal (
    const space_basic<T,M>&           Xh,
    const field_nonlinear_expr<Expr>& expr)
{
  interpolate_internal_check<T,M,Expr,Result> eval;
  return eval (Xh,expr);
}
// undeterminated-valued case
template<class T, class M, class Expr>
field_basic<T,M>
interpolate (const space_basic<T,M>& Xh, const field_nonlinear_expr<Expr>& expr)
{
  typedef typename Expr::value_type result_t;
  return interpolate_internal<T,M,Expr,result_t> (Xh, expr);
}
// --------------------------------------------------------------------------
// implementation details of the interpolate() function
// --------------------------------------------------------------------------
namespace details {

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);
  }
  return u;
}
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];
    }
  }
  return u;
}
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);
    }
  }
  return u;
}

} // namespace details

// --------------------------------------------------------------------------
// interface of 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:
// TODO: un peu general... utiliser des specialisation des classes-fonctions:
//	function<Float(const point&)>
//	function<point(const point&)>
//	function<tensor(const point&)>
template <class T, class M, class Function>
inline
field_basic<T,M>
interpolate (const space_basic<T,M>& Xh, const Function& f)
//>interpolate:
{
  typedef typename Function::result_type  result_t;
  return details::interpolate_tag (Xh, f, result_t());
}
// specialization for scalar-valued functions:
template <class T, class M>
inline
field_basic<T,M>
interpolate (const space_basic<T,M>& Xh, T (*f)(const point_basic<T>&))
{
  return details::interpolate_tag (Xh, f, T());
}
// specialization for vector-valued functions:
template <class T, class M>
inline
field_basic<T,M>
interpolate (const space_basic<T,M>& Xh, point_basic<T> (*f)(const point_basic<T>&))
{
  return details::interpolate_tag (Xh, f, point_basic<T>());
}
// specialization for re-interpoltion of fields (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);

}// namespace rheolef
#endif // _RHEOLEF_INTERPOLATE_H