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

# ifndef _RHEO_FIELD_H
# define _RHEO_FIELD_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
/// 
/// =========================================================================

namespace rheolef { 
struct field_on_domain;
struct field_component;
struct const_field_component;
}// namespace rheolef

#include "rheolef/vec.h"
#include "rheolef/space.h"
#include "rheolef/form_diag.h"
#include "rheolef/tensor.h"

namespace rheolef { 
//<field:  
/*Class:field
NAME:  @code{field} - piecewise polynomial finite element field
DESCRIPTION:       
  @noindent
  Store degrees of freedom associated to a mesh and
  a piecewise polynomial approximation, with respect
  to the numbering defined by the underlying @ref{space class}.

  @noindent
  This class contains two vectors, namely unknown and blocked
  degrees of freedoms, and the associated finite element space.
  Blocked and unknown degrees of freedom can be set by using
  domain name indexation:
  @example
        geo omega_h ("circle");
        space Vh (omega_h, "P1");
        Vh.block ("boundary");
        field uh (Vh);
        uh ["boundary"] = 0;
  @end example
  Interpolation of a function @code{u} in a field @code{uh} with respect to 
  the interpolation writes:
  @example
        Float u (const point&);
        uh = interpolate (Vh, u);
  @end example

EXAMPLE:
  Here is a complete example
  using the field class:
@example
        Float u (const point& x) @{ return x[0]*x[1]; @}
        main() @{
          geo omega_h ("square");
          space Vh (omega_h, "P1");
          field uh (Vh);
          uh = interpolate (Vh, u);
          cout << plotmtv << u;
        @}
@end example
FEATURES:
   Algebra, such as x+y, x*y, x/y, lambda*x, ...are supported
   
   Transformations applies to all values of a field:
@example
   field vh = compose(fabs, uh);
   field wh = compose(atan2, uh, vh);
@end example
   The composition supports also general unary and binary class-functions.

   Vector-valued and tensor-valued field support is yet partial only.
   This feature will be more documented in the future.

AUTHOR: 
    LMC-IMAG, 38041 Grenoble cedex 9, France
   | Pierre.Saramito@imag.fr
   | Jocelyn.Etienne@imag.fr
DATE:   2 july 1997  update: 22 may 2010
METHODS: @field
End:
*/


class field : public std::unary_function<point,Float> {
public :
// typedefs:

    typedef vec<Float>::size_type size_type;
    typedef geo::plot_options plot_options;

// allocator/deallocator:
  
    field ();
    explicit field (const space& V, const Float& init_value = std::numeric_limits<Float>::max());
    explicit field (const space& V, std::istream& is);
    explicit field (const std::string& filename);
    explicit field (const space& V, const std::string& filename);
    explicit field (const class form_diag& D);
             field (const field_component& ui);
             field (const const_field_component& ui);
	      //! copy field v into the space V with different blocked/unknown dofs.
	     field (const space& V, const field& v);

    // return the field on an optimized mesh
    friend field adapt (const class field& criteria, const Float& hcoef);
    friend field adapt (const class field& criteria, const adapt_option_type& = adapt_option_type());

// accessors:

    size_type  size() const;
    size_type  n_unknown() const;
    size_type  n_blocked() const;
    size_type  dimension() const;
    std::string coordinate_system() const;
    fem_helper::coordinate_type coordinate_system_type() const;
    std::string get_approx(size_type i_comp=0) const;
    const space&  get_space() const;
    const geo&    get_geo() const;
    Float min() const;
    Float max() const;
    Float max_abs() const;

// accessor: val := uh(x) and such

    Float  operator()      (const point& x, size_type i_comp = 0) const;
    point  vector_evaluate (const point& x) const;
    tensor tensor_evaluate (const point& x) const;

    void evaluate (const point& x, Float& result) const;
    void evaluate (const point& x, point& result) const;
    void evaluate (const point& x, tensor& result) const;

//! accessors & modifiers by components, or vector & tensor valued:
//!       Float  u = uh.at (i_dof);
//!       point  v = vh.vector_at (i_dof);
//!       tensor m = mh.tensor_at (i_dof);

