librheolef-dev 5.93-2 / usr / include / rheolef / interpolate.h

#ifndef _RHEO_INTERPOLATE_H
#define _RHEO_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 with a class-function: scalar- or vector-valuated
//
#include "rheolef/tensor.h"
#include "rheolef/fem_helper.h"
namespace rheolef { 

template <class Function>
field __interpolate_tag (const space& V, Function f, const Float&) {
  typedef field::size_type size_type;
  const geo&  g = V.get_geo() ;  
  const basis& b = V.get_basis() ;  
  field u;
  if (b.family()!=reference_element::Lagrange) u=field(V,0);
  else u=field(V);
  check_macro (V.n_component() == 1, "interpolate: expect scalar space, get "
	<< V.n_component() << "D vector-valued space.");
  std::vector<bool> marked (V.size(), false);
  vec<Float>::iterator xu = u.u.begin();
  vec<Float>::iterator xb = u.b.begin();
  bool Lagrange=(b.family()==reference_element::Lagrange);
  for (geo::const_iterator iter_K = g.begin(); iter_K != g.end(); iter_K++) {
      const geo_element& K = *iter_K;
      size_type nb_dof = b.size(K.type()) ;
      if (!Lagrange) nb_dof=b.size(K.type(), reference_element::Lagrange);
      tiny_vector<space::size_type> dof(nb_dof);
      if (Lagrange) V.set_dof(K, dof, 0) ;
      else V.set_dof(K, dof, 0, reference_element::Lagrange) ;
      for(size_type i = 0; i < nb_dof; i++) {
	  size_type i_dof = dof[i];
	  if (marked [i_dof]) continue; else marked [i_dof] = true;
	  point x = V.x_dof (K, i);
	  Float value = f(x);
	  size_type idx = V.index(i_dof);
	  if (V.is_blocked(i_dof))
	    xb [idx] = value;
	  else
	    xu [idx] = value;
      }
  }
  return u;
}
template <class Function>
field __interpolate_tag (const space& V, Function f, const point&) {
  typedef field::size_type size_type;
  field u(V);
  const geo&  g = V.get_geo() ;  
  const basis& b = V.get_basis() ;  
  size_type n_comp          = V.n_component();
  check_macro (n_comp > 0 && n_comp <= 3, "interpolate: expect a vector-valued functionnal space, get a "
	<< n_comp << "D vector-valued one.");
  check_macro (V.get_valued() == "vector", "interpolate: expect a vector-valued functionnal space, get a "
	<< V.get_valued() << " one");
  std::vector<bool> marked (V.size(), false);
  vec<Float>::iterator xu = u.u.begin();
  vec<Float>::iterator xb = u.b.begin();
  tiny_vector<space::size_type> dof [n_comp];
  for (geo::const_iterator iter_K = g.begin(); iter_K != g.end(); iter_K++) {
      const geo_element& K = *iter_K;
      size_type nb_dof = b.size(K.type()) ;
      for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
          dof[i_comp].resize(nb_dof);
          V.set_dof(K, dof[i_comp], i_comp);
      }
      for(size_type i = 0; i < nb_dof; i++) {
	  size_type i_dof_0 = dof[0][i];
	  if (marked [i_dof_0]) continue; else marked [i_dof_0] = true;
	  point x = V.x_dof (K, i);
	  point value = f(x);
	  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
	    size_type i_dof = dof[i_comp][i];
	    size_type idx = V.index(i_dof);
	    if (V.is_blocked(i_dof)) {
	      xb [idx] = value [i_comp];
	    } else {
	      xu [idx] = value [i_comp];
            }
          }
      }
  }
  return u;
}
template <class Function>
field __interpolate_tag (const space& V, Function f, const tensor&) {
  typedef field::size_type size_type;
  field u(V);
  const geo&  g = V.get_geo() ;  
  const basis& b = V.get_basis() ;  
  size_type n_comp          = V.n_component();
  check_macro (n_comp > 0 && n_comp <= 9,
	"interpolate: expect a tensor-valued functionnal space, get a "
	<< n_comp << "D vector-valued one.");
  std::string valued_name = V.get_valued();
  check_macro (valued_name == "tensor" || valued_name == "unsymmetric_tensor",
	"interpolate: expect a tensdor-valued functionnal space, get a "
	<< valued_name << " one");
  std::string sys_coord_name = V.coordinate_system();
  fem_helper::coordinate_type sys_coord = fem_helper::coordinate_system(sys_coord_name);
  fem_helper::valued_field_type valued = fem_helper::valued_field (valued_name);
  std::vector<bool> marked (V.size(), false);
  vec<Float>::iterator xu = u.u.begin();
  vec<Float>::iterator xb = u.b.begin();
  tiny_vector<space::size_type> dof [n_comp];
  fem_helper::pair_size_type ij_comp;
  for (geo::const_iterator iter_K = g.begin(); iter_K != g.end(); iter_K++) {
      const geo_element& K = *iter_K;
      size_type nb_dof = b.size(K.type()) ;
      for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
          dof[i_comp].resize(nb_dof);
          V.set_dof(K, dof[i_comp], i_comp);
      }
      for(size_type i = 0; i < nb_dof; i++) {
	  size_type i_dof_0 = dof[0][i];
	  if (marked [i_dof_0]) continue; else marked [i_dof_0] = true;
	  point x = V.x_dof (K, i);
	  tensor value = f(x);
	  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
	    size_type i_dof = dof[i_comp][i];
	    size_type idx = V.index(i_dof);
	    ij_comp = fem_helper::tensor_subscript (valued, sys_coord, i_comp);
	    if (V.is_blocked(i_dof)) {
	      xb [idx] = value (ij_comp.first, ij_comp.second);
	    } else {
	      xu [idx] = value (ij_comp.first, ij_comp.second);
	    }
          }
      }
  }
  return u;
}
template <class Func>
inline
field
interpolate (const space& V, Func f) {
    typedef typename Func::result_type  result_type ;
    return __interpolate_tag (V, f, result_type());
}
// full specialization for functions:
template <>
inline
field
interpolate (const space& V, Float (*f)(const point&)) {
    return __interpolate_tag (V, f, Float());
}
template <>
inline
field
interpolate (const space& V, point (*f)(const point&)) {
    return __interpolate_tag (V, f, point());
}
template <>
inline
field
interpolate (const space& V, tensor (*f)(const point&)) {
    return __interpolate_tag (V, f, tensor());
}

