Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature.h
Go to the documentation of this file.
1#ifndef _RHEO_QUADRATURE_H
2#define _RHEO_QUADRATURE_H
23#include "rheolef/smart_pointer.h"
24#include "rheolef/persistent_table.h"
25#include "rheolef/reference_element.h"
26#include "rheolef/point.h"
27#include "rheolef/integrate_option.h"
28#include "rheolef/reference_element_face_transformation.h"
29#include "rheolef/compiler_eigen.h"
30
31namespace rheolef {
32
33/*Class:quadrature
34NAME: @code{quadrature} - quadrature formulae on the reference lement
35@cindex quadrature formulae
36@clindex quadrature
37@clindex reference_element
38SYNOPSIS:
39@noindent
40The @code{quadrature} class defines a container for a quadrature
41formulae on the reference element (@pxref{reference_element iclass}).
42This container stores the nodes coordinates and the weights.
43The constructor takes two arguments: the reference element
44@math{K}
45and the order @math{r} of the quadrature formulae.
46The formulae is exact when computing the integral
47of a polynom @math{p} that degree is less or equal to order @math{r}.
48@example
49 n
50 / ___
51 | p(x) dx = \ p(x_q) w_q
52 / K /__
53 q=1
54@end example
55LIMITATIONS:
56@noindent
57The formulae is optimal when it uses a minimal number
58of nodes @math{n}.
59Optimal quadrature formula are hard-coded in this class.
60Not all reference elements and orders are yet
61implemented. This class will be completed in the future.
62
63AUTHORS:
64 LMC-IMAG, 38041 Grenoble cedex 9, France
65 | Pierre.Saramito@imag.fr
66DATE: 30 november 2003
67End:
68*/
69
70template<class T>
72 weighted_point () : x(), w(0) {}
73 weighted_point (const point_basic<T>& x1, const T& w1) : x(x1), w(w1) {}
76};
77
78// ----------------------------------------------------------------------------
79// quadrature on a specific ref element
80// ----------------------------------------------------------------------------
81template<class T>
82class quadrature_on_geo : public std::vector<weighted_point<T> > {
83public:
84
85// typedefs:
86
87 typedef std::vector<weighted_point<T> > base;
88 typedef typename base::size_type size_type;
90
91// alocators/deallocators:
92
94
97 base::operator= (q);
98 return *this;
99 }
100
101// accessor:
102
103 void get_nodes (Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const;
104
105// modifier:
106
108 void init_point (quadrature_option opt);
109 void init_edge (quadrature_option opt);
113 void init_prism (quadrature_option opt);
115
116 void wx (const point_basic<T>& x, const T& w) {
117 base::push_back (weighted_point<T>(x,w)); }
118
120
121 template<class U>
122 friend std::ostream& operator<< (std::ostream&, const quadrature_on_geo<U>&);
123};
124// ----------------------------------------------------------------------------
125// quadrature representation
126// ----------------------------------------------------------------------------
127template<class T>
129public:
130
131// typedefs:
132
134 typedef quadrature_option::family_type family_type;
135 typedef typename std::vector<weighted_point<T> >::const_iterator const_iterator;
137
138// allocators:
139
141 quadrature_rep (quadrature_option opt = quadrature_option());
142 quadrature_rep (const std::string& name);
145
146#ifdef TO_CLEAN
147// modifiers:
148
149 void set_order (size_type order);
150 void set_family (family_type ft);
151#endif // TO_CLEAN
152
153// accessors:
154
155 std::string name() const;
156 size_type get_order() const;
157 family_type get_family() const;
158 std::string get_family_name() const;
159 const quadrature_option& get_options() const;
160 size_type size (reference_element hat_K) const;
164 void get_nodes (reference_element hat_K, Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const;
165
166 template<class U>
167 friend std::ostream& operator<< (std::ostream&, const quadrature_rep<U>&);
168
169protected:
170// internal:
171 void _initialize (reference_element hat_K) const;
172// data:
173 quadrature_option _options;
174 mutable std::array<quadrature_on_geo<T>,reference_element::max_variant>
176 mutable std::vector<bool> _initialized;
177public:
178 static quadrature_rep* make_ptr (const std::string& name) { return new_macro(quadrature_rep(name)); }
179};
180// ----------------------------------------------------------------------------
181// quadrature class
182// ----------------------------------------------------------------------------
183//<quadrature:
184template<class T>
185class quadrature : public smart_pointer<quadrature_rep<T> >,
186 public persistent_table<quadrature<T>> {
187public:
188
189// typedefs:
190
193 typedef typename rep::size_type size_type;
197
198// allocators:
199
200
201 quadrature (const std::string& name = "");
202 quadrature (quadrature_option opt);
203
204// modifiers:
205
206 void set_order (size_type order);
207 void set_family (family_type ft);
208 void reset (const std::string& name);
209
210// accessors:
211
212 const quadrature_option& get_options() const;
213 std::string name() const { return get_options().name(); }
214 size_type get_order() const { return get_options().get_order();}
215 family_type get_family() const { return get_options().get_family();}
216 std::string get_family_name() const { return get_options().get_family_name();}
217 size_type size (reference_element hat_K) const { return base::data().size(hat_K); }
218 const_iterator begin (reference_element hat_K) const { return base::data().begin(hat_K); }
219 const_iterator end (reference_element hat_K) const { return base::data().end(hat_K); }
221 { return base::data().operator() (hat_K,q); }
222 void get_nodes (reference_element hat_K, Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const
223 { return base::data().get_nodes(hat_K,node); }
224};
225//>quadrature:
226template<class T>
227inline
228std::ostream& operator<< (std::ostream& os, const quadrature<T>& q)
229{
230 return os << q.data();
231}
232
233// ------------------------------------------------------------
234// inlined
235// ------------------------------------------------------------
236template<class T>
237inline
242template<class T>
243inline
246{
247 // when using n nodes : gauss quadrature formulae order is r=2*n-1
248 if (r == 0) return 1;
249 size_type n = (r % 2 == 0) ? r/2+1 : (r+1)/2;
250 return std::max(size_t(1), n);
251}
252template<class T>
253inline
256{
257 return _options.get_order();
258}
259template<class T>
260inline
263{
264 return _options.get_family();
265}
266template<class T>
267inline
268std::string
270{
271 return _options.get_family_name();
272}
273template<class T>
274inline
277{
278 return _options;
279}
280#ifdef TO_CLEAN
281template<class T>
282inline
283void
285{
286 if (get_order() != r) {
287 // do not re-initialize nodes-weights if unchanged
288 _options.set_order(r);
289 std::fill (_initialized.begin(), _initialized.end(), false);
290 }
291}
292template<class T>
293inline
294void
296{
297 if (get_family() != ft) {
298 // do not re-initialize nodes-weights if unchanged
299 _options.set_family(ft);
300 std::fill (_initialized.begin(), _initialized.end(), false);
301 }
302}
303#endif // TO_CLEAN
304
305}// namespace rheolef
306#endif // _RHEO_QUADRATURE_H
field::size_type size_type
Definition branch.cc:430
see the integrate_option page for the full documentation
see the persistent_table page for the full documentation
void init_tetrahedron(quadrature_option opt)
void initialize(reference_element hat_K, quadrature_option opt)
friend std::ostream & operator<<(std::ostream &, const quadrature_on_geo< U > &)
void init_prism(quadrature_option opt)
void init_point(quadrature_option opt)
base::size_type size_type
Definition quadrature.h:88
static size_type n_node_gauss(size_type r)
Definition quadrature.h:245
void get_nodes(Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
quadrature_on_geo & operator=(const quadrature_on_geo &q)
Definition quadrature.h:96
void init_hexahedron(quadrature_option opt)
void init_square(quadrature_option opt)
void init_edge(quadrature_option opt)
quadrature_on_geo(const quadrature_on_geo &q)
Definition quadrature.h:95
void init_triangle(quadrature_option opt)
std::vector< weighted_point< T > > base
Definition quadrature.h:87
void wx(const point_basic< T > &x, const T &w)
Definition quadrature.h:116
friend std::ostream & operator<<(std::ostream &, const quadrature_rep< U > &)
geo_element_indirect::orientation_type orientation_type
Definition quadrature.h:136
std::vector< weighted_point< T > >::const_iterator const_iterator
Definition quadrature.h:135
const quadrature_rep & operator=(const quadrature_rep< T > &q)
const_iterator begin(reference_element hat_K) const
quadrature_on_geo< T >::size_type size_type
Definition quadrature.h:133
std::array< quadrature_on_geo< T >, reference_element::max_variant > _quad
Definition quadrature.h:175
const_iterator end(reference_element hat_K) const
void get_nodes(reference_element hat_K, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
void set_family(family_type ft)
Definition quadrature.h:295
quadrature_option::family_type family_type
Definition quadrature.h:134
void set_order(size_type order)
Definition quadrature.h:284
std::string get_family_name() const
Definition quadrature.h:269
size_type get_order() const
Definition quadrature.h:255
std::vector< bool > _initialized
Definition quadrature.h:176
size_type size(reference_element hat_K) const
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
quadrature_option _options
Definition quadrature.h:173
void _initialize(reference_element hat_K) const
static quadrature_rep * make_ptr(const std::string &name)
Definition quadrature.h:178
const quadrature_option & get_options() const
Definition quadrature.h:276
std::string name() const
family_type get_family() const
Definition quadrature.h:262
rep::const_iterator const_iterator
Definition quadrature.h:195
quadrature_rep< T > rep
Definition quadrature.h:191
family_type get_family() const
Definition quadrature.h:215
std::string name() const
Definition quadrature.h:213
const_iterator end(reference_element hat_K) const
Definition quadrature.h:219
void reset(const std::string &name)
void get_nodes(reference_element hat_K, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
Definition quadrature.h:222
void set_family(family_type ft)
rep::orientation_type orientation_type
Definition quadrature.h:196
void set_order(size_type order)
rep::family_type family_type
Definition quadrature.h:194
std::string get_family_name() const
Definition quadrature.h:216
smart_pointer< rep > base
Definition quadrature.h:192
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
Definition quadrature.h:220
size_type size(reference_element hat_K) const
Definition quadrature.h:217
rep::size_type size_type
Definition quadrature.h:193
size_type get_order() const
Definition quadrature.h:214
const quadrature_option & get_options() const
const_iterator begin(reference_element hat_K) const
Definition quadrature.h:218
see the reference_element page for the full documentation
static const variant_type max_variant
see the smart_pointer page for the full documentation
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition catchmark.h:99
STL namespace.
weighted_point(const point_basic< T > &x1, const T &w1)
Definition quadrature.h:73
point_basic< T > x
Definition quadrature.h:74