Provided by: librheolef-dev_6.6-1build2_amd64 bug


       field - piecewise polynomial finite element field


       Store  degrees  of  freedom associated to a mesh and a piecewise polynomial approximation,
       with respect to the numbering defined by the underlying space(2).

       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:

               geo omega ("circle");
               space Xh (omega, "P1");
               Xh.block ("boundary");
               field uh (Xh);
               uh ["boundary"] = 0;


       Interpolation of a function u in a field uh with respect to the interpolation writes:

               Float u (const point& x) { return x[0]*x[1]; }
               field uh = interpolate (Xh, u);


       Linear algebra, such as uh+vh, uh-vh and lambda*uh + mu*vh, where lambda  and  mu  are  of
       type Float, are supported.  The duality product between two fields lh and vh writes simply
       dual(lh,vh): for discrete fields, it corresponds to a  simple  Euclidian  dot  product  in
       IR^n.   The application of a bilinear form (see form(2)) writes m(uh,vh) and is equivalent
       to dual(m*uh,vh).


       Non-linear operations, such as sqrt(uh) or 1/uh are  also  available.   Notice  that  non-
       linear  operations  do  not returns in general picewise polynomials: the value returned by
       sqrt(uh) may be filtered by interpolate,

               field vh = interpolate (Xh, sqrt(uh));

       the Lagrange interpolant, to becomes a piecewise polynomial: All standard unary and binary
       math  functions  abs,  cos,  sin... are available on fields. Also sqr(uh), the square of a
       field, and min(uh,vh), max(uh,vh) are provided. Binary functions can be used also  with  a
       scalar, as in

               field vh = interpolate (Xh, max (abs(uh), 0));
               field wh = interpolate (Xh, pow (abs(uh), 1./3));

       For applying a user-provided function to a field, use the compose function:

               field vh = interpolate(Xh, compose(f, uh));
               field wh = interpolate(Xh, compose(f, uh, vh));

       The  composition  supports  also  general  unary  and  binary  class-functions.  Also, the
       multiplication uh*vh and the division uh/vh returns a result  that  is  not  in  the  same
       discrete finite element space: its result may be filtered by the interpolate operator:

               field wh = interpolate(Xh, uh*vh);

       Any  function  or  class  function  can  be used in nonlinear expressions: the function is
       interpolated in the specified finite element space.

       There is a special predefined class-function named normal that represents the outer  unnit
       normal vector on a boundary domain or surfacic mesh:

               size_t k = omega.order();
               string n_approx = "P" + itos(k-1) + "d";
               space Nh (omega["boundary"], n_approx, "vector");
               field nh = interpolate(Nh, normal());

       The  normal()  function could appear in any nonlinear field expression: it is evaluated on
       the fly, based on the current mesh.  Notice that when using  isoparametric  elements,  the
       normal  vector  is  no  more constant along any face of the mesh.  Also, on general curved
       domains, the unit normal vector is discontinuous accross boundary element interfaces.


       The restriction of a field to a geometric domain, says "boundary"  writes  uh["boundary"]:
       it represents the trace of the field on the boundary:

               space Xh (omega, "P1");
               uh["boundary"] = 0;

       Extraction of the trace as a field is also possible:

               field wh = uh["boundary"];

       The   space   associated   to  the  trace  writes  wh.get_space()  and  is  equivalent  to
       space(omega["boundary"], "P1").  See see space(2).


       A vector-valued field contains several components, as:

               space Xh (omega, "P2", "vector");
               field uh (Xh);
               field vh = uh[0] - uh[1];
               field nh = norm (uh);

       The norm function returns the euclidian norm of the vector-valuated field at  each  degree
       of freedom: its result is a scalar field.


       A tensor-valued field can be constructed and used as:

               space Th (omega, "P1d", "tensor");
               field sigma_h (Xh);
               field trace_h = sigma(0,0) + sigma_h(1,1);
               field nh = norm (sigma_h);

       The  norm  function returns the euclidian norm of the tensor-valuated field at each degree
       of freedom: its result is a scalar  field.   Notice  that,  as  tensor-valued  fields  are
       symmetric, extra-diagonals are counted twice.


       A general multi-component field writes:

               space Th (omega, "P1d", "tensor");
               space Vh (omega, "P2", "vector");
               space Qh (omega, "P1");
               space Xh = Th*Vh*Qh;
               field xh (Xh);
               field tau_h = xh[0]; // tensor-valued
               field uh    = xh[1]; // vector-valued
               field qh    = xh[2]; // scalar

       Remark  the  hierarchical  multi-component field structure: the first-component is tensor-
       valued and the second-one is vector-valued.  There is no limitation upon the  hierarchical
       number of levels in use.

       For  any field xh, the string xh.valued() returns "scalar" for a scalar field and "vector"
       for a vector-valued one. Other possible valued are "tensor" and  "other".   The  xh.size()
       returns  the  number  of  field  components.  When the field is scalar, it returns zero by
       convention, and xh[0] is  undefined.   A  vector-valued  field  has  d  components,  where
       d=omega.dimension().     A   tensor-valued   field   has   d*(d+1)/2   components,   where


       The field class contains two arrays of degrees-of-freedom (dof) associated respectively to
       blocked  and  unknown  dofs. Blocked dofs corresponds to Dirichlet boundary conditions, as
       specified by space (See see space(2)).  For simpliity, direct public access to these array
       is allowed, as uh.b and uh.u: see see vec(2).


       The  field  class  provides  a  STL-like container interface for accessing the degrees-of-
       freedom (dof) of a finite element field uh. The number of dofs is uh.ndof()  and  any  dof
       can  be  accessed  via  uh.dof(idof).   A  non-local dof at the partition interface can be
       obtain via uh.dis_dof(dis_idof) where dis_idof is the (global) distribued index  assoiated
       to the distribution uh.ownership().

       For  performances,  a  STL-like  iterator  interface is available, with uh.begin_dof() and
       uh.end_dof() returns iterators to the arrays of dofs on the current  processor.   See  see
       array(2) for more about distributed arrays.

       For  convenience,  uh.max(),  uh.min()  and  uh.max_abs() retuns respectively the maximum,
       minimum and maximum of the absolute value of the degrees of freedom.




       The field expression use the expression template technics in order  to  avoid  temporaries
       when evaluating complex expressions.


       template <class T, class M = rheo_default_memory_model>
       class field_basic : public std::unary_function<point_basic<typename scalar_traits<T>::type>,T> {
       public :
       // typedefs:

           typedef typename std::size_t            size_type;
           typedef M                               memory_type;
           typedef T                               scalar_type;
           typedef typename float_traits<T>::type  float_type;
       //  typedef undeterminated_basic<T>         value_type; // TODO
           typedef T                               value_type; // TO_CLEAN
           typedef space_constant::valued_type     valued_type;
           typedef geo_basic  <float_type,M>       geo_type;
           typedef space_basic<float_type,M>       space_type;
           class iterator;
           class const_iterator;

       // allocator/deallocator:


           explicit field_basic (
               const space_type& V,
               const T& init_value = std::numeric_limits<T>::max());

           void resize (
               const space_type& V,
               const T& init_value = std::numeric_limits<T>::max());

           field_basic                   (const field_indirect<T,M>&);
           field_basic<T, M>& operator=  (const field_indirect<T,M>&);
           field_basic                   (const field_indirect_const<T,M>&);
           field_basic<T, M>& operator=  (const field_indirect_const<T,M>&);
           field_basic                   (const field_component<T,M>&);
           field_basic<T, M>& operator=  (const field_component<T,M>&);
           field_basic                   (const field_component_const<T,M>&);
           field_basic<T, M>& operator=  (const field_component_const<T,M>&);
           template <class Expr> field_basic                   (const field_expr<Expr>&);
           template <class Expr> field_basic<T, M>& operator=  (const field_expr<Expr>&);
           field_basic<T, M>& operator=  (const T&);

       // initializer list (c++ 2011):

           field_basic (const std::initializer_list<field_concat_value<T,M> >& init_list);
           field_basic<T,M>& operator= (const std::initializer_list<field_concat_value<T,M> >& init_list);

       // accessors:

           const space_type&  get_space()  const { return _V; }
           const geo_type&    get_geo()    const { return _V.get_geo(); }
           std::string        stamp()      const { return _V.stamp(); }
           std::string        get_approx() const { return _V.get_approx(); }
           valued_type        valued_tag() const { return _V.valued_tag(); }
           const std::string& valued()     const { return _V.valued(); }

       // accessors & modifiers to unknown & blocked parts:

           const vec<T,M>&     u() const { dis_dof_update_needed(); return _u; }
           const vec<T,M>&     b() const { dis_dof_update_needed(); return _b; }
                 vec<T,M>& set_u()       { return _u; }
                 vec<T,M>& set_b()       { return _b; }

       // accessors to extremas:

           T min() const;
           T max() const;
           T max_abs() const;
           T min_abs() const;

       // accessors by domains:

           field_indirect<T,M>        operator[] (const geo_basic<T,M>& dom);
           field_indirect_const<T,M>  operator[] (const geo_basic<T,M>& dom) const;
           field_indirect<T,M>        operator[] (std::string dom_name);
           field_indirect_const<T,M>  operator[] (std::string dom_name) const;

       // accessors by components:

           size_type size() const { return _V.size(); }
           field_component<T,M>       operator[] (size_type i_comp);
           field_component_const<T,M> operator[] (size_type i_comp) const;
           field_component<T,M>       operator() (size_type i_comp, size_type j_comp);
           field_component_const<T,M> operator() (size_type i_comp, size_type j_comp) const;

       // accessors by degrees-of-freedom (dof):

           const distributor& ownership() const { return get_space().ownership(); }
           const communicator& comm() const { return ownership().comm(); }
           size_type     ndof() const { return ownership().size(); }
           size_type dis_ndof() const { return ownership().dis_size(); }
                 T& dof (size_type idof);
           const T& dof (size_type idof) const;
           const T& dis_dof (size_type dis_idof) const;
           iterator begin_dof();
           iterator end_dof();
           const_iterator begin_dof() const;
           const_iterator end_dof() const;

       // input/output:

           idiststream& get (idiststream& ips);
           odiststream& put (odiststream& ops) const;
           odiststream& put_field (odiststream& ops) const;

       // evaluate uh(x) where x is given locally as hat_x in K:

           T dis_evaluate (const point_basic<T>& x, size_type i_comp = 0) const;
           T operator()   (const point_basic<T>& x) const { return dis_evaluate (x,0); }
           point_basic<T> dis_vector_evaluate (const point_basic<T>& x) const;

       // internals:

           // evaluate uh(x) where x is given locally as hat_x in K:
           // requiers to call field::dis_dof_upgrade() before.
           T evaluate (const geo_element& K, const point_basic<T>& hat_xq, size_type i_comp = 0) const;

           // propagate changed values shared at partition boundaries to others procs
           void dis_dof_update() const;

           template <class Expr>
           void assembly_internal (
               const geo_basic<T,M>&         dom,
               const geo_basic<T,M>&         band,
               const band_basic<T,M>&        gh,
               const Expr&                   expr,
               const quadrature_option_type& qopt,
               bool                          is_on_band);
           template <class Expr>
           void assembly (
               const geo_basic<T,M>&         domain,
               const Expr&                   expr,
               const quadrature_option_type& qopt);
           template <class Expr>
           void assembly (
               const band_basic<T,M>&        gh,
               const Expr&                   expr,
               const quadrature_option_type& qopt);

           void dis_dof_update_internal() const;
           void dis_dof_update_needed() const;

       // data:
           space_type   _V;
           vec<T,M>     _u;
           vec<T,M>     _b;
           mutable bool _dis_dof_update_needed;
       template <class T, class M>
       idiststream& operator >> (odiststream& ips, field_basic<T,M>& u);

       template <class T, class M>
       odiststream& operator << (odiststream& ops, const field_basic<T,M>& uh);

       typedef field_basic<Float> field;
       typedef field_basic<Float,sequential> field_sequential;


       space(2), form(2), space(2), space(2), vec(2), array(2)