// spline interpolation
void calculate_spline(field& ih);
inline
field
spline (const space& V, const field& uh) {
    check_macro(V.get_basis().family()==reference_element::Hermite,
    	"Only Hermite spaces allowed for spline interpoltion");
    // interpolate Lagrange dofs
    const space& W=uh.get_space();
    space S(V.get_geo(), V.get_approx(), V.get_valued());
    V.freeze();
    if (V.n_blocked()!=0) error_macro("Boundary conditions not yet supported for spline interpolation");
    S.block_Lagrange();
    field ih;
    if ((S.n_component() == 1) &&
  	 (W.n_component() == 1) &&
  	 (S.get_geo().familyname()==W.get_geo().familyname()) &&
	 (S.get_geo().version()==W.get_geo().version()) &&
  	 (S.get_approx()==W.get_approx()) )
	   ih = field(S, uh);
    else
      // interpolate Lagrange dofs
      field ih=interpolate (S, uh);
    // calculate the Hermite ones and transfer the field
    calculate_spline(ih);
    // back to original space (without blocked Lagrange dofs)
    return field(V, ih);
}
template <class T>
inline
field
spline (const space& V, T (*f)(const point&)) {
    check_macro(V.get_basis().family()==reference_element::Hermite,
    	"Only Hermite spaces allowed for spline interpoltion");
    // interpolate Lagrange dofs
    space S(V.get_geo(), V.get_approx(), V.get_valued());
    V.freeze();
    if (V.n_blocked()!=0) error_macro("Boundary conditions not yet supported for spline interpolation");
    S.block_Lagrange();
    // interpolate Lagrange dofs
    field ih=interpolate (S, f);
    // calculate the Hermite ones and transfer the field
    calculate_spline(ih);
    // back to original space (without blocked Lagrange dofs)
    return field(V, ih);
}
/*
template <>
inline
field
spline (const space& V, point (*f)(const point&)) {
    return __spline_tag (V, f, point());
}
template <>
inline
field
spline (const space& V, tensor (*f)(const point&)) {
    return __spline_tag (V, f, tensor());
}
*/
}// namespace rheolef
#endif // _RHEO_INTERPOLATE_H