/usr/include/rheolef/riesz

#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/characteristic.h"
#include "rheolef/field_element.h"
namespace rheolef { 

/*Class:riesz_representer
NAME: @code{riesz_representer} - integrate a function by using quadrature formulae
@findex  riesz_representer
@cindex  riesz representer
@cindex  quadrature formulae
@clindex space
@clindex field

DESCRIPTION:
 The function @code{riesz_representer} implements the
 approximation of an integral by using quadrature formulae.
 This is an expreimental implementation: please do not use 
 yet for practical usage.
SYNOPSYS:
 template <class Function>
 field riesz_representer (const space& Vh, const Function& f);
    
 template <class Function>
 field riesz_representer (const space& Vh, const Function& f,
    	quadrature_option_type qopt);
EXAMPLE:
@noindent
 The following code compute the Riesz representant, denoted
 by @code{mfh} of f(x), and the integral of f over the domain omega:
@example
  Float f(const point& x);
  ...
  space Vh (omega_h, "P1");
  field mfh = riesz_representer(Vh, f);
  Float int_f = dot(mfh, field(Vh,1.0));
@end example
 The Riesz representer is the @code{mfh} vector of values:
@example
	mfh(i) = integrate f(x) phi_i(x) dx
@end example
 where phi_i is the i-th basis function in Vh
 and the integral is evaluated by using a quadrature formulae.
 By default the quadrature formule is the Gauss one with
 the order equal to the polynomial order of Vh.
 Alternative quadrature formulae and order is available
 by passing an optional variable to riesz_representer.
End: */

template <class ElementIterator, class Function>
inline
void
assembly (
    ElementIterator               iter_element,
    ElementIterator               last_element,
    const field_element&          field_e,
    const Function&               fct,
    field&                        mfh)
{
    typedef size_t size_type;
    const space& X = mfh.get_space();
    std::vector<Float> uk;
    tiny_vector<size_type> idx;

    for (; iter_element != last_element; iter_element++) {
      const geo_element& K = *iter_element;
      // --------------------------------
      // compute elementary matrix ak
      // --------------------------------
      field_e (K, fct, uk);
      // ------------------------------
      // assembly ak in sparse matrix a
      // ------------------------------
      X.set_global_dof (K, idx);
      size_type nx = uk.size();
      for (size_type i = 0; i < nx; i++) {
	    size_type i_dof = X.domain_dof(idx[i]);
	    size_type ii    = X.index(i_dof);
	    if (X.is_blocked(i_dof)) mfh.b.at(ii) += uk[i];
	    else                     mfh.u.at(ii) += uk[i];
      }
    }
}
template <class Function>
inline
void
field_assembly (
   field&                        mfh, 
   const field_element&          field_e,
   const Function&               fct)
{
  const space& X = mfh.get_space();
  const geo&   g = mfh.get_geo();
  assembly (g.begin(), g.end(), field_e, fct, mfh);
}
//<riesz_representer:
template <class Function>
field
riesz_representer (
    const space& Vh,
    const Function& f,
    quadrature_option_type qopt
       = quadrature_option_type(quadrature_option_type::max_family,0))
//>riesz_representer:
{
  field_element field_e;
  if (qopt.get_family() == quadrature_option_type::max_family) {
    qopt.set_family(quadrature_option_type::gauss_lobatto);
    qopt.set_order (Vh.degree());
  }
  field_e.initialize (Vh, qopt);
  field mfh (Vh, 0.);
  field_assembly (mfh, field_e, f);
  return mfh;
}
// --------------------------------------------------------------------
// integrate:
//   Float value = integrate (omega, f, quadrature_option);
// where f is a template function
// --------------------------------------------------------------------
template <class Function>
inline
Float
integrate (const geo& g, Function f, quadrature_option_type qopt) {
    space Xh (g, "P0");
    field int_K_f_dx = riesz_representer (Xh, f, qopt);
    // TODO: optimize with a sum !
    return dot (field(Xh,1.0), int_K_f_dx);
}
}// namespace rheolef
#endif // _RHEO_RIESZ_H
librheolef-dev 5.93-2 / usr / include / rheolef / riesz_representer.h
