librheolef-dev 6.7-6 / 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
/// 
/// =========================================================================
//
// SUMMARY:
// 1. numeric integration
//    1.1. general integration of a nonlinear expression
//    1.2. measure of the domain
//    1.3. when the valued result type is undetermined
// 2. field-result integration of a variational expression
//    2.1. general call
//    2.2. missing domain
//    2.3. subdomain by its name
// 3. form-result integration of a variational expression
//    3.1. general call
//    3.2. missing domain
//    3.3. subdomain by its name
// 4. variational integration: on a band
//
#include "rheolef/field_expr_v2_nonlinear.h"
#include "rheolef/field_expr_v2_variational.h"
#include "rheolef/form_expr_v2_variational.h"

#include "rheolef/field_expr_v2_value_assembly.h"
#include "rheolef/field_vf_assembly.h"
#include "rheolef/form_vf_assembly.h"

#include "rheolef/functor.h" // used to convert functions to functors

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
 Float integrate (geo domain);
 Float integrate (geo domain, quadrature_option_type qopt);
 Value integrate (geo domain, Expression, quadrature_option_type qopt);

 field integrate (Expression);
 field integrate (Expression, quadrature_option_type qopt);
 field integrate (geo domain, Expression);
 field integrate (geo domain, Expression, quadrature_option_type qopt);

 form integrate (Expression);
 form integrate (Expression, form_option_type qopt);
 form integrate (geo domain, Expression);
 form integrate (geo domain, Expression, form_option_type qopt);
@end example
EXAMPLE:
@noindent
  For computing the measure of a domain:
@example
  Float meas_omega = integrate (omega);
@end example
  For computing the integral of a function:
@example
  Float f (const point& x);
  ...
  quadrature_option_type qopt;
  qopt.set_order (3);
  Float int_f = integrate (omega, f, qopt);
@end example
  The last argument specifies the quadrature formulae 
  (see @ref{quadrature_option_type class})
  used for the computation of the integral.
  The function can be replaced by any field-valued expression (see @ref{field class}).
  For computing a right-hand-side of a variational formulation
  with the previous function @code{f}:
@example
  space Xh (omega, "P1");
  test v (Xh);
  field lh = integrate (f*v);
@end example
  For computing a bilinear form:
@example
  trial u (Xh);
  test v (Xh);
  form m = integrate (u*v);
@end example
  The expression @code{u*v} can be replaced by any bilinear expression (see @ref{field class}).
DEFAULT ARGUMENTS:
  In the case of a linear or bilinear form, the domain is optional: by default it is
  the full domain definition of the test function.
  Also, the quadrature formulae is optional: by default, its order
  is @code{2*k+1} where @code{k} is the polynomial degree of the
  @code{Xh} space associated to the test function @code{v}.
  When both a test @code{u} and trial @code{v} 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.
  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);
  field l3h = integrate ("boundary", f*v);
@end example
  For convenience, only the domain name can be supplied.
End: */

// ---------------------------------------------------
// 1. numeric integration
// ---------------------------------------------------
// 1.1. general integration of a nonlinear expression
// ---------------------------------------------------
template <class T, class M, class Expr,
    class Result = typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>
inline
typename std::enable_if<
     details::is_field_expr_v2_nonlinear_arg<Expr>::value
  && ! is_undeterminated<Result>::value,
  Result
>::type
integrate (const geo_basic<T,M>& omega, const Expr& expr, const quadrature_option_type& qopt,
	   Result dummy = Result())
{
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type  wrap_t;
  if (omega.map_dimension() < omega.get_background_geo().map_dimension()) {
    omega.get_background_geo().neighbour_guard();
  }
  Result result(0);
  field_expr_v2_value_assembly (omega, wrap_t(expr), qopt, result);
  return result;
}
// ---------------------------------------------------
// 1.2. measure of the domain
// ---------------------------------------------------
template <class T, class M>
T
integrate (const geo_basic<T,M>& omega, quadrature_option_type&& qopt = quadrature_option_type())
{
  if (qopt.get_order() == std::numeric_limits<quadrature_option_type::size_type>::max()) {
    qopt.set_order(0);
  }
  details::f_constant <point_basic<T>,T> one(1);
  return integrate (omega, one, qopt);
}
// ---------------------------------------------------
// 1.3. when the valued result type is undetermined
// ---------------------------------------------------
// TODO: return a overdetermined<T> value that is an union of all possibilities with a valued_tag
template<class T, class M, class Expr>
inline
typename std::enable_if<
     details::is_field_expr_v2_nonlinear_arg<Expr>::value
  && is_undeterminated<typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>::value,
  typename scalar_traits<typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>::type
