Loading...
Searching...
No Matches
quadrature.h
Go to the documentation of this file.
#ifndef _RHEO_QUADRATURE_H
#define _RHEO_QUADRATURE_H
#include "rheolef/smart_pointer.h"
#include "rheolef/persistent_table.h"
#include "rheolef/reference_element.h"
#include "rheolef/point.h"
#include "rheolef/integrate_option.h"
#include "rheolef/reference_element_face_transformation.h"
#include "rheolef/compiler_eigen.h"
namespace rheolef {
/*Class:quadrature
NAME: @code{quadrature} - quadrature formulae on the reference lement
@cindex quadrature formulae
@clindex quadrature
@clindex reference_element
SYNOPSIS:
@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;
};
// ----------------------------------------------------------------------------
// quadrature on a specific ref element
// ----------------------------------------------------------------------------
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;
}
// accessor:
void get_nodes (Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const;
// modifier:
void initialize (reference_element hat_K, quadrature_option opt);
void init_point (quadrature_option opt);
void init_edge (quadrature_option opt);
void init_triangle (quadrature_option opt);
void init_square (quadrature_option opt);
void init_tetrahedron (quadrature_option opt);
void init_prism (quadrature_option opt);
void init_hexahedron (quadrature_option 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>&);
};
// ----------------------------------------------------------------------------
// quadrature representation
// ----------------------------------------------------------------------------
template<class T>
class quadrature_rep {
public:
// typedefs:
typedef typename quadrature_on_geo<T>::size_type size_type;
typedef quadrature_option::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_rep (quadrature_option opt = quadrature_option());
quadrature_rep (const std::string& name);
quadrature_rep (const quadrature_rep<T>& q);
const quadrature_rep& operator= (const quadrature_rep<T>& q);
#ifdef TO_CLEAN
// modifiers:
void set_order (size_type order);
void set_family (family_type ft);
#endif // TO_CLEAN
// accessors:
std::string name() const;
size_type get_order() const;
family_type get_family() const;
std::string get_family_name() const;
const quadrature_option& 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;
void get_nodes (reference_element hat_K, Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const;
template<class U>
friend std::ostream& operator<< (std::ostream&, const quadrature_rep<U>&);
protected:
// internal:
void _initialize (reference_element hat_K) const;
// data:
quadrature_option _options;
mutable std::array<quadrature_on_geo<T>,reference_element::max_variant>
_quad;
mutable std::vector<bool> _initialized;
public:
static quadrature_rep* make_ptr (const std::string& name) { return new_macro(quadrature_rep(name)); }
};
// ----------------------------------------------------------------------------
// quadrature class
// ----------------------------------------------------------------------------
//<quadrature:
template<class T>
class quadrature : public smart_pointer<quadrature_rep<T> >,
public persistent_table<quadrature<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 (const std::string& name = "");
quadrature (quadrature_option opt);
// modifiers:
void set_order (size_type order);
void set_family (family_type ft);
void reset (const std::string& name);
// accessors:
const quadrature_option& get_options() const;
std::string name() const { return get_options().name(); }
size_type get_order() const { return get_options().get_order();}
family_type get_family() const { return get_options().get_family();}
std::string get_family_name() const { return get_options().get_family_name();}
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); }
void get_nodes (reference_element hat_K, Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const
{ return base::data().get_nodes(hat_K,node); }
};
//>quadrature:
template<class T>
inline
std::ostream& operator<< (std::ostream& os, const quadrature<T>& q)
{
return os << q.data();
}
// ------------------------------------------------------------
// 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&
quadrature_rep<T>::get_options () const
{
return _options;
}
#ifdef TO_CLEAN
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);
}
}
#endif // TO_CLEAN
}// namespace rheolef
#endif // _RHEO_QUADRATURE_H
size_type
field::size_type size_type
Definition branch.cc:430
rheolef::geo_element_indirect::orientation_type
short int orientation_type
Definition geo_element_indirect.h:38
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition integrate_option.h:125
rheolef::integrate_option::set_order
void set_order(size_t r)
Definition integrate_option.h:254
rheolef::persistent_table
see the persistent_table page for the full documentation
Definition persistent_table.h:84
rheolef::point_basic
Definition point.h:86
rheolef::quadrature_on_geo
Definition quadrature.h:82
rheolef::quadrature_on_geo::init_tetrahedron
void init_tetrahedron(quadrature_option opt)
Definition quadrature_Te.cc:28
rheolef::quadrature_on_geo::x
point_basic< T > x
Definition quadrature.h:89
rheolef::quadrature_on_geo::initialize
void initialize(reference_element hat_K, quadrature_option opt)
Definition quadrature_rep.cc:105
rheolef::quadrature_on_geo::operator<<
friend std::ostream & operator<<(std::ostream &, const quadrature_on_geo< U > &)
rheolef::quadrature_on_geo::quadrature_on_geo
quadrature_on_geo()
Definition quadrature.h:238
rheolef::quadrature_on_geo::init_prism
void init_prism(quadrature_option opt)
Definition quadrature_Pr.cc:28
rheolef::quadrature_on_geo::init_point
void init_point(quadrature_option opt)
Definition quadrature_p.cc:27
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition quadrature.h:88
rheolef::quadrature_on_geo::n_node_gauss
static size_type n_node_gauss(size_type r)
Definition quadrature.h:245
rheolef::quadrature_on_geo::get_nodes
void get_nodes(Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
Definition quadrature_rep.cc:123
rheolef::quadrature_on_geo::operator=
quadrature_on_geo & operator=(const quadrature_on_geo &q)
Definition quadrature.h:96
rheolef::quadrature_on_geo::init_hexahedron
void init_hexahedron(quadrature_option opt)
Definition quadrature_H.cc:28
rheolef::quadrature_on_geo::init_square
void init_square(quadrature_option opt)
Definition quadrature_q.cc:46
rheolef::quadrature_on_geo::init_edge
void init_edge(quadrature_option opt)
Definition quadrature_e.cc:38
rheolef::quadrature_on_geo::quadrature_on_geo
quadrature_on_geo(const quadrature_on_geo &q)
Definition quadrature.h:95
rheolef::quadrature_on_geo::init_triangle
void init_triangle(quadrature_option opt)
Definition quadrature_t.cc:50
rheolef::quadrature_on_geo::base
std::vector< weighted_point< T > > base
Definition quadrature.h:87
rheolef::quadrature_on_geo::wx
void wx(const point_basic< T > &x, const T &w)
Definition quadrature.h:116
rheolef::quadrature_rep
Definition quadrature.h:128
rheolef::quadrature_rep::operator<<
friend std::ostream & operator<<(std::ostream &, const quadrature_rep< U > &)
rheolef::quadrature_rep::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition quadrature.h:136
rheolef::quadrature_rep::const_iterator
std::vector< weighted_point< T > >::const_iterator const_iterator
Definition quadrature.h:135
rheolef::quadrature_rep::~quadrature_rep
~quadrature_rep()
Definition quadrature_rep.cc:246
rheolef::quadrature_rep::operator=
const quadrature_rep & operator=(const quadrature_rep< T > &q)
Definition quadrature_rep.cc:157
rheolef::quadrature_rep::begin
const_iterator begin(reference_element hat_K) const
Definition quadrature_rep.cc:186
rheolef::quadrature_rep::size_type
quadrature_on_geo< T >::size_type size_type
Definition quadrature.h:133
rheolef::quadrature_rep::_quad
std::array< quadrature_on_geo< T >, reference_element::max_variant > _quad
Definition quadrature.h:175
rheolef::quadrature_rep::end
const_iterator end(reference_element hat_K) const
Definition quadrature_rep.cc:194
rheolef::quadrature_rep::get_nodes
void get_nodes(reference_element hat_K, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
Definition quadrature_rep.cc:179
rheolef::quadrature_rep::set_family
void set_family(family_type ft)
Definition quadrature.h:295
rheolef::quadrature_rep::family_type
quadrature_option::family_type family_type
Definition quadrature.h:134
rheolef::quadrature_rep::set_order
void set_order(size_type order)
Definition quadrature.h:284
rheolef::quadrature_rep::get_family_name
std::string get_family_name() const
Definition quadrature.h:269
rheolef::quadrature_rep::get_order
size_type get_order() const
Definition quadrature.h:255
rheolef::quadrature_rep::_initialized
std::vector< bool > _initialized
Definition quadrature.h:176
rheolef::quadrature_rep::size
size_type size(reference_element hat_K) const
Definition quadrature_rep.cc:202
rheolef::quadrature_rep::operator()
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
Definition quadrature_rep.cc:210
rheolef::quadrature_rep::_options
quadrature_option _options
Definition quadrature.h:173
rheolef::quadrature_rep::_initialize
void _initialize(reference_element hat_K) const
Definition quadrature_rep.cc:171
rheolef::quadrature_rep::make_ptr
static quadrature_rep * make_ptr(const std::string &name)
Definition quadrature.h:178
rheolef::quadrature_rep::get_options
const quadrature_option & get_options() const
Definition quadrature.h:276
rheolef::quadrature_rep::name
std::string name() const
Definition quadrature_rep.cc:165
rheolef::quadrature_rep::get_family
family_type get_family() const
Definition quadrature.h:262
rheolef::quadrature
Definition quadrature.h:186
rheolef::quadrature::const_iterator
rep::const_iterator const_iterator
Definition quadrature.h:195
rheolef::quadrature::rep
quadrature_rep< T > rep
Definition quadrature.h:191
rheolef::quadrature::get_family
family_type get_family() const
Definition quadrature.h:215
rheolef::quadrature::name
std::string name() const
Definition quadrature.h:213
rheolef::quadrature::end
const_iterator end(reference_element hat_K) const
Definition quadrature.h:219
rheolef::quadrature::reset
void reset(const std::string &name)
Definition quadrature_rep.cc:254
rheolef::quadrature::get_nodes
void get_nodes(reference_element hat_K, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
Definition quadrature.h:222
rheolef::quadrature::set_family
void set_family(family_type ft)
Definition quadrature_rep.cc:292
rheolef::quadrature::orientation_type
rep::orientation_type orientation_type
Definition quadrature.h:196
rheolef::quadrature::set_order
void set_order(size_type order)
Definition quadrature_rep.cc:284
rheolef::quadrature::family_type
rep::family_type family_type
Definition quadrature.h:194
rheolef::quadrature::get_family_name
std::string get_family_name() const
Definition quadrature.h:216
rheolef::quadrature::base
smart_pointer< rep > base
Definition quadrature.h:192
rheolef::quadrature::operator()
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
Definition quadrature.h:220
rheolef::quadrature::size
size_type size(reference_element hat_K) const
Definition quadrature.h:217
rheolef::quadrature::size_type
rep::size_type size_type
Definition quadrature.h:193
rheolef::quadrature::get_order
size_type get_order() const
Definition quadrature.h:214
rheolef::quadrature::get_options
const quadrature_option & get_options() const
Definition quadrature_rep.cc:264
rheolef::quadrature::begin
const_iterator begin(reference_element hat_K) const
Definition quadrature.h:218
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition reference_element.h:82
rheolef::smart_pointer_base< T, details::constructor_copy< T > >::data
const T & data() const
Definition smart_pointer.h:266
rheolef::smart_pointer
see the smart_pointer page for the full documentation
Definition smart_pointer.h:351
T
Expr1::float_type T
Definition field_expr.h:230
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition catchmark.h:99
std
STL namespace.
rheolef::weighted_point
Definition quadrature.h:71
rheolef::weighted_point::w
T w
Definition quadrature.h:75
rheolef::weighted_point::weighted_point
weighted_point(const point_basic< T > &x1, const T &w1)
Definition quadrature.h:73
rheolef::weighted_point::weighted_point
weighted_point()
Definition quadrature.h:72
rheolef::weighted_point::x
point_basic< T > x
Definition quadrature.h:74