form_lazy_terminal.h

Go to the documentation of this file.

# ifndef _RHEOLEF_FORM_LAZY_TERMINAL_H
# define _RHEOLEF_FORM_LAZY_TERMINAL_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
// 
// =========================================================================
// form_lazy_expr = result of a lazy integrate() that are combined in exprs
// AUTHOR: Pierre.Saramito@imag.fr
// DATE:   30 march 1920
//
// SUMMARY:
// 1. concept & base class
// 2. terminal
//    2.1. integrate
//    2.2. integrate on a band
// see also "form_lazy_expr.h"
 
#include "rheolef/form.h"
#include "rheolef/field_lazy_form_mult.h"
 
namespace rheolef {
 
// -------------------------------------------------------------------
// 1. concept & base class
// -------------------------------------------------------------------
namespace details {
// Define a trait type for detecting form expression valid arguments
//     template<class T> struct is_form_lazy: std::false_type {};
// => should be defined in form.h for compilation reason,
//    otherwise Sfinae is always failed in form.h
 
// Define a base class
template<class Derived>
class form_lazy_base {
public:
// a.trans_mult(u)
 
  template<class FieldExpr>
  typename std::enable_if<
    details::is_field_lazy<FieldExpr>::value
   ,field_lazy_trans_mult_form<Derived,FieldExpr>
  >::type
  trans_mult (const FieldExpr& u_expr) const {
    return field_lazy_trans_mult_form<Derived,FieldExpr>(derived(), u_expr);
  }
  // TODO: how to obtain T and M from Derived ?
  //       typename Derived::scalar_type => compilation failed
  //       see eg EigenBase in eigen library
  template<class T, class M>
  field_lazy_trans_mult_form<Derived,field_lazy_terminal_field<T,M>>
  trans_mult (const field_basic<T,M>& uh) const {
    return field_lazy_trans_mult_form<Derived,field_lazy_terminal_field<T,M>>
             (derived(), field_lazy_terminal_field<T,M>(uh));
  }
 
protected:
  Derived&       derived()       { return *static_cast<      Derived*>(this); }
  const Derived& derived() const { return *static_cast<const Derived*>(this); }
};
 
} // namespace details
// -------------------------------------------------------------------
// 2. terminal
// -------------------------------------------------------------------
// 2.1. integrate
// -------------------------------------------------------------------
// form_lazy_terminal_integrate is returned by the integrate() function
// It contains the domain of integration
// and integration options (quadrature)
namespace details {
 
template<class Expr>
class form_lazy_terminal_integrate_rep {
public :
// definitions:
 
  using size_type   = geo_element::size_type;
  using memory_type = typename Expr::memory_type;
  using scalar_type = typename Expr::scalar_type;
  using space_type  = typename Expr::space_type;
  using geo_type    = typename Expr::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
 
// allocator:
 
  form_lazy_terminal_integrate_rep (
    const geo_type&         domain,
    const Expr&             expr,
    const integrate_option& iopt)
  : _domain(domain),
    _expr(expr),
    _iopt(iopt),
    _prev_omega_K(),
    _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
    _prev_ak()
  {}
 
// accessors:
 
  const geo_type&   get_geo()         const { return _domain; }
  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
  const space_type& get_test_space()  const { return _expr.get_test_space(); }
  bool  is_on_band()                  const { return false; }
  band_type get_band()                const { return band_type(); }
 
  void initialize (const geo_type& omega_K) const;
  bool is_diagonal() const { return false; }
 
  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            matrix_element_type& ak) const;
// data:
protected:
  geo_type                     _domain;
  mutable Expr                 _expr;
  mutable integrate_option     _iopt;
  mutable geo_type             _prev_omega_K;
  mutable size_type            _prev_K_dis_ie;
  mutable matrix_element_type  _prev_ak;
};
// inlined;
template<class Expr>
void
form_lazy_terminal_integrate_rep<Expr>::initialize (const geo_type& omega_K) const
{
  // TODO: sub-domain
  // _domain could be a subdomain of omega_K:
  // in that case evaluate(omega_K,K) should return a zero ak matrix
  // when K do not belongs to _domain
  _iopt._is_on_interface          = false;
  _iopt._is_inside_on_local_sides = false;
  _expr.initialize (_domain, _iopt);
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
}
template<class Expr>
void
form_lazy_terminal_integrate_rep<Expr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        matrix_element_type& ak) const
{
  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
    trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
    ak = _prev_ak;
    return;
  }
  _expr.evaluate (omega_K, K, ak);
  _prev_ak       = ak; // expensive to compute, so memorize it for common subexpressions
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = K.dis_ie();
  trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
}
 