    Float  at            (size_type i_dof) const;
    Float& at            (size_type i_dof);

          Float  at      (size_type i_dof, size_type i_comp) const;
    const Float& set_at  (size_type i_dof, size_type i_comp, const Float& value);
          void   incr_at (size_type i_dof, size_type i_comp, const Float& value);

          point  vector_at      (size_type i_dof) const;
    const point& set_vector_at  (size_type i_dof, const point& value);
          void incr_vector_at (size_type i_dof, const point& value);

          tensor  tensor_at      (size_type i_dof) const;
    const tensor& set_tensor_at  (size_type i_dof, const tensor& value);
          void    incr_tensor_at (size_type i_dof, const tensor& value);

// assignment: u = lambda; u = v[i];

    field operator = (Float lambda);
    field operator = (const field_component&);
    field operator = (const const_field_component&);

// computed assignment: u op= lambda;

    field operator *= (Float lambda);
    field operator /= (Float lambda);
    field operator += (const field& f);
    field operator -= (const field& f);

// assignment to a domain:  u [dom] := value;

    field_on_domain operator [] (const domain& dom);
    field_on_domain operator [] (const std::string& dom_name);

// multi-field accessors: field[i], field(i,j), field(i,j,k,l)

    size_type n_component() const;
    field_component operator [] (size_type i_comp);
    field_component operator () (size_type i_comp, size_type j_comp);
    field_component operator () (size_type i, size_type j, size_type k, size_type l);
    const_field_component operator [] (size_type i_comp) const;
    const_field_component operator () (size_type i_comp, size_type j_comp) const;
    const_field_component operator () (size_type i, size_type j, size_type k, size_type l) const;

    std::string get_valued() const; // "scalar", "vector", "tensor", "tensor4"...
    fem_helper::valued_field_type get_valued_type() const;

// input/output:

    friend std::ostream& operator << (std::ostream& s, const field& x);
    friend std::istream& operator >> (std::istream& s, field& x);
    void save (std::string basename = std::string(), std::string username = std::string()) const;
    void put_vtk (std::ostream& s,
	const std::string& name = "scalars", bool put_header = true) const;
    std::ostream& put_vtk_vector (std::ostream& vtk,
    	const std::string& name = "vectors", bool put_header = true) const;
    std::ostream& put_vtk_tensor (std::ostream& vtk,
    	const std::string& name = "tensors", bool put_header = true) const;
    void check () const;
    int gnuplot2d_velocity (const std::string& basename, plot_options& opt) const;
 
// specific:

    	//! Transforms a P_k continuous on subdomains to a P_k discontinuous (currently k=1 only)
    field piecewise_to_discontinuous() const;

// implementation:
    	//! change to use tangential-normal representation locally 
    void set_locked_boundaries ();   
        //! change to use cartesian/cylindrical representation everywhere
    void unset_locked_boundaries (); 
    bool locked_boundaries_set () const { return _locked_boundaries_set; }

protected:

    friend std::istream& load_scalar_field (std::istream&, field&);
    friend std::istream& load_multi_field  (std::istream&, field&);

    std::ostream& put (std::ostream& s) const;
    std::ostream& put_scalar (std::ostream& s) const;
    std::ostream& put_multi (std::ostream& s) const;
    int plotmtv2d_P2 (const std::string& basename, bool execute, bool clean, bool verbose, 
		      bool grid = false, bool domains=false,
		      bool fill = false, bool elevation = false,
		      bool label = false,
		      size_type isovalue_table_zero_index = std::numeric_limits<size_type>::max()) const;
    int plotmtv2d_P1d(const std::string& basename, bool execute, bool clean, bool verbose, 
		      bool grid = false, bool domains=false,
		      bool fill = false, bool elevation = false,
		      bool label = false,
		      size_type isovalue_table_zero_index = std::numeric_limits<size_type>::max()) const;
    int plotmtv2d_P1 (const std::string& basename, bool execute, bool clean, bool verbose,
		      bool grid = false, bool domains=false,
		      bool fill = false, bool elevation = false,
		      bool label = false,
		      size_type isovalue_table_zero_index = std::numeric_limits<size_type>::max()) const;
    int plotmtv2d_P0 (const std::string& basename, bool execute, bool clean, bool verbose,
		      bool grid = false, bool domains=false,
		      bool fill = false, bool elevation = false) const;

