GNU.WIKI: The GNU/Linux Knowledge Base

#### 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  arrays  of  degrees-of-freedom  (dof)
associated respectively to  blocked  and  unknown  dofs.  Blocked  dofs
corresponds  to  Dirichlet  boundary  conditions, as specified by space
is allowed, as uh.b and uh.u: 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 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.
```

#### 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;
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                   (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):

#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 { 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:
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,
bool                          is_on_band);
template <class Expr>
void assembly (
const geo_basic<T,M>&         domain,
const Expr&                   expr,
template <class Expr>
void assembly (
const band_basic<T,M>&        gh,
const Expr&                   expr,

protected:
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)`