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 (See see
space(2)). For
simpliity, direct public access to these array 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.
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,
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_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;
SEE ALSO¶
space(2),
form(2),
space(2),
space(2),
vec(2),
array(2)