librheolef-dev 5.93-2 / 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
/// 
/// =========================================================================

/*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 internal}).
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:
*/

#include "rheolef/reference_element.h"
#include "rheolef/point.h"

namespace rheolef { 
struct weighted_point {
    weighted_point ();
    weighted_point (const point& x, Float w);
    point x;
    Float w;
};

class quadrature_option_type {
public:
    typedef enum {
	gauss		= 0,
	gauss_lobatto	= 1,
	gauss_radau	= 2,
	middle_edge	= 3,
	superconvergent	= 4,
	max_family	= 5
    } family_type; // update also family_name[] in quatrature.cc

    typedef size_t size_type;
    quadrature_option_type(
        family_type ft = quadrature_option_type::gauss,
        size_type k = std::numeric_limits<size_type>::max());
    size_t         get_order() const;
    family_type    get_family() const;
    std::string    get_family_name() const;
    void set_order (size_t r);
    void set_family (family_type ft);
protected:
    family_type   _family;
    size_t        _order;
};
class quadrature_on_geo : public std::vector<weighted_point> {
public:

// typedefs:

    typedef std::vector<weighted_point>:: size_type size_type;
    typedef point                                   x;

// alocators/deallocators:

    quadrature_on_geo ();

// 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& x, const Float& w) {
	    push_back (weighted_point(x,w)); }

    static size_type n_node_gauss (size_t r);

    friend std::ostream& operator<< (std::ostream&, const quadrature_on_geo&);
private:
    quadrature_on_geo (const quadrature_on_geo&);
    quadrature_on_geo operator= (const quadrature_on_geo&);
};
//<quadrature:
class quadrature {
public:

// typedefs:

    typedef quadrature_on_geo::size_type size_type;
    typedef quadrature_option_type::family_type family_type;
    typedef std::vector<weighted_point>::const_iterator const_iterator;

// allocators:

    quadrature (quadrature_option_type opt = quadrature_option_type());

// 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;
    size_type      size  (reference_element hat_K) const;
    const_iterator begin (reference_element hat_K) const;
    const_iterator end   (reference_element hat_K) const;
    friend std::ostream& operator<< (std::ostream&, const quadrature&);
protected:
    quadrature_option_type     _options;
    mutable quadrature_on_geo  _quad        [reference_element::max_size];
    mutable std::vector<bool>  _initialized;
    void _initialize (reference_element hat_K) const;
private:
    quadrature (const quadrature&);
    quadrature operator= (const quadrature&);
};
//>quadrature:
// ------------------------------------------------------------
// inlined
// ------------------------------------------------------------
inline 
quadrature_option_type::quadrature_option_type (family_type ft, size_type k)
 : _family(ft),
   _order(k)
{
}
inline
quadrature_option_type::size_type
quadrature_option_type::get_order () const
{
    return _order;
}
inline
quadrature_option_type::family_type
quadrature_option_type::get_family () const
{
    return _family;
}
inline
void
quadrature_option_type::set_order (size_t r)
{
    _order = r;
}
inline
void
quadrature_option_type::set_family (family_type ft)
{
    _family = ft;
}
inline
weighted_point::weighted_point ()
 : x(), w(0) {}

inline
weighted_point::weighted_point (const point& x1, Float w1)
 : x(x1), w(w1) {}

inline 
quadrature_on_geo::quadrature_on_geo ()
 : std::vector<weighted_point>()
{
}
inline 
quadrature::size_type
quadrature::get_order () const
{
    return _options.get_order();
}
inline 
quadrature::family_type
quadrature::get_family () const
{
    return _options.get_family();
}
inline 
std::string
quadrature::get_family_name () const
{
    return _options.get_family_name();
}
inline 
quadrature::quadrature (quadrature_option_type opt)
 : _options(opt),
   _quad(),
   _initialized(reference_element::max_size, false)
{
}

inline
void
quadrature::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);
  }
}
inline
void
quadrature::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);
  }
}
inline
quadrature_on_geo::size_type
quadrature_on_geo::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);
}
}// namespace rheolef
#endif // _RHEO_QUADRATURE_H