librheolef-dev 6.6-1build2 / usr / include / rheolef / form.h

# ifndef _RHEOLEF_FORM_H
# define _RHEOLEF_FORM_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

#include "rheolef/csr.h"
#include "rheolef/field.h"
#include "rheolef/quadrature.h"
#include "rheolef/form_option_type.h"
namespace rheolef { 

// these classes are used for allocator from the std::initializer_list
template <class T, class M> class form_concat_value;
template <class T, class M> class form_concat_line;

// forward declaration:
template <class T, class M> class band_basic;

/*Class:form
NAME:  @code{form} - representation of a finite element bilinear form
DESCRIPTION:       
  The form class groups four sparse matrix, associated to 
  a bilinear form on two finite element spaces:
  @example
       a: U*V   ----> IR
         (u,v)  |---> a(u,v)
  @end example
  The operator @code{A} associated to the bilinear form is
  defined by:
  @example
       A: U  ----> V'
          u  |---> A(u)
  @end example
  where @code{u} and @code{v} are fields (@pxref{field class}),
  and @code{A(u)} is such that @code{a(u,v)=<A(u),v>}
  for all u in U and v in V and where @code{<.,.>} denotes
  the duality product between V and V'.
  Since V is a finite dimensional spaces, the duality product is
  the euclidian product in IR^dim(V).

  Since both U and V are finite dimensional spaces,
  the linear operator can be represented by a matrix.
  The @code{form} class is represented by four sparse matrix
  in @code{csr} format (@pxref{csr class}),
  associated to unknown and blocked degrees of freedom 
  of origin and destination spaces (@pxref{space class}).
EXAMPLE:
  
  The operator A associated to a bilinear form a(.,.) by
  the relation (Au,v) = a(u,v) could be applied by using
  a sample matrix notation A*u, as shown by the following code:
  @example
      geo omega("square");
      space V (omega,"P1");
      form a (V,V,"grad_grad");
      field uh = interpolate (fct, V);
      field vh = a*uh;
      cout << v;
  @end example

  The form-field @code{vh=a*uh} operation is equivalent to
  the following matrix-vector operations:
  @example
     vh.set_u() = a.uu()*uh.u() + a.ub()*uh.b();
     vh.set_b() = a.bu()*uh.u() + a.bb()*uh.b();
  @end example

ALGEBRA:
  Forms, as matrices (see @ref{csr class}), support linear algebra:
  Adding or substracting two forms writes @code{a+b} and @code{a-b}, respectively,
  and multiplying a form by a field @code{uh} writes @code{a*uh}.
  Thus, any linear combination of forms is available.

@cindex quarature formula
WEIGHTED FORM:
  A weighted form is a form with an extra weight function @code{w(x)}, e.g.:
  @example
                   /
                  |
       a(uh,vh) = |       grad(uh).grad(vh) w(x) dx
                  |
                 / Omega
  @end example
  In the present implementation, @code{w} can be any field, function or class-function
  or any nonlinear field expression (@pxref{field class}).
  As the integration cannot be performed exactly in general, a quadrature formula
  can be supplied.
  This feature is extensively used when solving nonlinear problems.
 
SEE ALSO: "space"(3), "field"(3), "csr"(3)
AUTHOR: Pierre.Saramito@imag.fr
DATE:   2 july 1997
METHODS: @form
End:
*/

//<form:
template<class T, class M>
class form_basic {
public :
// typedefs:

    typedef typename csr<T,M>::size_type    size_type;
    typedef T                               value_type;
    typedef typename scalar_traits<T>::type float_type;
    typedef geo_basic<float_type,M>         geo_type;
    typedef space_basic<float_type,M>       space_type;

// allocator/deallocator:

    form_basic ();
    form_basic (const form_basic<T,M>&);


// allocators from initializer list (c++ 2011):

#ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
    form_basic (const std::initializer_list<form_concat_value<T,M> >& init_list);
    form_basic (const std::initializer_list<form_concat_line <T,M> >& init_list);
#endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST

// accessors:

    const space_type& get_first_space() const;
    const space_type& get_second_space() const;
    const geo_type&   get_geo() const;

    const communicator& comm() const;

// linear algebra:

    form_basic<T,M>  operator+  (const form_basic<T,M>& b) const;
    form_basic<T,M>  operator-  (const form_basic<T,M>& b) const;
    form_basic<T,M>  operator*  (const form_basic<T,M>& b) const;
    form_basic<T,M>& operator*= (const T& lambda);
    field_basic<T,M> operator*  (const field_basic<T,M>& xh) const;
    field_basic<T,M> trans_mult (const field_basic<T,M>& yh) const;
    float_type operator () (const field_basic<T,M>& uh, const field_basic<T,M>& vh) const;

// io:

    odiststream& put (odiststream& ops, bool show_partition = true) const;
    void dump (std::string name) const;

// accessors & modifiers to unknown & blocked parts:

    const csr<T,M>&     uu() const { return _uu; }
    const csr<T,M>&     ub() const { return _ub; }
    const csr<T,M>&     bu() const { return _bu; }
    const csr<T,M>&     bb() const { return _bb; }
          csr<T,M>& set_uu()       { return _uu; }
          csr<T,M>& set_ub()       { return _ub; }
          csr<T,M>& set_bu()       { return _bu; }
          csr<T,M>& set_bb()       { return _bb; }

// data
protected:
    space_type  _X;
    space_type  _Y;
    csr<T,M>    _uu;
    csr<T,M>    _ub;
    csr<T,M>    _bu;
    csr<T,M>    _bb;

// internals:
public: 
    // with vf expression arg
    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 form_option_type&       fopt,
        bool                          is_on_band);
    template <class Expr>
    void assembly (
        const geo_basic<T,M>&         domain,
        const Expr&                   expr,
        const form_option_type&       fopt);
    template <class Expr>
    void assembly (
        const band_basic<T,M>&        gh,
        const Expr&                   expr,
        const form_option_type&       fopt);

