librheolef-dev 6.7-6 / usr / include / rheolef / riesz.h

#ifndef _RHEO_RIESZ_H
#define _RHEO_RIESZ_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/field.h"
#include "rheolef/band.h"
#include "rheolef/test.h"
#include "rheolef/integrate.h"

namespace rheolef { 

// --------------------------------------------------------------------------
// implementation
// --------------------------------------------------------------------------
namespace details { 
// --------------------------------------------------------------------------
// generic implementation
// --------------------------------------------------------------------------
// scalar-valued function
template <class T, class M, class Geo, class Function>
field_basic<T,M>
riesz_function_internal (
    const space_basic<T,M>&         Xh,
    const Function&                 f,
    const quadrature_option_type&   qopt,
    const Geo&                      dom,
    std::false_type)
{
  test_basic<T,M,vf_tag_01> v (Xh);
  return integrate (dom, f*v, qopt);
}
// vector-valued function
template <class T, class M, class Geo, class Function>
field_basic<T,M>
riesz_function_internal (
    const space_basic<T,M>&         Xh,
    const Function&                 f,
    const quadrature_option_type&   qopt,
    const Geo&                      dom,
    std::true_type)
{
  test_basic<T,M,vf_tag_01> v (Xh);
  return integrate (dom, dot(f,v), qopt);
}
// --------------------------------------------------------------------------
// switch between compile-time return type:
// --------------------------------------------------------------------------
template <class T, class M, class Geo, class Function>
field_basic<T,M>
riesz_tag2 (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const quadrature_option_type& qopt,
    const Geo&                    dom,
    const T&)
{
  return riesz_function_internal (Xh, f, qopt, dom, std::false_type());
}
template <class T, class M, class Geo, class Function>
field_basic<T,M>
riesz_tag2 (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const quadrature_option_type& qopt,
    const Geo&                    dom,
    const point_basic<T>&)
{
  return riesz_function_internal (Xh, f, qopt, dom, std::true_type());
}
template <class T, class M, class Geo, class Function>
field_basic<T,M>
riesz_tag2 (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const quadrature_option_type& qopt,
    const Geo&                    dom,
    const undeterminated_basic<T>&)
{
  if (Xh.valued_tag() == space_constant::scalar) {
    return riesz_function_internal (Xh, f, qopt, dom, std::false_type());
  } else {
    return riesz_function_internal (Xh, f, qopt, dom, std::true_type());
  }
}
// --------------------------------------------------------------------------
// f: switch between constant or function
// --------------------------------------------------------------------------
// f is a function:
template <class T, class M, class Geo, class Function>
typename
std::enable_if<
  is_field_function<Function>::value
 ,field_basic<T,M>
>::type
riesz_tag (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const quadrature_option_type& qopt,
    const Geo&                    dom,
    std::false_type)
{
  using result_type = typename field_function_traits<Function>::result_type;
  return riesz_tag2 (Xh, f, qopt, dom, result_type());
}
// f is an expression
template <class T, class M, class Geo, class Function>
typename
std::enable_if<
  ! is_field_function<Function>::value
 ,field_basic<T,M>
>::type
riesz_tag (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const quadrature_option_type& qopt,
    const Geo&                    dom,
    std::false_type)
{
  using result_type = typename Function::result_type;
  return riesz_tag2 (Xh, f, qopt, dom, result_type());
}
// f is a constant:
// thanks to a wrapper class, we go to the first case
template <class T, class M, class Geo, class Constant>
field_basic<T,M>
riesz_tag (
    const space_basic<T,M>&       Xh,
    const Constant&               value,
    const quadrature_option_type& qopt,
    const Geo&                    dom,
    std::true_type)
{
  typedef typename field_basic<T,M>::float_type float_type;
  typedef typename constant_promote<Constant>::type constant_type; // int to Float, but conserve point_basic<T> & tensor
  return riesz_tag (Xh, f_constant<point_basic<float_type>,constant_type>(value),
			     qopt, dom, std::false_type());
}

} // namespace details

/*Class:riesz
NAME: @code{riesz} - approximate a Riesz representer
@findex  riesz
@cindex  riesz representer
@cindex  quadrature formula
@clindex space
@clindex field

SYNOPSYS:
  The @code{riesz} function is now obsolete: it has been now suppersetted by the @code{integrate} function
  @pxref{integrate algorithm}.
@example
 template <class Expr>
 field riesz (space, Expr expr);
 field riesz (space, Expr expr, quadrature_option_type);
 field riesz (space, Expr expr, domain);
 field riesz (space, Expr expr, domain, quadrature_option_type);
@end example
  The domain can be also provided by its name as a string.
  The old-fashioned code:
NOTE:
  The @code{riesz} function is now obsolete: it has been now suppersetted by the @code{integrate} function
  @pxref{integrate algorithm}.
  The old-fashioned code:
  @example
     field l1h = riesz (Xh, f);
     field l2h = riesz (Xh, f, "boundary");
  @end example
  writes now:
  @example
     test v (Xh);
     field l1h = integrate (f*v);
     field l2h = integrate ("boundary", f*v);
  @end example
  The @code{riesz} function is still present in the library for backward compatibility purpose.
DESCRIPTION:
  Let @code{f} be any continuous function, its Riesz representer in the finite
  element space @code{Xh} on the domain @code{Omega} is defind by:
@example
                /
                |
  dual(lh,vh) = |      f(x) vh(x) dx
                |
               / Omega
@end example
  for all @code{vh} in @code{Xh}, where @code{dual} denotes the duality
  between @code{Xh} and its dual. As @code{Xh} is a finite dimensional space,
  its dual is identified as @code{Xh} and the duality product as the Euclidian one.
  The Riesz representer is thus the @code{lh} field of @code{Xh} where 
  its i-th degree of freedom is:
@example
                /
                |
  dual(lh,vh) = |      f(x) phi_i(x) dx
                |
               / Omega
@end example
  where phi_i is the i-th basis function in @code{Xh}.
  The integral is evaluated by using a quadrature formula.
  By default the quadrature formule is the Gauss one with
  the order equal to @code{2*k-1} where $@code{k} is the polynomial degree in @code{Xh}.
  Alternative quadrature formula and order is available
  by passing an optional variable to riesz.

@noindent
  The function @code{riesz} implements the
  approximation of the Riesz representer 
  by using some quadrature formula for the evaluation of the integrals.
  Its argument can be any function, class-function or linear or nonlinear
  expressions mixing fields and continuous functions.

@cindex  integrate
@findex  integrate
EXAMPLE:
@noindent
  The following code compute the Riesz representant, denoted
  by @code{lh} of f(x), and the integral of f over the domain omega:
@example
  Float f(const point& x);
  ...
  space Xh (omega_h, "P1");
  field lh = riesz (Xh, f);
  Float int_f = dual(lh, 1);
@end example
OPTIONS:
  An optional argument specifies the quadrature formula used
  for the computation of the integral.
  The domain of integration is by default the mesh associated to
  the finite element space.
  An alternative domain @code{dom}, e.g. a part of 
  the boundary can be supplied as an extra argument.
  This domain can be also a @code{band} associated to
  the banded level set method.
End: */

//<riesz:
template <class T, class M, class Function>
inline
field_basic<T,M>
riesz (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const quadrature_option_type& qopt
       = quadrature_option_type())
//>riesz:
{
  return details::riesz_tag (Xh, f, qopt, Xh.get_geo(),
	typename details::is_constant<Function>::type());
}
//<riesz:
template <class T, class M, class Function>
field_basic<T,M>
riesz (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const geo_basic<T,M>&         dom,
    const quadrature_option_type& qopt
       = quadrature_option_type())
//>riesz:
{
  return details::riesz_tag (Xh, f, qopt, dom,
	typename details::is_constant<Function>::type());
}
//<riesz:
template <class T, class M, class Function>
field_basic<T,M>
riesz (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    std::string                   dom_name,
    const quadrature_option_type& qopt
       = quadrature_option_type())
//>riesz:
{
  const geo_basic<T,M>& omega = Xh.get_geo();
  return riesz (Xh, f, omega[dom_name], qopt);
}
//<riesz:
template <class T, class M, class Function>
field_basic<T,M>
riesz (
    const space_basic<T,M>&       Xh,
    const Function&               f,
    const band_basic<T,M>&        gh,
    const quadrature_option_type& qopt
       = quadrature_option_type())
//>riesz:
{
  return details::riesz_tag (Xh, f, qopt, gh,
	typename details::is_constant<Function>::type());
}

}// namespace rheolef
#endif // _RHEO_RIESZ_H