form.h

Go to the documentation of this file.

# 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
// 
// =========================================================================
// AUTHOR: Pierre.Saramito@imag.fr
// DATE:   2 july 1997
 
namespace rheolef {
} // namespace rheolef
 
#include "rheolef/csr.h"
#include "rheolef/field.h"
 
namespace rheolef { 
 
namespace details { 
// 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;
 
// Define traits for un-assembled forms:
template<class Expr> struct is_form_lazy: std::false_type {};
} // namespace details
 
// forward declaration:
template <class T, class M> class band_basic;
 
 
// [verbatim_form_basic]
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>&);
    form_basic<T,M>& operator= (const form_basic<T,M>&);
 
    template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
    form_basic (const Expr&);
 
    template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
    form_basic<T,M>& operator= (const Expr&);
 
// allocators from initializer list (c++ 2011):
 
    form_basic (const std::initializer_list<details::form_concat_value<T,M> >& init_list);
    form_basic (const std::initializer_list<details::form_concat_line <T,M> >& init_list);
 
// accessors:
 
    const space_type& get_first_space() const;
    const space_type& get_second_space() const;
    const geo_type&   get_geo() const;
    bool is_symmetric() const;
    void set_symmetry (bool is_symm = true) const;
    bool is_definite_positive() const;
    void set_definite_positive (bool is_dp = true) const;
    bool is_symmetric_definite_positive() const;
    void set_symmetric_definite_positive() 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; }
// [verbatim_form_basic]
// internals:
    int constraint_process_rank() const { return uu().constraint_process_rank(); }
// 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: 
    // assembly from an un-assembled form_lazy expr
    template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
    void convert_from_form_lazy (const Expr& expr);
 
    // with vf expression arg
    template <class Expr>
    void do_integrate_internal (
        const geo_basic<T,M>&         dom,
        const geo_basic<T,M>&         band,
        const band_basic<T,M>&        gh,
        const Expr&                   expr,
        const integrate_option&       fopt,
        bool                          is_on_band);
    template <class Expr>
    void do_integrate (
        const geo_basic<T,M>&         domain,
        const Expr&                   expr,
        const integrate_option&       fopt);
    template <class Expr>
    void do_integrate (
        const band_basic<T,M>&        gh,
        const Expr&                   expr,
        const integrate_option&       fopt);
 
    // backward compat: named forms
    form_basic (const space_type& X, const space_type& Y,
        const std::string& name = "",
        const quadrature_option& qopt = quadrature_option());
 
    form_basic (const space_type& X, const space_type& Y,
        const std::string& name,
        const field_basic<T,M>& weight, 
        const quadrature_option& qopt = quadrature_option());
 
    template<class Function>
    form_basic (const space_type& X, const space_type& Y,
        const std::string& name,
        Function weight,
        const quadrature_option& qopt = quadrature_option());
 
    form_basic (const space_type& X, const space_type& Y,
        const std::string& name,
        const geo_basic<T,M>& gamma,
        const quadrature_option& qopt = quadrature_option());
 
    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& qopt = quadrature_option());
 
    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& qopt = quadrature_option());
protected: 
    // backward compat: named forms (cont.)
    template<class WeightFunction>
    void form_init (
                   const std::string&      name,
                   bool                    has_weight,
                   WeightFunction          weight,
                   const quadrature_option& 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& qopt);
// [verbatim_form_basic_cont]
};
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);
// [verbatim_form_basic_cont]
 
// [verbatim_form]
typedef form_basic<Float,rheo_default_memory_model> form;
// [verbatim_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
form_basic<T,M>& 
form_basic<T,M>::operator= (const form_basic<T,M>& a)
{
  _X.operator=  (a._X);
  _Y.operator=  (a._Y);
  _uu.operator= (a._uu);
  _ub.operator= (a._ub);
  _bu.operator= (a._bu);
  _bb.operator= (a._bb);
  return *this;
}
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;
}
template<class T, class M>
inline
bool
form_basic<T,M>::is_symmetric() const
{
  return _uu.is_symmetric() && _bb.is_symmetric();
}
template<class T, class M>
inline
void
form_basic<T,M>::set_symmetry (bool is_symm) const
{
  _uu.set_symmetry(is_symm);
  _bb.set_symmetry(is_symm);
}
template<class T, class M>
inline
bool
form_basic<T,M>::is_definite_positive() const
{
  return _uu.is_definite_positive() && _bb.is_definite_positive();
}
template<class T, class M>
inline
void
form_basic<T,M>::set_definite_positive (bool is_dp) const
{
  _uu.set_definite_positive(is_dp);
  _bb.set_definite_positive(is_dp);
}
template<class T, class M>
inline
bool
form_basic<T,M>::is_symmetric_definite_positive() const
{
  return is_symmetric() && is_definite_positive();
}
template<class T, class M>
inline
void
form_basic<T,M>::set_symmetric_definite_positive() const
{
  set_symmetry();
  set_definite_positive();
}
 
}// namespace rheolef
# endif /* _RHEOLEF_FORM_H */