template<class Expr>
class form_lazy_terminal_integrate: public smart_pointer_nocopy<form_lazy_terminal_integrate_rep<Expr>>,
                                    public form_lazy_base      <form_lazy_terminal_integrate<Expr>> {
public :
// definitions:
 
  using rep    = form_lazy_terminal_integrate_rep<Expr>;
  using base1  = smart_pointer_nocopy<rep>;
  using base2  = form_lazy_base<form_lazy_terminal_integrate<Expr>>;
  using size_type   = typename rep::size_type;
  using memory_type = typename rep::memory_type;
  using scalar_type = typename rep::scalar_type;
  using space_type  = typename rep::space_type;
  using geo_type    = typename rep::geo_type;
  using band_type   = typename rep::band_type;
  using matrix_element_type = typename rep::matrix_element_type;
 
// allocator:
 
  form_lazy_terminal_integrate (
    const geo_type&         domain,
    const Expr&             expr,
    const integrate_option& iopt)
   : base1(new_macro(rep(domain,expr,iopt))),
     base2()
   {}
 
// accessors:
 
  const geo_type&   get_geo()         const { return base1::data().get_geo(); }
  const space_type& get_trial_space() const { return base1::data().get_trial_space(); }
  const space_type& get_test_space()  const { return base1::data().get_test_space(); }
  bool  is_on_band()                  const { return base1::data().is_on_band(); }
  band_type get_band()                const { return base1::data().get_band(); }
 
  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
  bool is_diagonal() const                        { return base1::data().is_diagonal(); }
 
  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            matrix_element_type& ak) const
      { base1::data().evaluate (omega_K, K, ak); }
};
// concept;
template<class Expr> struct is_form_lazy <form_lazy_terminal_integrate<Expr> > : std::true_type {};
 
}// namespace details
// ----------------------------------------------
// 2.2. integrate on a band
// ----------------------------------------------
namespace details {
 
template<class Expr>
class form_lazy_terminal_integrate_band_rep {
public :
// definitions:
 
  using size_type   = geo_element::size_type;
  using memory_type = typename Expr::memory_type;
  using scalar_type = typename Expr::scalar_type;
  using space_type  = typename Expr::space_type;
  using geo_type    = typename Expr::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
 
// allocators:
 
  form_lazy_terminal_integrate_band_rep (
    const band_type&        gh,
    const Expr&             expr,
    const integrate_option& iopt);
 
// accessors:
 
  const geo_type&   get_geo()         const { return _gh.level_set(); }
  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
  const space_type& get_test_space()  const { return _expr.get_test_space(); }
  bool  is_on_band()                  const { return true; }
  band_type get_band()                const { return _gh; }
 
  void initialize (const geo_type& omega_K) const;
 
  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            matrix_element_type& ak) const;
// data:
protected:
  band_type                    _gh;
  mutable Expr                 _expr;
  mutable integrate_option     _iopt;
  mutable geo_type             _prev_omega_K;
  mutable size_type            _prev_K_dis_ie;
  mutable matrix_element_type  _prev_ak;
};
// inlined;
template<class Expr>
form_lazy_terminal_integrate_band_rep<Expr>::form_lazy_terminal_integrate_band_rep (
  const band_type&        gh,
  const Expr&             expr,
  const integrate_option& iopt)
