field_lazy_form_mult.h

Go to the documentation of this file.

# ifndef _RHEOLEF_FIELD_LAZY_FORM_MULT_H
# define _RHEOLEF_FIELD_LAZY_FORM_MULT_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*field_lazy and such
// AUTHOR: Pierre.Saramito@imag.fr
// DATE:   6 april 1920
//
// SUMMARY:
// 1. form_lazy*field_lazy -> field_lazy
// 2. form_lazy.trans_mult (field_lazy) -> field_lazy
//
#include "rheolef/field_lazy_node.h"
#include "rheolef/form.h"
 
namespace rheolef {
 
// -------------------------------------------------------------------
// 1. form_lazy*field_lazy -> field_lazy
// -------------------------------------------------------------------
namespace details {
 
template<class FormExpr, class FieldExpr>
class field_lazy_mult_form_rep {
public :
// definitions:
 
  using size_type = geo_element::size_type;
  using memory_type = typename FieldExpr::memory_type;
  using scalar_type = typename FieldExpr::scalar_type;
  using space_type  = typename FieldExpr::space_type;
  using geo_type    = typename FieldExpr::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
 
// allocator:
 
  field_lazy_mult_form_rep (const FormExpr& a_expr, const FieldExpr& u_expr);
 
// accessors:
 
  const geo_type&   get_geo()   const { return _a_expr.get_geo(); }
  const space_type& get_space() const { return _a_expr.get_test_space(); }
  bool  is_on_band()            const { return _u_expr.is_on_band(); }
  band_type get_band()          const { return _u_expr.get_band(); }
 
  void initialize (const geo_type& omega_K) const;
 