    // backward compat: named forms
    form_basic (const space_type& X, const space_type& Y,
	const std::string& name = "",
	const quadrature_option_type& qopt = quadrature_option_type(quadrature_option_type::max_family));

    form_basic (const space_type& X, const space_type& Y,
	const std::string& name,
	const field_basic<T,M>& weight, 
	const quadrature_option_type& qopt = quadrature_option_type(quadrature_option_type::max_family));

    template<class Function>
    form_basic (const space_type& X, const space_type& Y,
	const std::string& name,
	Function weight,
	const quadrature_option_type& qopt = quadrature_option_type(quadrature_option_type::max_family));

    form_basic (const space_type& X, const space_type& Y,
	const std::string& name,
	const geo_basic<T,M>& gamma,
	const quadrature_option_type& qopt = quadrature_option_type(quadrature_option_type::max_family));

    form_basic (const space_type& X, const space_type& Y,
	const std::string& name,
	const geo_basic<T,M>& gamma,
	const field_basic<T,M>& weight,
	const quadrature_option_type& qopt = quadrature_option_type(quadrature_option_type::max_family));

    template<class Function>
    form_basic (
        const space_type& X, 
        const space_type& Y,
        const std::string& name,
        const geo_basic<T,M>& gamma,
        Function weight,
	const quadrature_option_type& qopt = quadrature_option_type(quadrature_option_type::max_family));
protected: 
    // backward compat: named forms (cont.)
    template<class WeightFunction>
    void form_init (
                   const std::string&      name,
                   bool                    has_weight,
                   WeightFunction          weight,
                   const quadrature_option_type& qopt);
    template<class WeightFunction>
    void form_init_on_domain (
    		   const std::string&      name,
                   const geo_basic<T,M>&   gamma,
                   bool                    has_weight,
    		   WeightFunction          weight,
    		   const geo_basic<T,M>&   w_omega, // the domain where the fct weight is defined
    		   const quadrature_option_type& qopt);
};
template<class T, class M> form_basic<T,M> trans (const form_basic<T,M>& a);
template<class T, class M> field_basic<T,M> diag (const form_basic<T,M>& a);
template<class T, class M> form_basic<T,M>  diag (const field_basic<T,M>& dh);
typedef form_basic<Float,rheo_default_memory_model> form;
//>form:

// ------------ inline'd -----------------------------------

template<class T, class M>
inline
form_basic<T,M>::form_basic ()
: _X(), _Y(), _uu(), _ub(), _bu(), _bb()
{
}
template<class T, class M>
inline
form_basic<T,M>::form_basic (const form_basic<T,M>& a)
: _X(a._X), _Y(a._Y), _uu(a._uu), _ub(a._ub), _bu(a._bu), _bb(a._bb)
{
}
template<class T, class M>
inline
const typename form_basic<T,M>::space_type&
form_basic<T,M>::get_first_space() const
{
  return _X;
}
template<class T, class M>
inline
const typename form_basic<T,M>::space_type&
form_basic<T,M>::get_second_space() const
{
  return _Y;
}
template<class T, class M>
inline
const typename form_basic<T,M>::geo_type&
form_basic<T,M>::get_geo() const
{
  return _X.get_geo();
}
template<class T, class M>
inline
const communicator&
form_basic<T,M>::comm() const
{
  return get_geo().comm();
}
// ----------------
// linear albebra
// ----------------
template<class T, class M>
inline
form_basic<T,M>
form_basic<T,M>::operator+ (const form_basic<T,M>& b) const
{
  form_basic<T,M> c (get_first_space(), get_second_space());
  c._uu = _uu + b._uu;
  c._ub = _ub + b._ub;
  c._bu = _bu + b._bu;
  c._bb = _bb + b._bb;
  return c;
}
template<class T, class M>
inline
form_basic<T,M>
form_basic<T,M>::operator- (const form_basic<T,M>& b) const
{
  form_basic<T,M> c (get_first_space(), get_second_space());
  c._uu = _uu - b._uu;
  c._ub = _ub - b._ub;
  c._bu = _bu - b._bu;
  c._bb = _bb - b._bb;
  return c;
}
template<class T, class M>
inline
form_basic<T,M>
form_basic<T,M>::operator* (const form_basic<T,M>& b) const
{
  form_basic<T,M> c (b.get_first_space(), get_second_space());
  c._uu = _uu*b._uu + _ub*b._bu;
  c._ub = _uu*b._ub + _ub*b._bb;
  c._bu = _bu*b._uu + _bb*b._bu;
  c._bb = _bu*b._ub + _bb*b._bb;
  return c;
}
template<class T, class M>
inline
form_basic<T,M>&
form_basic<T,M>::operator*= (const T& lambda)
{
  _uu *= lambda;
  _ub *= lambda;
  _bu *= lambda;
  _bb *= lambda;
  return *this;
}
template<class T, class M>
inline
form_basic<T,M>
operator* (const T& lambda, const form_basic<T,M>& a)
{
  form_basic<T,M> b = a;
  b *= lambda;
  return b;
}
template<class T, class M>
inline
form_basic<T,M>
operator- (const form_basic<T,M>& a)
{
  return T(-1)*a;
}

}// namespace rheolef
# endif /* _RHEOLEF_FORM_H */