: _gh(gh),
  _expr(expr),
  _iopt(iopt),
  _prev_omega_K(),
  _prev_K_dis_ie(std::numeric_limits<size_type>::max()), 
  _prev_ak()
{
}
template<class Expr>
void
form_lazy_terminal_integrate_band_rep<Expr>::initialize (const geo_type& omega_K) const
{
  const geo_type& band  = _gh.band(); // the 3D intersected tetras
  const geo_type& bgd_omega = get_test_space().get_constitution().get_background_geo();
  check_macro (omega_K == get_geo(),
        "integrate on band: unexpected integration domain");
  check_macro (band.get_background_geo() == bgd_omega,
        "do_integrate: incompatible integration domain "<<_gh.level_set().name() << " and test function based domain "
        << bgd_omega.name());
  _expr.initialize (_gh, _iopt);
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
}
template<class Expr>
void
form_lazy_terminal_integrate_band_rep<Expr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        matrix_element_type& ak) const
{
  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
    trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
    ak = _prev_ak;
    return;
  }
  // TODO: memorize _prev
  _expr.evaluate (omega_K, K, ak);
  _prev_ak       = ak; // expensive to compute, so memorize it for common subexpressions
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = K.dis_ie();
}
template<class Expr>
class form_lazy_terminal_integrate_band: public smart_pointer_nocopy<form_lazy_terminal_integrate_band_rep<Expr>>,
                                          public form_lazy_base     <form_lazy_terminal_integrate_band<Expr>> {
public :
// definitions:
 
  using rep    = form_lazy_terminal_integrate_band_rep<Expr>;
  using base1  = smart_pointer_nocopy<rep>;
  using base2  = form_lazy_base<form_lazy_terminal_integrate_band<Expr>>;
  using size_type   = typename rep::size_type;
  using memory_type = typename rep::memory_type;
  using scalar_type = typename rep::scalar_type;
  using space_type  = typename rep::space_type;
  using geo_type    = typename rep::geo_type;
  using band_type   = typename rep::band_type;
  using matrix_element_type = typename rep::matrix_element_type;
 
// allocator:
 
  form_lazy_terminal_integrate_band (
    const band_type&        gh, 
    const Expr&             expr,
    const integrate_option& iopt)
   : base1(new_macro(rep(gh,expr,iopt))),
     base2()
   {}
 
// accessors:
 
  const geo_type&   get_geo()         const { return base1::data().get_geo(); }
  const space_type& get_trial_space() const { return base1::data().get_trial_space(); }
  const space_type& get_test_space()  const { return base1::data().get_test_space(); }
  bool  is_on_band()                  const { return base1::data().is_on_band(); }
  band_type get_band()                const { return base1::data().get_band(); }
 
  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
 
  void evaluate (
      const geo_type&            omega_K,
      const geo_element&         K,
            matrix_element_type& ak) const
      { base1::data().evaluate (omega_K, K, ak); }
};
// concept;
template<class Expr> struct is_form_lazy <form_lazy_terminal_integrate_band<Expr> > : std::true_type {};
 
}// namespace details
 
template <class Expr>
typename
std::enable_if<
  details::is_form_expr_quadrature_arg<Expr>::value
 ,details::form_lazy_terminal_integrate_band <Expr>
>::type
lazy_integrate (
  const band_basic<typename float_traits<typename Expr::scalar_type>::type,
                   typename Expr::memory_type>& gh, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
  return details::form_lazy_terminal_integrate_band<Expr> (gh, expr, iopt);
}
template <class Expr>
typename
std::enable_if<
  details::is_form_expr_v2_variational_arg<Expr>::value
 ,details::form_lazy_terminal_integrate_band <details::form_expr_quadrature_on_element<Expr>>
>::type
lazy_integrate (
  const band_basic<typename float_traits<typename Expr::scalar_type>::type,
                   typename Expr::memory_type>& gh, 
  const Expr& expr,
  const integrate_option& iopt = integrate_option())
{
 
  details::form_expr_quadrature_on_element<Expr> expr_quad(expr);
  return lazy_integrate (gh, expr_quad, iopt);
}
 
}// namespace rheolef
# endif /* _RHEOLEF_FORM_LAZY_TERMINAL_H */