librheolef-dev 6.5-1+b1 / usr / include / rheolef / integrate.h

#ifndef _RHEO_INTEGRATE_H
#define _RHEO_INTEGRATE_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/test.h"
#include "rheolef/field_vf_expr.h"
#include "rheolef/form_vf_expr.h"
#include "rheolef/form_option_type.h"
#include "rheolef/integrate_numeric.icc"
#include "rheolef/riesz.h" // for details::f_constant
namespace rheolef { 

/*Class:integrate
NAME: @code{integrate} - integrate a function or an expression
@findex  integrate
DESCRIPTION:
  Integrate an expression over a domain by using a quadrature formulae.
  There are three main usages of the integrate function, depending upon the 
  type of the expression.
  (i) When the expression is a numerical one, it leads to a numerical value.
  (ii) When the expression involves a symbolic test-function @pxref{test class},
  the result is a linear form, represented by the @code{field} class. 
  (iii) When the expression involves both  symbolic trial- and test-functions @pxref{test class},
  the result is a bilinear form, represented by the @code{field} class. 
SYNOPSYS:
@example
 template <class T, class M, class Expr>
 T integrate (const geo_basic<T,M>& omega, Expr expr,
   quadrature_option_type qopt)

 template <class T, class M, class Expr>
 field integrate (const geo_basic<T,M>& omega, VFExpr expr,
   quadrature_option_type qopt)

 template <class T, class M, class Expr>
 form integrate (const geo_basic<T,M>& omega, VFExpr expr,
   form_option_type fopt)
@end example
EXAMPLE:
@noindent
@example
  Float f (const point& x);
  ...
  quadrature_option_type qopt;
  Float value = integrate (omega, f, qopt);
  field lh = integrate (omega, f*v, qopt);
@end example
  The last argument specifies the quadrature formulae used
  for the computation of the integral.
  The expression can be any function, classs-function or 
  any linear or nonlinear field expression @pxref{field class}.
DEFAULT ARGUMENTS:
  In the case of a linear form, the domain is optional: by default it is
  the full domain definition of the test function.
@example
  field l1h = integrate (f*v, qopt);
@end example
  When the integration is perfomed on a subdomain, this subdomain
  simply replace the first argument and a domain name could also be used:
@example
  field l2h = integrate (omega["boundary"], f*v, qopt);
  field l3h = integrate ("boundary", f*v, qopt);
@end example
  The quadrature formulae is required, except when a test and/or trial
  function is provided in the expression to integrate.
  In that case, the quadrature formulae is deduced from the space
  containing the test (or trial) function. 
  When a test function is suppied, let k be its polynomial degree.
  Then the default quadrature is
  choosen to be exact at least for 2*k+1 polynoms.
  When both a test and trial functions are suppied, let k1 and k2 be their polynomial degrees.
  Then the default quadrature is
  choosen to be exact at least for k1+k2+1 polynoms.
  Also, when the expression is a constant, the quadrature function is optional:
  in that case, the constant is also optional and the following call:
@example
  Float meas = integrate (omega);
@end example
  is valid and returns the measure of the domain.
End: */

// ----------------------------------------------
// variational integration
// ----------------------------------------------
template <class T, class M>
inline
field_basic<T,M>
integrate (const geo_basic<T,M>& domain, const test_basic<T,M,details::vf_tag_01>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (domain, expr, qopt);
  return lh;
}
template <class T, class M, class Expr>
inline
field_basic<T,M>
integrate (const geo_basic<T,M>& domain, const field_vf_expr<Expr>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (domain, expr, qopt);
  return lh;
}
template <class T, class M, class Expr>
inline
form_basic<T, M>
integrate (const geo_basic<T,M>& domain, const form_vf_expr<Expr>& expr,
	const form_option_type& fopt = form_option_type())
{
  form_basic<T,M> a;
  a.assembly (domain, expr, fopt);
  return a;
}
// ----------------------------------------------
// variational integration: missing domain
// ----------------------------------------------
template <class T, class M>
inline
field_basic<T,M>
integrate (const test_basic<T,M,details::vf_tag_01>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (expr.get_vf_space().get_geo(), expr, qopt);
  return lh;
}
template <class Expr>
inline
field_basic<typename field_vf_expr<Expr>::scalar_type, typename field_vf_expr<Expr>::memory_type>
integrate (const field_vf_expr<Expr>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  typedef typename field_vf_expr<Expr>::scalar_type T;
  typedef typename field_vf_expr<Expr>::memory_type M;
  field_basic<T,M> lh;
  lh.assembly (expr.get_vf_space().get_geo(), expr, qopt);
  return lh;
}
template <class Expr>
inline
form_basic<typename form_vf_expr<Expr>::scalar_type, typename form_vf_expr<Expr>::memory_type>
integrate (const form_vf_expr<Expr>& expr,
	const form_option_type& fopt = form_option_type())
{
  // TODO: prendre le min de expr.get_test_space().get_geo() expr.get_trial_space().get_geo()
  typedef typename field_vf_expr<Expr>::scalar_type T;
  typedef typename field_vf_expr<Expr>::memory_type M;
  form_basic<T,M> a;
  a.assembly (expr.get_test_space().get_geo(), expr, fopt);
  return a;
}
// ----------------------------------------------
// variational integration: subdomain by its name
// ----------------------------------------------
template <class T, class M>
inline
field_basic<T,M>
integrate (const std::string& domname, const test_basic<T,M,details::vf_tag_01>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (expr.get_vf_space().get_geo()[domname], expr, qopt);
  return lh;
}
template <class Expr>
inline
field_basic<typename field_vf_expr<Expr>::scalar_type, typename field_vf_expr<Expr>::memory_type>
integrate (const std::string& domname, const field_vf_expr<Expr>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  typedef typename field_vf_expr<Expr>::scalar_type T;
  typedef typename field_vf_expr<Expr>::memory_type M;
  field_basic<T,M> lh;
  lh.assembly (expr.get_vf_space().get_geo()[domname], expr, qopt);
  return lh;
}
template <class Expr>
inline
form_basic<typename form_vf_expr<Expr>::scalar_type, typename form_vf_expr<Expr>::memory_type>
integrate (const std::string& domname, const form_vf_expr<Expr>& expr,
	const form_option_type& fopt = form_option_type())
{
  // TODO: prendre le max de expr.get_test_space().get_geo() expr.get_trial_space().get_geo() puis [domain]
  typedef typename field_vf_expr<Expr>::scalar_type T;
  typedef typename field_vf_expr<Expr>::memory_type M;
  form_basic<T,M> a;
  a.assembly (expr.get_test_space().get_geo()[domname], expr, fopt);
  return a;
}
// ----------------------------------------------
// variational integration: on a band
// ----------------------------------------------
template <class T, class M>
inline
field_basic<T,M>
integrate (const band_basic<T,M>& gh, const test_basic<T,M,details::vf_tag_01>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (gh, expr, qopt);
  return lh;
}
template <class T, class M, class Expr>
inline
field_basic<T,M>
integrate (const band_basic<T,M>& gh, const field_vf_expr<Expr>& expr,
	const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (gh, expr, qopt);
  return lh;
}
template <class T, class M, class Expr>
inline
form_basic<T, M>
integrate (const band_basic<T,M>& gh, const form_vf_expr<Expr>& expr,
	const form_option_type& fopt = form_option_type())
{
  form_basic<T,M> a;
  a.assembly (gh, expr, fopt);
  return a;
}
// ----------------------------------------------
// numeric integration
// ----------------------------------------------
template <class T, class M, class Function>
T integrate (const geo_basic<T,M>& dom, const field_functor<Function,T>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_function<Function>  arg_t;
  typedef field_nonlinear_expr<arg_t>             expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f.get_ref())), qopt);
}
// Float r = integrate(omega); ICI
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, const T& f = 1, const quadrature_option_type& qopt = quadrature_option_type())
{
  typedef details::f_constant<point_basic<T>,T> f_cte_t;
  typedef field_expr_terminal_function<f_cte_t> arg_t;
  typedef field_nonlinear_expr<arg_t>           expr_t;
  f_cte_t f_cte (f);
  return details::integrate_numeric (dom, expr_t(arg_t(f_cte)), qopt);
}
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, T (*f)(const point_basic<T>&), const quadrature_option_type& qopt)
{
  typedef T (*fun_t)(const point_basic<T>&);
  typedef field_expr_terminal_function<fun_t>     arg_t;
  typedef field_nonlinear_expr<arg_t>             expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, const field_basic<T,M>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_field<T,M>  arg_t;
  typedef field_nonlinear_expr<arg_t>     expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, const field_indirect<T,M>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_field<T,M>  arg_t;
  typedef field_nonlinear_expr<arg_t>     expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, const field_indirect_const<T,M>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_field<T,M>  arg_t;
  typedef field_nonlinear_expr<arg_t>     expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, const field_component<T,M>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_field<T,M>  arg_t;
  typedef field_nonlinear_expr<arg_t>     expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M>
T integrate (const geo_basic<T,M>& dom, const field_component_const<T,M>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_field<T,M>  arg_t;
  typedef field_nonlinear_expr<arg_t>     expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M, class E>
T integrate (const geo_basic<T,M>& dom, const field_expr<E>& f, const quadrature_option_type& qopt)
{
  typedef field_expr_terminal_field<T,M>  arg_t;
  typedef field_nonlinear_expr<arg_t>     expr_t;
  return details::integrate_numeric (dom, expr_t(arg_t(f)), qopt);
}
template <class T, class M, class E>
T integrate (const geo_basic<T,M>& dom, const field_nonlinear_expr<E>& f, const quadrature_option_type& qopt)
{
  return details::integrate_numeric (dom, f, qopt);
}

}// namespace rheolef
#endif // _RHEO_INTEGRATE_H