librheolef-dev 6.5-1build1 / usr / include / rheolef / quadrature.h

#ifndef _RHEO_QUADRATURE_H
#define _RHEO_QUADRATURE_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/smart_pointer.h"
#include "rheolef/reference_element.h"
#include "rheolef/point.h"
#include "rheolef/quadrature_option_type.h"
#include "rheolef/reference_element_face_transformation.h"

namespace rheolef { 

/*Class:quadrature
NAME: @code{quadrature} - quadrature formulae on the reference lement
@cindex  quadrature formulae
@clindex quadrature
@clindex reference_element
SYNOPSYS:
@noindent
The @code{quadrature} class defines a container for a quadrature
formulae on the reference element (@pxref{reference_element iclass}).
This container stores the nodes coordinates and the weights.
The constructor takes two arguments: the reference element 
@math{K}
and the order @math{r} of the quadrature formulae.
The formulae is exact when computing the integral
of a polynom @math{p} that degree is less or equal to order @math{r}.
@example
                  n
    /            ___
    | p(x) dx =  \    p(x_q) w_q
    / K          /__
                 q=1
@end example
LIMITATIONS:
@noindent
The formulae is optimal when it uses a minimal number
of nodes @math{n}.
Optimal quadrature formula are hard-coded in this class.
Not all reference elements and orders are yet 
implemented. This class will be completed in the future.

AUTHORS:
    LMC-IMAG, 38041 Grenoble cedex 9, France
   | Pierre.Saramito@imag.fr
DATE:   30 november 2003
End:
*/

template<class T>
struct weighted_point {
  weighted_point () : x(), w(0) {}
  weighted_point (const point_basic<T>& x1, const T& w1) : x(x1), w(w1) {}
  point_basic<T> x;
  T              w;
};

template<class T>
class quadrature_on_geo : public std::vector<weighted_point<T> > {
public:

// typedefs:

    typedef std::vector<weighted_point<T> > base;
    typedef typename base::size_type        size_type;
    typedef point_basic<T>                  x;

// alocators/deallocators:

    quadrature_on_geo ();

    quadrature_on_geo (const quadrature_on_geo& q) : base(q) {}
    quadrature_on_geo& operator= (const quadrature_on_geo& q) {
	base::operator= (q);
        return *this;
    }

// modifier:

    void initialize       (reference_element hat_K, quadrature_option_type opt);
    void init_point       (quadrature_option_type opt);
    void init_edge        (quadrature_option_type opt);
    void init_triangle    (quadrature_option_type opt);
    void init_square      (quadrature_option_type opt);
    void init_tetrahedron (quadrature_option_type opt);
    void init_prism       (quadrature_option_type opt);
    void init_hexahedron  (quadrature_option_type opt);

    void wx (const point_basic<T>& x, const T& w) {
	    base::push_back (weighted_point<T>(x,w)); }

    static size_type n_node_gauss (size_type r);
    
    template<class U>
    friend std::ostream& operator<< (std::ostream&, const quadrature_on_geo<U>&);
};
template<class T>
class quadrature_rep {
public:

// typedefs:

    typedef typename quadrature_on_geo<T>::size_type                 size_type;
    typedef quadrature_option_type::family_type                      family_type;
    typedef typename std::vector<weighted_point<T> >::const_iterator const_iterator;
    typedef geo_element_indirect::orientation_type 		     orientation_type;

// allocators:

    quadrature_rep (quadrature_option_type opt = quadrature_option_type());
    quadrature_rep (const quadrature_rep<T>& q);
    const quadrature_rep& operator= (const quadrature_rep<T>& q);

// modifiers:

    void set_order  (size_type order);
    void set_family (family_type ft);

// accessors:

    size_type      get_order() const;
    family_type    get_family() const;
    std::string    get_family_name() const;
    const quadrature_option_type& get_options() const;
    size_type      size  (reference_element hat_K) const;
    const_iterator begin (reference_element hat_K) const;
    const_iterator end   (reference_element hat_K) const;
    const weighted_point<T>& operator() (reference_element hat_K, size_type q) const;