  void evaluate (
        const geo_type&            omega_K,
        const geo_element&         K,
              vector_element_type& uk) const;
// data:
protected:
  mutable FormExpr            _a_expr;
  mutable FieldExpr           _u_expr;
  mutable geo_type            _prev_omega_K;
  mutable size_type           _prev_K_dis_ie;
  mutable vector_element_type _prev_vk;
};
// inlined;
template<class FormExpr, class FieldExpr>
field_lazy_mult_form_rep<FormExpr,FieldExpr>::field_lazy_mult_form_rep (
  const FormExpr& a_expr,
  const FieldExpr& u_expr)
  : _a_expr(a_expr),
    _u_expr(u_expr),
    _prev_omega_K(),
    _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
    _prev_vk()
{
}
template<class FormExpr, class FieldExpr>
void
field_lazy_mult_form_rep<FormExpr,FieldExpr>::initialize (const geo_type& omega_K) const
{
  // TODO: subdomains e.g. robin
  check_macro (_a_expr.get_trial_space().get_geo() == _u_expr.get_geo(),
    "lazy_multiply: different domain not yet supported");
 
  check_macro (_a_expr.get_trial_space() == _u_expr.get_space(),
    "lazy_multiply: incompatible spaces \""
    <<_a_expr.get_trial_space().name()<<"\" and \""
    <<_u_expr.get_space().name()<<"\"");
 
  if (! _a_expr. get_test_space().get_constitution().have_compact_support_inside_element() ||
      ! _a_expr.get_trial_space().get_constitution().have_compact_support_inside_element()) {
    warning_macro("lazy_multiply: requires compact support inside elements (e.g. discontinuous or bubble)");
    warning_macro("lazy_multiply: form spaces was "
        << "[\"" <<_a_expr.get_trial_space().name()<<"\", \"" <<_a_expr.get_test_space().name()<<"\"]");
    fatal_macro("lazy_multiply: HINT: convert to \"form\" before to do the product");
  }
  _a_expr.initialize (omega_K);
  _u_expr.initialize (omega_K);
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max(); 
}
template<class FormExpr, class FieldExpr>
void
field_lazy_mult_form_rep<FormExpr,FieldExpr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& vk) const
{
  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
    trace_macro("amulx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
    vk = _prev_vk;
    return;
  }
  trace_macro("amulx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
  matrix_element_type ak;
  vector_element_type uk;
  _a_expr.evaluate (omega_K, K, ak);
  _u_expr.evaluate (omega_K, K, uk);
#ifdef _RHEOLEF_PARANO
  check_macro (ak.cols() == uk.size(), "a*u: invalid sizes");
#endif // _RHEOLEF_PARANO
  vk = ak*uk;
  _prev_omega_K  = omega_K;
  _prev_vk       = vk; // expensive to compute, so memorize it for common subexpressions
  _prev_K_dis_ie = K.dis_ie();
}
template<class FormExpr, class FieldExpr>
class field_lazy_mult_form: public smart_pointer_nocopy<field_lazy_mult_form_rep<FormExpr,FieldExpr>>,
                            public field_lazy_base     <field_lazy_mult_form    <FormExpr,FieldExpr>> {
public :
// definitions:
 
  using rep   = field_lazy_mult_form_rep<FormExpr,FieldExpr>;
  using base1 = smart_pointer_nocopy<rep>;
  using base2 = field_lazy_base<field_lazy_mult_form<FormExpr,FieldExpr>>;
  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 vector_element_type = typename rep::vector_element_type;
 
// allocator:
 
  field_lazy_mult_form (const FormExpr& a_expr, const FieldExpr& u_expr)
   : base1(new_macro(rep(a_expr,u_expr))),
     base2()
   {}
 
// accessors:
 
  const geo_type&   get_geo()   const { return base1::data().get_geo(); }
  const space_type& get_space() const { return base1::data().get_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 { base1::data().initialize (omega_K); }
 
   void evaluate (
        const geo_type&            omega_K,
        const geo_element&         K,
              vector_element_type& uk) const
        { base1::data().evaluate (omega_K, K, uk); }
};
// concept;
template<class FormExpr, class FieldExpr>
struct is_field_lazy <field_lazy_mult_form<FormExpr,FieldExpr> > : std::true_type {};
 
}// namespace details

template<class FormExpr, class FieldExpr,
    class Sfinae1 = typename std::enable_if<details:: is_form_lazy<FormExpr> ::value, FormExpr>::type,
    class Sfinae2 = typename std::enable_if<details::is_field_lazy<FieldExpr>::value, FieldExpr>::type>
details::field_lazy_mult_form<FormExpr,FieldExpr>
operator* (const FormExpr& a, const FieldExpr& u)
{
  return details::field_lazy_mult_form<FormExpr,FieldExpr> (a,u);
}

template<
  class FormExpr
 ,class Sfinae = typename std::enable_if<details::is_form_lazy<FormExpr>::value, FormExpr>::type
>
details::field_lazy_mult_form<
  FormExpr
 ,details::field_lazy_terminal_field<
    typename FormExpr::scalar_type
   ,typename FormExpr::memory_type
  >
>
operator* (
  const FormExpr& a,
  const field_basic<typename FormExpr::scalar_type,typename FormExpr::memory_type>& u)
{
  using wrap_type
    = details::field_lazy_terminal_field<
        typename FormExpr::scalar_type
       ,typename FormExpr::memory_type>;
  return details::field_lazy_mult_form<FormExpr,wrap_type> (a,wrap_type(u));
}
// -------------------------------------------------------------------
// 2. form_lazy.trans_mult (field_lazy) -> field_lazy
// -------------------------------------------------------------------
namespace details {
 
template<class FormExpr, class FieldExpr>
class field_lazy_trans_mult_form_rep {
public :
// definitions:
 
  using size_type   = geo_element::size_type;
  using memory_type = typename FieldExpr::memory_type;
  using scalar_type = typename FieldExpr::scalar_type;
  using space_type  = typename FieldExpr::space_type;
  using geo_type    = typename FieldExpr::geo_type;
  using float_type  = typename float_traits<scalar_type>::type;
  using band_type   = band_basic<float_type,memory_type>;
  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
 
// allocator:
 
  field_lazy_trans_mult_form_rep (const FormExpr& a_expr, const FieldExpr& u_expr);
 
// accessors:
 
  const geo_type&   get_geo()   const { return _a_expr.get_geo(); }
  const space_type& get_space() const { return _a_expr.get_trial_space(); }
  bool  is_on_band()            const { return _u_expr.is_on_band(); }
  band_type get_band()          const { return _u_expr.get_band(); }
 
  void initialize (const geo_type& omega_K) const;
 
  void evaluate (
        const geo_type&      omega_K,
        const geo_element&   K,
        vector_element_type& uk) const;
// data:
protected:
  mutable FormExpr            _a_expr;
  mutable FieldExpr           _u_expr;
  mutable geo_type            _prev_omega_K;
  mutable size_type           _prev_K_dis_ie;
  mutable vector_element_type _prev_vk;
};
// inlined;
template<class FormExpr, class FieldExpr>
field_lazy_trans_mult_form_rep<FormExpr,FieldExpr>::field_lazy_trans_mult_form_rep (
  const FormExpr& a_expr,
  const FieldExpr& u_expr)
  : _a_expr(a_expr),
    _u_expr(u_expr),
    _prev_omega_K(),
    _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
    _prev_vk()
{
}
template<class FormExpr, class FieldExpr>
void
field_lazy_trans_mult_form_rep<FormExpr,FieldExpr>::initialize (const geo_type& omega_K) const
{
  // TODO: subdomains e.g. robin
  check_macro (_a_expr.get_test_space().get_geo() == _u_expr.get_geo(),
    "lazy_trans_mult: different domain not yet supported");
 
  check_macro (_a_expr.get_test_space() == _u_expr.get_space(),
    "lazy_trans_mult: incompatible spaces \""
    <<_a_expr.get_test_space().name()<<"\" and \""
    <<_u_expr.get_space().name()<<"\"");
 
  if (! _a_expr. get_test_space().get_constitution().have_compact_support_inside_element() ||
      ! _a_expr.get_trial_space().get_constitution().have_compact_support_inside_element()) {
    warning_macro("lazy_trans_mult: requires compact support inside elements (e.g. discontinuous or bubble)");
    warning_macro("lazy_trans_mult: form spaces was "
        << "[\"" <<_a_expr.get_trial_space().name()<<"\", \"" <<_a_expr.get_test_space().name()<<"\"]");
    fatal_macro("lazy_trans_mult: HINT: convert to \"form\" before to do the product");
  }
  _a_expr.initialize (omega_K);
  _u_expr.initialize (omega_K);
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = std::numeric_limits<size_type>::max(); 
}
template<class FormExpr, class FieldExpr>
void
field_lazy_trans_mult_form_rep<FormExpr,FieldExpr>::evaluate (
  const geo_type&            omega_K,
  const geo_element&         K,
        vector_element_type& vk) const
{
  if (_prev_K_dis_ie == K.dis_ie()) {
    trace_macro("atansx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
    vk = _prev_vk;
    return;
  }
  trace_macro("atansx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
  matrix_element_type ak;
  vector_element_type uk;
  _a_expr.evaluate (omega_K, K, ak);
  _u_expr.evaluate (omega_K, K, uk);
#ifdef _RHEOLEF_PARANO
  check_macro (ak.rows() == uk.size(), "a*u: invalid sizes");
#endif // _RHEOLEF_PARANO
  vk = uk.transpose()*ak; 
  _prev_vk       = vk; // expensive to compute, so memorize it for common subexpressions
  _prev_omega_K  = omega_K;
  _prev_K_dis_ie = K.dis_ie();
}
template<class FormExpr, class FieldExpr>
class field_lazy_trans_mult_form: public smart_pointer_nocopy<field_lazy_trans_mult_form_rep<FormExpr,FieldExpr>>,
                                  public field_lazy_base     <field_lazy_trans_mult_form    <FormExpr,FieldExpr>> {
public :
// definitions:
 
  using rep   = field_lazy_trans_mult_form_rep<FormExpr,FieldExpr>;
  using base1 = smart_pointer_nocopy<rep>;
  using base2 = field_lazy_base<field_lazy_trans_mult_form<FormExpr,FieldExpr>>;
  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 vector_element_type = typename rep::vector_element_type;
 
// allocator:
 
  field_lazy_trans_mult_form (const FormExpr& a_expr, const FieldExpr& u_expr)
   : base1(new_macro(rep(a_expr,u_expr))),
     base2()
   {}
 
// accessors:
 
  const geo_type&   get_geo()   const { return base1::data().get_geo(); }
  const space_type& get_space() const { return base1::data().get_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 { base1::data().initialize (omega_K); }
 
   void evaluate (
        const geo_type&            omega_K,
        const geo_element&         K,
              vector_element_type& uk) const
        { base1::data().evaluate (omega_K, K, uk); }
};
// concept;
template<class FormExpr, class FieldExpr>
struct is_field_lazy <field_lazy_trans_mult_form<FormExpr,FieldExpr> > : std::true_type {};
 
}// namespace details
}// namespace rheolef
# endif /* _RHEOLEF_FIELD_LAZY_FORM_MULT_H */