>::type
integrate (const geo_basic<T,M>& omega, const Expr& expr, const quadrature_option_type& qopt)
{
  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type undef_t;
  typedef typename scalar_traits<undef_t>::type scalar_type;
  switch (expr.valued_tag()) {
    case space_constant::scalar: {
        return integrate (omega, expr, qopt, scalar_type());
    }
    // others type: problem on how to return a run-type type ?
    // TODO: return an overdetermined union type that convert to one of scalar, point, tensor, etc ?
    default:
        warning_macro ("Expr="<<pretty_typename_macro(Expr));
        error_macro ("integrate: not yet for `"
	  << space_constant::valued_name (expr.valued_tag())
          << "' valued expression");
        return 0;
  }
}
// -------------------------------------------------------
// 2. field-result integration of a variational expression
// -------------------------------------------------------
// 2.1. general call
// -------------------------------------------------------
template <class T, class M, class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,field_basic<T,M>
>::type
integrate (
  const geo_basic<T,M>& domain, 
  const Expr& expr,
  const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic<T,M> lh;
  lh.assembly (domain, expr, qopt);
  return lh;
}
// ----------------------------------------------
// 2.2. missing domain
// ----------------------------------------------
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const Expr& expr,
  const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic <typename Expr::scalar_type, typename Expr::memory_type> lh;
  lh.assembly (expr.get_vf_space().get_geo(), expr, qopt);
  return lh;
}
// ----------------------------------------------
// 2.3. subdomain by its name
// ----------------------------------------------
template <class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const std::string& domname, 
  const Expr& expr,
  const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic <typename Expr::scalar_type, typename Expr::memory_type> lh;
  lh.assembly (expr.get_vf_space().get_geo()[domname], expr, qopt);
  return lh;
}
// -------------------------------------------------------
// 3. form-result integration of a variational expression
// -------------------------------------------------------
// 3.1. general call
// -------------------------------------------------------
template <class T, class M, class Expr>
inline
typename
std::enable_if<
  details::is_form_expr_v2_variational_arg<Expr>::value
 ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const geo_basic<T,M>& domain, 
  const Expr& expr,
  const form_option_type& fopt = form_option_type())
{
  form_basic<T,M> a;
  a.assembly (domain, expr, fopt);
  return a;
}
// ----------------------------------------------
// 3.2. missing domain
// ----------------------------------------------
template <class Expr>
inline
typename
std::enable_if<
  details::is_form_expr_v2_variational_arg<Expr>::value
 ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const Expr& expr,
  const form_option_type& fopt = form_option_type())
{
  form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
  a.assembly (expr.get_test_space().get_geo(), expr, fopt);
  return a;
}
// ----------------------------------------------
// 3.3. subdomain by its name
// ----------------------------------------------
template <class Expr>
inline
typename
std::enable_if<
  details::is_form_expr_v2_variational_arg<Expr>::value
 ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const std::string& domname, 
  const Expr& expr,
  const form_option_type& fopt = form_option_type())
{
  form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
  a.assembly (expr.get_test_space().get_geo()[domname], expr, fopt);
  return a;
}
// ----------------------------------------------
// 4. variational integration: on a band
// ----------------------------------------------
template <class T, class M, class Expr>
inline
typename
std::enable_if<
  details::is_field_expr_v2_variational_arg<Expr>::value
 ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const band_basic<T,M>& gh, 
  const Expr& expr,
  const quadrature_option_type& qopt = quadrature_option_type())
{
  field_basic <typename Expr::scalar_type, typename Expr::memory_type> lh;
  lh.assembly (gh, expr, qopt);
  return lh;
}
template <class T, class M, class Expr>
inline
typename
std::enable_if<
  details::is_form_expr_v2_variational_arg<Expr>::value
 ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
>::type
integrate (
  const band_basic<T,M>& gh, 
  const Expr& expr,
  const form_option_type& fopt = form_option_type())
{
  form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
  a.assembly (gh, expr, fopt);
  return a;
}

}// namespace rheolef
#endif // _RHEO_INTEGRATE_H