bionic (2) field.2rheolef.gz

Provided by: librheolef-dev_6.7-6_amd64 bug

NAME

       field - piecewise polynomial finite element field

DESCRIPTION

       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

       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

       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 ALGEBRA

       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.

ACCESS BY DOMAIN

       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).

VECTOR VALUED FIELD

       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.

TENSOR VALUED 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.

GENERAL MULTI-COMPONENT INTERFACE

       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 d=omega.dimension().

BLOCKED AND UNBLOCKED ARRAYS

       The  field  class contains two vectors 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)).  Access to these vectors is allowed via some accessors: a read-only one, as uh.u() and uh.b(),
       and a read-and-write one, as uh.set_u() and uh.set_b(), see see vec(2).

LOW-LEVEL DEGREE-OF-FREEDOM ACCESS

       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 dofs on the current processor.  See see disarray(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.

FILE FORMAT

       TODO

IMPLEMENTATION NOTE

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

IMPLEMENTATION

       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;
           typedef typename vec<T,M>::dis_reference dis_reference;
           class iterator;
           class const_iterator;

       // allocator/deallocator:

           field_basic();

           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<T,M>& operator=  (const T&);
           field_basic<T,M>& operator=  (const field_basic<T,M>&);

       #ifdef TO_CLEAN
           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>&);
       #endif // TO_CLEAN

           // linear expressions:
           template <class Expr, class Sfinae
                           = typename std::enable_if<
                                    details::is_field_expr_v2_linear_arg<Expr>::value
                               && ! details::is_field<Expr>::value
                             >::type>
           field_basic (const Expr& expr);

           template <class Expr, class Sfinae
                           = typename std::enable_if<
                                    details::is_field_expr_v2_linear_arg<Expr>::value
                               && ! details::is_field<Expr>::value
                             >::type>
           field_basic<T, M>& operator=  (const Expr&);

       // initializer list (c++ 2011):

       #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
           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);
       #endif // _RHEOLEF_HAVE_STD_INITIALIZER_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 { return _u; }
           const vec<T,M>&     b() const { return _b; }
                 vec<T,M>& set_u()       { dis_dof_indexes_requires_update(); return _u; }
                 vec<T,M>& set_b()       { dis_dof_indexes_requires_update(); 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;
           // write access to non-local dis_idof changes to others procs
           dis_reference dis_dof_entry (size_type dis_idof);

           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:
       public:

           // 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);

       protected:
           void dis_dof_indexes_requires_update() const;
           void dis_dof_assembly_requires_update() const;

       // data:
                   space_type   _V;
           mutable vec<T,M>     _u;
           mutable vec<T,M>     _b;
           mutable bool         _dis_dof_indexes_requires_update;
           mutable bool         _dis_dof_assembly_requires_update;
       };
       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;

SEE ALSO

       space(2), form(2), space(2), space(2), vec(2), disarray(2)