.\" .de Id .. .de Sp .if n .sp .if t .sp 0.4 .. .TH field 2rheolef "rheolef-6.7" "rheolef-6.7" "rheolef-6.7" .\" label: /*Class:field .SH NAME \fBfield\fP - piecewise polynomial finite element field .SH 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). .PP 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: .\" begin_example .Sp .nf geo omega ("circle"); space Xh (omega, "P1"); Xh.block ("boundary"); field uh (Xh); uh ["boundary"] = 0; .Sp .fi .\" end_example .PP .\" skip: @findex interpolate .SH INTERPOLATION Interpolation of a function \fBu\fP in a field \fBuh\fP with respect to the interpolation writes: .\" begin_example .Sp .nf Float u (const point& x) { return x[0]*x[1]; } ... field uh = interpolate (Xh, u); .Sp .fi .\" end_example .PP .\" skip: @findex dual .SH LINEAR ALGEBRA Linear algebra, such as \fBuh+vh\fP, \fBuh-vh\fP and \fBlambda*uh + mu*vh\fP, where \fBlambda\fP and \fBmu\fP are of type \fBFloat\fP, are supported. The duality product between two fields \fBlh\fP and \fBvh\fP writes simply \fBdual(lh,vh)\fP: for discrete fields, it corresponds to a simple Euclidian dot product in \fBIR^n\fP. The application of a bilinear form (see form(2)) writes \fBm(uh,vh)\fP and is equivalent to \fBdual(m*uh,vh)\fP. .PP .SH NON-LINEAR ALGEBRA Non-linear operations, such as \fBsqrt(uh)\fP or \fB1/uh\fP are also available. Notice that non-linear operations do not returns in general picewise polynomials: the value returned by \fBsqrt(uh)\fP may be filtered by \fBinterpolate\fP, .\" begin_example .Sp .nf field vh = interpolate (Xh, sqrt(uh)); .Sp .fi .\" end_example the Lagrange interpolant, to becomes a piecewise polynomial: All standard unary and binary math functions \fBabs, cos, sin\fP... are available on fields. Also \fBsqr(uh)\fP, the square of a field, and \fBmin(uh,vh)\fP, \fBmax(uh,vh)\fP are provided. Binary functions can be used also with a scalar, as in .\" begin_example .Sp .nf field vh = interpolate (Xh, max (abs(uh), 0)); field wh = interpolate (Xh, pow (abs(uh), 1./3)); .Sp .fi .\" end_example For applying a user-provided function to a field, use the \fBcompose\fP function: .\" begin_example .Sp .nf field vh = interpolate(Xh, compose(f, uh)); field wh = interpolate(Xh, compose(f, uh, vh)); .Sp .fi .\" end_example The composition supports also general unary and binary class-functions. Also, the multiplication \fBuh*vh\fP and the division \fBuh/vh\fP returns a result that is not in the same discrete finite element space: its result may be filtered by the \fBinterpolate\fP operator: .\" begin_example .Sp .nf field wh = interpolate(Xh, uh*vh); .Sp .fi .\" end_example Any function or class function can be used in nonlinear expressions: the function is interpolated in the specified finite element space. .PP .\" skip: @findex normal There is a special predefined class-function named \fBnormal\fP that represents the outer unnit normal vector on a boundary domain or surfacic mesh: .\" begin_example .Sp .nf 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()); .Sp .fi .\" end_example 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. .PP .SH ACCESS BY DOMAIN The restriction of a field to a geometric domain, says \fB"boundary"\fP writes \fBuh["boundary"]\fP: it represents the trace of the field on the boundary: .\" begin_example .Sp .nf space Xh (omega, "P1"); uh["boundary"] = 0; .Sp .fi .\" end_example Extraction of the trace as a field is also possible: .\" begin_example .Sp .nf field wh = uh["boundary"]; .Sp .fi .\" end_example The space associated to the trace writes \fBwh.get_space()\fP and is equivalent to \fBspace(omega["boundary"], "P1")\fP. See see space(2). .PP .SH VECTOR VALUED FIELD A vector-valued field contains several components, as: .\" begin_example .Sp .nf space Xh (omega, "P2", "vector"); field uh (Xh); field vh = uh[0] - uh[1]; field nh = norm (uh); .Sp .fi .\" end_example The \fBnorm\fP function returns the euclidian norm of the vector-valuated field at each degree of freedom: its result is a scalar field. .PP .SH TENSOR VALUED FIELD A tensor-valued field can be constructed and used as: .\" begin_example .Sp .nf space Th (omega, "P1d", "tensor"); field sigma_h (Xh); field trace_h = sigma(0,0) + sigma_h(1,1); field nh = norm (sigma_h); .Sp .fi .\" end_example The \fBnorm\fP 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. .PP .SH GENERAL MULTI-COMPONENT INTERFACE A general multi-component field writes: .\" begin_example .Sp .nf 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 .Sp .fi .\" end_example 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. .PP For any field \fBxh\fP, the string \fBxh.valued()\fP returns \fB"scalar"\fP for a scalar field and \fB"vector"\fP for a vector-valued one. Other possible valued are \fB"tensor"\fP and \fB"other"\fP. The \fBxh.size()\fP returns the number of field components. When the field is scalar, it returns zero by convention, and \fBxh[0]\fP is undefined. A vector-valued field has \fBd\fP components, where \fBd=omega.dimension()\fP. A tensor-valued field has \fBd*(d+1)/2\fP components, where \fBd=omega.dimension()\fP. .PP .SH 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 \fBuh.u()\fP and \fBuh.b()\fP, and a read-and-write one, as \fBuh.set_u()\fP and \fBuh.set_b()\fP, see see vec(2). .PP .SH 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 \fBuh\fP. The number of dofs is \fBuh.ndof()\fP and any dof can be accessed via \fBuh.dof(idof)\fP. A non-local dof at the partition interface can be obtain via \fBuh.dis_dof(dis_idof)\fP where \fBdis_idof\fP is the (global) distribued index assoiated to the distribution \fBuh.ownership()\fP. .PP For performances, a STL-like iterator interface is available, with \fBuh.begin_dof()\fP and \fBuh.end_dof()\fP returns iterators to the dofs on the current processor. See see disarray(2) for more about distributed arrays. .PP For convenience, \fBuh.max()\fP, \fBuh.min()\fP and \fBuh.max_abs()\fP retuns respectively the maximum, minimum and maximum of the absolute value of the degrees of freedom. .PP .SH FILE FORMAT TODO .PP .SH IMPLEMENTATION NOTE The field expression use the expression template technics in order to avoid temporaries when evaluating complex expressions. .\" skip start:AUTHOR: .\" skip start:DATE: .\" END .SH IMPLEMENTATION .\" begin_example .Sp .nf template class field_basic : public std::unary_function::type>,T> { public : // typedefs: typedef typename std::size_t size_type; typedef M memory_type; typedef T scalar_type; typedef typename float_traits::type float_type; // typedef undeterminated_basic value_type; // TODO typedef T value_type; // TO_CLEAN typedef space_constant::valued_type valued_type; typedef geo_basic geo_type; typedef space_basic space_type; typedef typename vec::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::max()); void resize ( const space_type& V, const T& init_value = std::numeric_limits::max()); field_basic& operator= (const T&); field_basic& operator= (const field_basic&); #ifdef TO_CLEAN field_basic (const field_indirect&); field_basic& operator= (const field_indirect&); field_basic (const field_indirect_const&); field_basic& operator= (const field_indirect_const&); field_basic (const field_component&); field_basic& operator= (const field_component&); field_basic (const field_component_const&); field_basic& operator= (const field_component_const&); #endif // TO_CLEAN // linear expressions: template ::value && ! details::is_field::value >::type> field_basic (const Expr& expr); template ::value && ! details::is_field::value >::type> field_basic& operator= (const Expr&); // initializer list (c++ 2011): #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST field_basic (const std::initializer_list >& init_list); field_basic& operator= (const std::initializer_list >& 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& u() const { return _u; } const vec& b() const { return _b; } vec& set_u() { dis_dof_indexes_requires_update(); return _u; } vec& 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 operator[] (const geo_basic& dom); field_indirect_const operator[] (const geo_basic& dom) const; field_indirect operator[] (std::string dom_name); field_indirect_const operator[] (std::string dom_name) const; // accessors by components: size_type size() const { return _V.size(); } field_component operator[] (size_type i_comp); field_component_const operator[] (size_type i_comp) const; field_component operator() (size_type i_comp, size_type j_comp); field_component_const 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& x, size_type i_comp = 0) const; T operator() (const point_basic& x) const { return dis_evaluate (x,0); } point_basic dis_vector_evaluate (const point_basic& 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& hat_xq, size_type i_comp = 0) const; // propagate changed values shared at partition boundaries to others procs void dis_dof_update() const; template void assembly_internal ( const geo_basic& dom, const geo_basic& band, const band_basic& gh, const Expr& expr, const quadrature_option_type& qopt, bool is_on_band); template void assembly ( const geo_basic& domain, const Expr& expr, const quadrature_option_type& qopt); template void assembly ( const band_basic& 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 _u; mutable vec _b; mutable bool _dis_dof_indexes_requires_update; mutable bool _dis_dof_assembly_requires_update; }; template idiststream& operator >> (odiststream& ips, field_basic& u); template odiststream& operator << (odiststream& ops, const field_basic& uh); typedef field_basic field; typedef field_basic field_sequential; .Sp .fi .\" end_example .\" LENGTH = 6 .SH SEE ALSO space(2), form(2), space(2), space(2), vec(2), disarray(2)