    int plotmtv1d_P2 (const std::string& basename, bool execute, bool clean, bool verbose
		     ) const;
    int plotmtv1d_P1 (const std::string& basename, bool execute, bool clean, bool verbose
		     ) const;

    std::ostream& put_vtk_2d_elevation_P1 (std::ostream& vtk) const;

    int plotmtv_deformation_P1 (const std::string& basename, bool execute, bool clean, 
		     bool verbose, const Float& vscale) const;

    int vtk_deformation (const std::string& basename, bool execute, bool clean,
		     bool verbose, const Float& vscale) const;

    int plotmtv_velocity (const std::string& basename, bool execute, bool clean, 
		     bool verbose, const Float& vscale, bool put_domains) const;
    int vtk_velocity (const std::string& basename, bool execute, bool clean, 
		     bool verbose, const Float& vscale) const;

public: 
    // some low-level outputs:
    void extract_line_path (const std::vector<size_t>& mask, std::vector<point>& x, std::vector<Float>& y) const;
    void put_gnuplot1d_data (std::ostream& gdat) const;
protected: 
    int gnuplot1d    (const std::string& basename, bool execute, bool clean, bool verbose
	    	     ) const;
    int gnuplot2d_elevation (const std::string& basename, plot_options& opt) const;
    int vtk2d_Pk     (const std::string& basename, bool execute, bool clean, bool verbose,
		      bool grid = false, bool fill = false) const;
    int vtk3d_Pk     (const std::string& basename, bool execute, bool clean, bool verbose
		     ) const;
    int vtk_elevation   (const std::string& basename, bool execute, bool clean, bool verbose,
		      bool grid = false, bool fill = false) const;

    field cut        (const std::string& basename, bool execute, bool clean, bool verbose,
    		      point normal, point origin) const ;
    geo iso          (const std::string& basename, bool execute, bool clean, bool verbose,
    		      const Float value) const;

    // file format conversion
    std::istream& get_bamg_mbb (std::istream& in);
    std::ostream& put_bamg_mbb (std::ostream& out) const;

    std::istream& get_bamg_metric (std::istream& in);
    std::ostream& put_bamg_metric (std::ostream& out) const;

    std::istream& get_cemagref (std::istream& in);
    std::ostream& put_cemagref (std::ostream& out) const;

    std::ostream& put_gmsh (std::ostream& out) const;

    std::istream& get_mmg3d_metric (std::istream& in);
    std::ostream& put_mmg3d_metric (std::ostream& out) const;
    std::ostream& put_yams_metric (std::ostream& out) const;
public :

// low level accessors:

    void get_dof_value (const basis& b, const geo_element& K,
		tiny_vector<Float>& dofv, size_type i_comp = 0) const;
    void get_dof_value_from_global (const basis& b, const geo_element& K,
		tiny_vector<Float>& dofv, size_type i_comp = 0) const;
     //! same but with K defined in the global mesh
    Float evaluate (const meshpoint& S, size_type i_comp = 0) const;
    void  evaluate (const meshpoint& S, point& result) const;
    Float evaluate_d_dxi (size_type i, const meshpoint& S, size_type i_comp=0) const;
    Float evaluate (const point& x_hat, size_type e, size_type i_comp = 0) const
		{ return evaluate(meshpoint(x_hat,e), i_comp); }

// obsolete:

    void boundary_val( const field&  y,
                       size_type     num_cmp,
                       const domain& d);
    void boundary_val( const field&  y, 
		       size_type     num_cmp,
		       const std::string& dom_name);
    void boundary_val( Float         lambda,
		       size_type     num_cmp,
		       const domain& d);
    void boundary_val( Float         lambda,
		       size_type     num_cmp,
		       const std::string& dom_name) ;
    void boundary_val( Float         (*f)(const point& x),
	 	       const space&  V0, 
		       size_type     num_cmp,
		       const std::string& dom_name);
    void from_boundary_val( const field&  y, 
		       size_type     i_comp,
		       const domain& d);
    field get_comp(const space& X_h, size_type i_comp=0) const; 