    template<class U>
    friend std::ostream& operator<< (std::ostream&, const quadrature_rep<U>&);

#ifdef TO_CLEAN
// side accessor:
    void side_initialize (
        reference_element             tilde_K,
        const side_information_type&  sid) const;
#endif // TO_CLEAN

protected:
// internal:
    void _initialize (reference_element hat_K) const;
// data:
    quadrature_option_type        _options;
    mutable boost::array<quadrature_on_geo<T>,reference_element::max_variant>
				  _quad;
    mutable std::vector<bool>     _initialized;
private:
};
//<quadrature:
template<class T>
class quadrature : public smart_pointer<quadrature_rep<T> > {
public:

// typedefs:

    typedef quadrature_rep<T>              rep;
    typedef smart_pointer<rep>             base;
    typedef typename rep::size_type        size_type;
    typedef typename rep::family_type      family_type;
    typedef typename rep::const_iterator   const_iterator;
    typedef typename rep::orientation_type orientation_type;

// allocators:

    quadrature (quadrature_option_type opt = quadrature_option_type())
     : base(new_macro(rep(opt))) {}

// modifiers:

    void set_order  (size_type order) { base::data().set_order(order); }
    void set_family (family_type ft)  { base::data().set_family(ft); }

// accessors:

    size_type      get_order() const { return base::data().get_order();}
    family_type    get_family() const { return base::data().get_family();}
    std::string    get_family_name() const { return base::data().get_family_name();}
    const quadrature_option_type& get_options() const { return base::data().get_options(); }
    size_type      size  (reference_element hat_K) const { return base::data().size(hat_K); }
    const_iterator begin (reference_element hat_K) const { return base::data().begin(hat_K); }
    const_iterator end   (reference_element hat_K) const { return base::data().end(hat_K); }
    const weighted_point<T>& operator() (reference_element hat_K, size_type q) const
    							 { return base::data().operator() (hat_K,q); }
    template<class U>
    friend std::ostream& operator<< (std::ostream& os, const quadrature<U>& q) {
	return os << q.data(); }

// side accessor:

    void side_initialize (
        reference_element             tilde_K,
        size_type                     loc_isid,
        reference_element             hat_S,
        size_type                     shift,
        orientation_type              orient) const {
    	    base::data().side_initialize (tilde_K, loc_isid, hat_S, shift, orient);
    }
};
//>quadrature:
// ------------------------------------------------------------
// inlined
// ------------------------------------------------------------
template<class T>
inline 
quadrature_on_geo<T>::quadrature_on_geo ()
 : std::vector<weighted_point<T> >()
{
}
template<class T>
inline
typename quadrature_on_geo<T>::size_type
quadrature_on_geo<T>::n_node_gauss (size_type r)
{
    // when using n nodes : gauss quadrature formulae order is r=2*n-1
    if (r == 0) return 1;
    size_type n = (r % 2 == 0) ? r/2+1 : (r+1)/2;
    return std::max(size_t(1), n);
}
template<class T>
inline 
typename quadrature_rep<T>::size_type
quadrature_rep<T>::get_order () const
{
    return _options.get_order();
}
template<class T>
inline 
typename quadrature_rep<T>::family_type
quadrature_rep<T>::get_family () const
{
    return _options.get_family();
}
template<class T>
inline 
std::string
quadrature_rep<T>::get_family_name () const
{
    return _options.get_family_name();
}
template<class T>
inline 
const quadrature_option_type&
quadrature_rep<T>::get_options () const
{
    return _options;
}
template<class T>
inline
void
quadrature_rep<T>::set_order (size_type r)
{
  if (get_order() != r) {
    // do not re-initialize nodes-weights if unchanged
    _options.set_order(r);
    std::fill (_initialized.begin(), _initialized.end(), false);
  }
}
template<class T>
inline
void
quadrature_rep<T>::set_family (family_type ft)
{
  if (get_family() != ft) {
    // do not re-initialize nodes-weights if unchanged
    _options.set_family(ft);
    std::fill (_initialized.begin(), _initialized.end(), false);
  }
}

}// namespace rheolef
#endif // _RHEO_QUADRATURE_H