    friend class branch;

// data
protected:
    space       _V;
    bool	_locked_boundaries_set;
public :
    vec<Float> u, b;
};
// interpolation:
template <class Function> field interpolate (const space& V, Function f);
template <> field interpolate (const space& V, field uh);

// linear algebra (partial)
field operator - (const field& x);
field operator + (const field& x, const field& y);
field operator + (const Float& lambda, const field& x);
field operator + (const field& x, const Float& lambda);
field operator - (const field& x, const field& y);
field operator - (const Float& lambda, const field& x);
field operator - (const field& x, const Float& lambda);
field operator * (const field& x, const field& y);
field operator * (const Float& lambda, const field& x);
field operator * (const field& x, const Float& lambda);
field operator / (const field& x, const field& y);
field operator / (const Float& lambda, const field& x);
field operator / (const field& x, const Float& lambda);
field sqr (const field& x);
field sqrt (const field& x);
field abs (const field& x);
field pow (const field& x, Float alpha);
Float dot (const field& x, const field& y);
Float norm (const field& x); 
Float max_abs (const field& x); 

//! for a multi-field, returns field = sqrt(sum u(x_i)^2)
field euclidian_norm (const field& u);
field euclidian_norm2 (const field& u);

//! for a tensor field, returns field: sqrt(tau(0,0)^2+2*tau(0,1)^2+..)
field tensor_norm (const field& tau);

// general field & function composition:
// v(x) := f(u(x)) <==> v := compose(f,u)
//                 <==> transform(u,f,v)  : STL-like style
template<class Function>
field compose (Function f, const field& u);
field compose (Float (*f)(const Float&), const field& u);
template<class Function>
void transform (const field& u, Function f, field& v);
void transform (const field& u, Float (*f)(const Float&), field& v);
// supported for backward compatibility purpose, use compose instead:
template<class Function>
field transform (const field& u, Function f);
field transform (const field& u, Float (*f)(const Float&));

// w(x) := f(u(x),v(x)) <==> w := compose(f,u,v)
//                      <==> transform(u,v,f,w)  : STL-like style
template<class Function>
field compose (Function f, const field& u, const field& v);
field compose (Float (*f)(const Float&, const Float&), const field& u, const field& v);
template<class Function>
void transform (const field& u, const field& v, Function f, field& w);
void transform (const field& u, const field& v, Float (*f)(const Float&), field& w);

field trans (const field& tau); // tensor field transposition
field tensor_symmetrize (const space& Sh, const field& tau); 
// when tau is unsymmetric_tensor and Sh is symmetric tensor

// obsolete: concatenation on a product space
field fld_cat2(const field& f1, const field& f2) ;
field fld_cat3(const field& f1, const field& f2, const field& f3) ;

// vector_field & tensor_field:
struct vector_field : std::unary_function<point,point> {
  vector_field(const field& uh);
  point operator() (const point& x) const;
protected:
  field    _uh;
};
struct tensor_field : std::unary_function<point,tensor> {
  tensor_field(const field& mh);
  tensor operator() (const point& x) const;
protected:
  field    _mh;
};
//>field:
// ------------------------------------------------------------------
// vector_field & tensor_field:
// ------------------------------------------------------------------
inline
vector_field::vector_field (const field& uh)
  : std::unary_function<point,point>(), _uh(uh)
{
}
inline
tensor_field::tensor_field (const field& mh)
  : std::unary_function<point,tensor>(), _mh(mh)
{
}
inline
point
vector_field::operator() (const point& x) const
{
    return _uh.vector_evaluate(x);
}
inline
tensor
tensor_field::operator() (const point& x) const 
{
    return _mh.tensor_evaluate(x);
}

}// namespace rheolef

// ------------ inline'd -----------------------------------

#include "rheolef/interpolate.h" // templatized interpolate(V,class_fct)

namespace rheolef { 

inline 
field::field ()
: _V(), _locked_boundaries_set(false), u(), b()
{
}

inline 
field::field (const space& V, const Float& init_value)
: _V(V), _locked_boundaries_set(false), u(), b()
{
   _V.freeze();
   u.resize(V.n_unknown(), init_value);
   b.resize(V.n_blocked(), init_value);
}
inline
field::size_type
field::dimension() const 
{
    return _V.dimension();
}
inline
std::string
field::coordinate_system() const 
{
    return _V.coordinate_system();
}
inline
fem_helper::coordinate_type
field::coordinate_system_type() const 
{
    return _V.coordinate_system_type();
}
inline
field::size_type
field::n_component() const 
{
    return _V.n_component();
}
inline
const space&
field::get_space() const 
{
    return _V;
}

inline
const geo&
field::get_geo() const 
{
    return _V.get_geo();
}

inline
field::size_type
field::size() const
{
    return _V.size();
}

inline
field::size_type 
field::n_unknown() const
{
    return u.size();
}

inline
field::size_type 
field::n_blocked() const
{
    return b.size();
}

inline
Float 
field::min() const
{
    return xmin(u.min(), b.min());
}
inline
Float 
field::max() const
{
    return xmax(u.max(), b.max());
}
inline
Float 
field::max_abs() const
{
    return xmax(u.max_abs(), b.max_abs());
}
inline
std::string
field::get_approx(size_type i_comp) const
{
    return _V.get_approx(i_comp);
}
template<class Function>
inline
void
transform (const field& x, Function f, field& y)
{
    transform(x.u.begin(), x.u.end(), y.u.begin(), f);
    transform(x.b.begin(), x.b.end(), y.b.begin(), f);
}
inline
void
transform (const field& x, Float (*f)(const Float&), field& y)
{
    transform(x.u.begin(), x.u.end(), y.u.begin(), f);
    transform(x.b.begin(), x.b.end(), y.b.begin(), f);
}
template <class Function>
inline
field
compose (Function f, const field& x)
{
    field y (x.get_space());
    transform (x, f, y);
    return y;
}
inline
field
compose (Float (*f)(const Float&), const field& x)
{
    field y (x.get_space());
    transform (x, f, y);
    return y;
}
template <class Function>
inline
field
transform (const field& x, Function f)
{
    warning_macro("vh = transform(uh,f) is obsolete: please, use vh = compose(f,uh) instead.");
    return compose(f,x);
}
inline
field
transform (const field& x, Float (*f)(const Float&))
{
    warning_macro("vh = transform(uh,f) is obsolete: please, use vh = compose(f,uh) instead.");
    return compose(f,x);
}
template<class Function>
inline
void
transform (const field& x, const field& y, Function f, field& z)
{
    transform(x.u.begin(), x.u.end(), y.u.begin(), z.u.begin(), f);
    transform(x.b.begin(), x.b.end(), y.b.begin(), z.b.begin(), f);
}
inline
void
transform (const field& x, const field& y, Float (*f)(const Float&, const Float&), field& z)
{
    transform(x.u.begin(), x.u.end(), y.u.begin(), z.u.begin(), f);
    transform(x.b.begin(), x.b.end(), y.b.begin(), z.b.begin(), f);
}
template <class Function>
inline
field
compose (Function f, const field& x, const field& y)
{
    field z (x.get_space());
    transform (x, y, f, z);
    return z;
}
inline
field
compose (Float (*f)(const Float&, const Float&), const field& x, const field& y)
{
    field z (x.get_space());
    transform (x, y, f, z);
    return z;
}
inline
void
field::boundary_val(const field&   y,
                     size_type     num_cmp,
                     const std::string& dom_name)
{
    const geo& g = get_geo();
    boundary_val (y, num_cmp, g[dom_name]);
}
inline
void
field::boundary_val( Float 	   lambda,
                     size_type     num_cmp,
                     const std::string& dom_name)
{
    const geo& g = get_geo();
    boundary_val (lambda, num_cmp, g[dom_name]);
}
inline 
Float 
norm (const field& x)
{
    return ::sqrt(dot(x,x));
}
inline 
Float 
max_abs (const field& x)
{
    return x.max_abs();
}
inline 
field 
euclidian_norm (const field& x)
{
    return sqrt(euclidian_norm2(x));
}
inline
field
adapt (const class field& c,
        const Float& hcoef)
{
  adapt_option_type opt;
  opt.hcoef = hcoef;
  return adapt (c, opt);
}
inline
fem_helper::valued_field_type
field::get_valued_type() const
{
	return _V.get_valued_type();
}
inline
std::string
field::get_valued() const
{
	return _V.get_valued();
}
inline
point
field::vector_evaluate (const point& x) const
{
    point val;
    evaluate (x, val);
    return val;
}
inline
tensor
field::tensor_evaluate (const point& x) const
{
    tensor val;
    evaluate (x, val);
    return val;
}
// -----------------------------------------------------------
// accessors & modifiers by components or vector/tensor valued
// -----------------------------------------------------------
inline
Float
field::at (size_type i_dof) const
{
    size_type idx = _V.index(i_dof);
    if (! _locked_boundaries_set) {
      if (! _V.is_blocked(i_dof))
        return u (idx);
      else
        return b (idx);
    }
    // in that case, we may have to reconstruct cartesian components 
    // from tangential-normal representation
    if (!_V.is_locked(i_dof)) {
    	if (! _V.is_blocked(i_dof))
	  return u (idx);
	else
	  return b (idx);
    } else { 
        // return blocked * n[i_comp] + unknown * t[i_comp], where n=locked_dir 
      	// and t its direct orthogonal vector
        warning_macro ("Calculating value for i_dof " << i_dof);
	if (_V.is_blocked(i_dof)) // <=> i_comp==1
	    return b(idx)*_V.locked_component(i_dof, 1)
	    	+u(_V.index_locked_with(i_dof))*_V.unlocked_component(i_dof, 1);
	else // then i_comp==0
	    return b(_V.index_locked_with(i_dof))*_V.locked_component(i_dof,0)
	    	+u(idx)*_V.unlocked_component(i_dof,0);
    }
}
inline
Float&
field::at (size_type i_dof)
{
    check_macro (!_locked_boundaries_set,
    	"Direct access to field not permitted: use unset_locked_boundaries() first"); 
    size_type idx = _V.index(i_dof);
    if (! _V.is_blocked(i_dof))
        return u (idx);
    else
        return b (idx);
}
inline
Float
field::at (size_type i_dof, size_type i_comp) const
{
  return at (get_space().start(i_comp) + i_dof);
}
inline
const Float&
field::set_at (size_type i_dof, size_type i_comp, const Float& value)
{
  at (get_space().start(i_comp) + i_dof) = value;
  return value;
}
inline
void
field::incr_at (size_type i_dof, size_type i_comp, const Float& value)
{
  at (get_space().start(i_comp) + i_dof) += value;
}
inline
point
field::vector_at (size_type i_dof) const
{
  point value;
  for (size_type i_comp = 0; i_comp < dimension(); i_comp++)
    value[i_comp] = at (i_dof, i_comp);
  return value;
}
inline
const point&
field::set_vector_at (size_type i_dof, const point& value)
{
  check_macro ( get_valued_type() == fem_helper::vectorial ||
	       (get_valued_type() == fem_helper::scalar && dimension() == 1),
	"field::at: vector field expected, get " << get_valued() << " one.");
  for (size_type i_comp = 0; i_comp < n_component(); i_comp++)
    set_at (i_dof, i_comp, value[i_comp]);
  return value;
}
inline
void
field::incr_vector_at (size_type i_dof, const point& value)
{
  for (size_type i_comp = 0; i_comp < dimension(); i_comp++)
    incr_at (i_dof, i_comp, value[i_comp]);
}
inline
tensor
field::tensor_at (size_type i_dof) const
{
  tensor value;
  fem_helper::coordinate_type   sys_coord = coordinate_system_type();
  fem_helper::valued_field_type valued    = get_valued_type();
  check_macro ( get_valued_type() == fem_helper::tensorial ||
	        get_valued_type() == fem_helper::unsymmetric_tensorial ||
	       (get_valued_type() == fem_helper::scalar && dimension() == 1),
	"field.at: tensor field expected, get " << get_valued() << " one.");
  for (size_type i_comp = 0; i_comp < n_component(); i_comp++) {
    fem_helper::pair_size_type ij_comp = fem_helper::tensor_subscript (valued, sys_coord, i_comp);
    value (ij_comp.first, ij_comp.second) = at (i_dof, i_comp);
    if (ij_comp.first != ij_comp.second) {
      value (ij_comp.second, ij_comp.first) = value (ij_comp.first, ij_comp.second);
    }
  }
  return value;
}
inline
const tensor&
field::set_tensor_at (size_type i_dof, const tensor& value)
{
  fem_helper::coordinate_type   sys_coord = coordinate_system_type();
  fem_helper::valued_field_type valued    = get_valued_type();
  check_macro ( get_valued_type() == fem_helper::tensorial ||
	        get_valued_type() == fem_helper::unsymmetric_tensorial ||
	       (get_valued_type() == fem_helper::scalar && dimension() == 1),
	"field.at: tensor field expected, get " << get_valued() << " one.");
  for (size_type i_comp = 0; i_comp < n_component(); i_comp++) {
    fem_helper::pair_size_type ij_comp = fem_helper::tensor_subscript (valued, sys_coord, i_comp);
    set_at (i_dof, i_comp, value(ij_comp.first,ij_comp.second));
  }
  return value;
}
inline
void
field::incr_tensor_at (size_type i_dof, const tensor& value)
{
  fem_helper::coordinate_type   sys_coord = coordinate_system_type();
  fem_helper::valued_field_type valued    = get_valued_type();
  check_macro ( get_valued_type() == fem_helper::tensorial ||
	        get_valued_type() == fem_helper::unsymmetric_tensorial ||
	       (get_valued_type() == fem_helper::scalar && dimension() == 1),
	"field.at: tensor field expected, get " << get_valued() << " one.");
  for (size_type i_comp = 0; i_comp < n_component(); i_comp++) {
    fem_helper::pair_size_type ij_comp = fem_helper::tensor_subscript (valued, sys_coord, i_comp);
    incr_at (i_dof, i_comp, value(ij_comp.first,ij_comp.second));
  }
}
// ---------------------------------------------------------
// u["boundary"] = 0; 		fields on domains...
// ---------------------------------------------------------
struct field_on_domain {
    typedef field::size_type size_type;
    field_on_domain();
    field_on_domain(field& x, const domain& d, size_type i_comp = std::numeric_limits<size_type>::max());
    const Float& operator = (const Float& lambda);
    const field& operator = (const field& y);
    field*    _px;
    domain    _d;
    size_type _i_comp;
};
inline
field_on_domain::field_on_domain()
 : _px(0), _d(), _i_comp(std::numeric_limits<size_type>::max())
{
}
inline
field_on_domain::field_on_domain(field& x, const domain& d, size_type i_comp) 
 : _px(&x), _d(d), _i_comp(i_comp)
{
}
inline
field_on_domain
field::operator [] (const domain& d)
{
    return field_on_domain(*this,d);
}
inline
field_on_domain
field::operator [] (const std::string& dom_name)
{
    const geo& g = get_geo();
    return operator [] (g[dom_name]);
}
// ---------------------------------------------------------
// u[1]["boundary"] = 0; 	multicomponent fields...
// ---------------------------------------------------------
struct field_component {
    typedef field::size_type size_type;
    field_component();
    field_component(field& x, size_type i);
    const Float& operator = (const Float& lambda);
    const field& operator = (const field& y);
    field operator = (const field_component& y);
    field_component& operator += (const field& y);
    field_on_domain operator [] (const std::string& dom_name);
    field_on_domain operator [] (const domain& dom);
    friend std::ostream& operator << (std::ostream& s, const field_component& x);
    friend field operator + (const field_component&, const field_component&);
    friend field operator - (const field_component&, const field_component&);
    friend field operator * (const field_component&, const field_component&);
    friend field operator / (const field_component&, const field_component&);
    friend field compose (Float (*f)(const Float&), const field_component&);
    friend field transform (const field_component&, Float (*f)(const Float&));
    friend void transform (const field_component&, Float (*f)(const Float&), field&);
#ifdef TODO
    template<class Function>
    friend field compose (Function f, const field_component&);
    template<class Function>
    friend field transform (const field_component&, Function f);
    template<class Function>
    friend void transform (const field_component&, Function f, field&);
#endif // TODO
    friend field sqr (const field_component&);
#ifdef TODO
    friend field sqrt(const field_component&);
    friend field abs(const field_component&);
#endif // TODO
    size_type u_size() const;
    size_type b_size() const;
    vec<Float>::const_iterator u_begin() const;
    vec<Float>::const_iterator b_begin() const;
    vec<Float>::const_iterator u_end() const;
    vec<Float>::const_iterator b_end() const;
    vec<Float>::iterator u_begin();
    vec<Float>::iterator b_begin();
    vec<Float>::iterator u_end();
    vec<Float>::iterator b_end();
    field*          _px;
    size_type       _i_comp;
};
inline
field_component::field_component()
 : _px(0), _i_comp(0)
{
}
inline
field_component::field_component(field& x, size_type i)
 : _px(&x), _i_comp(i)
{
}
inline
field_component 
field::operator [] (size_type i_comp)
{
    return field_component (*this, i_comp);
}
inline
field_on_domain
field_component::operator [] (const domain& d)
{
    return field_on_domain(*_px,d,_i_comp);
}
inline
field_on_domain
field_component::operator [] (const std::string& dom_name)
{
    const geo& g = (*_px).get_geo();
    return operator [] (g[dom_name]);
}
// ---------------------------------------------------------
// field w = u[1]*v[0]
// ---------------------------------------------------------
struct const_field_component {
    typedef field::size_type size_type;
    const_field_component();
    const_field_component(const field& x, size_type i);
    const_field_component(const field_component&);
    friend std::ostream& operator << (std::ostream& s, const const_field_component& x);
    friend field operator + (const const_field_component&, const const_field_component&);
    friend field operator - (const const_field_component&, const const_field_component&);
    friend field operator * (const const_field_component&, const const_field_component&);
    friend field operator / (const const_field_component&, const const_field_component&);
    friend field compose (Float (*f)(const Float&), const const_field_component&);
    friend field transform (const const_field_component&, Float (*f)(const Float&));
    friend void transform (const const_field_component&, Float (*f)(const Float&), field&);
#ifdef TODO
    template<class Function>
    friend field compose (Function f, const const_field_component&);
    template<class Function>
    friend void transform (const const_field_component&, Function f, field&);
    template<class Function>
    friend field transform (const const_field_component&, Function f);
#endif // TODO
    friend field sqr (const const_field_component&);
#ifdef TODO
    friend field sqrt(const const_field_component&);
    friend field abs(const const_field_component&);
#endif // TODO
    size_type u_size() const;
    size_type b_size() const;
    vec<Float>::const_iterator u_begin() const;
    vec<Float>::const_iterator b_begin() const;
    vec<Float>::const_iterator u_end() const;
    vec<Float>::const_iterator b_end() const;
    const field*    _px;
    size_type       _i_comp;
};
inline
const_field_component::const_field_component()
 : _px(0), _i_comp(0)
{
}
inline
const_field_component::const_field_component(const field& x, size_type i)
 : _px(&x), _i_comp(i)
{
}
inline
const_field_component::const_field_component(const field_component& x)
 : _px(x._px), _i_comp(x._i_comp)
{
}
inline
const_field_component 
field::operator [] (size_type i_comp) const
{
    return const_field_component (*this, i_comp);
}
inline
field
field::operator += (const field& f)
{
    u += f.u;
    b += f.b;
    return *this;
}
inline
field
field::operator -= (const field& f)
{
    u -= f.u;
    b -= f.b;
    return *this;
}
}// namespace rheolef
# endif /* _RHEO_FIELD_H */