Rheolef  7.1
an efficient C++ finite element environment
quadrature.h
Go to the documentation of this file.
1 #ifndef _RHEO_QUADRATURE_H
2 #define _RHEO_QUADRATURE_H
3 #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 
31 namespace rheolef {
32 
33 /*Class:quadrature
34 NAME: @code{quadrature} - quadrature formulae on the reference lement
35 @cindex quadrature formulae
36 @clindex quadrature
37 @clindex reference_element
38 SYNOPSIS:
39 @noindent
40 The @code{quadrature} class defines a container for a quadrature
41 formulae on the reference element (@pxref{reference_element iclass}).
42 This container stores the nodes coordinates and the weights.
43 The constructor takes two arguments: the reference element
44 @math{K}
45 and the order @math{r} of the quadrature formulae.
46 The formulae is exact when computing the integral
47 of 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
55 LIMITATIONS:
56 @noindent
57 The formulae is optimal when it uses a minimal number
58 of nodes @math{n}.
59 Optimal quadrature formula are hard-coded in this class.
60 Not all reference elements and orders are yet
61 implemented. This class will be completed in the future.
62 
63 AUTHORS:
64  LMC-IMAG, 38041 Grenoble cedex 9, France
65  | Pierre.Saramito@imag.fr
66 DATE: 30 november 2003
67 End:
68 */
69 
70 template<class T>
72  weighted_point () : x(), w(0) {}
73  weighted_point (const point_basic<T>& x1, const T& w1) : x(x1), w(w1) {}
75  T w;
76 };
77 
78 // ----------------------------------------------------------------------------
79 // quadrature on a specific ref element
80 // ----------------------------------------------------------------------------
81 template<class T>
82 class quadrature_on_geo : public std::vector<weighted_point<T> > {
83 public:
84 
85 // typedefs:
86 
87  typedef std::vector<weighted_point<T> > base;
88  typedef typename base::size_type size_type;
89  typedef point_basic<T> x;
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);
110  void init_triangle (quadrature_option opt);
111  void init_square (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 
119  static size_type n_node_gauss (size_type r);
120 
121  template<class U>
122  friend std::ostream& operator<< (std::ostream&, const quadrature_on_geo<U>&);
123 };
124 // ----------------------------------------------------------------------------
125 // quadrature representation
126 // ----------------------------------------------------------------------------
127 template<class T>
129 public:
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 
140  ~quadrature_rep();
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;
161  const_iterator begin (reference_element hat_K) const;
162  const_iterator end (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 
169 protected:
170 // internal:
171  void _initialize (reference_element hat_K) const;
172 // data:
174  mutable std::array<quadrature_on_geo<T>,reference_element::max_variant>
176  mutable std::vector<bool> _initialized;
177 public:
178  static quadrature_rep* make_ptr (const std::string& name) { return new_macro(quadrature_rep(name)); }
179 };
180 // ----------------------------------------------------------------------------
181 // quadrature class
182 // ----------------------------------------------------------------------------
183 //<quadrature:
184 template<class T>
185 class quadrature : public smart_pointer<quadrature_rep<T> >,
186  public persistent_table<quadrature<T>> {
187 public:
188 
189 // typedefs:
190 
193  typedef typename rep::size_type size_type;
194  typedef typename rep::family_type family_type;
197 
198 // allocators:
199 
200 
201  quadrature (const std::string& name = "");
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:
226 template<class T>
227 inline
228 std::ostream& operator<< (std::ostream& os, const quadrature<T>& q)
229 {
230  return os << q.data();
231 }
232 
233 // ------------------------------------------------------------
234 // inlined
235 // ------------------------------------------------------------
236 template<class T>
237 inline
239  : std::vector<weighted_point<T> >()
240 {
241 }
242 template<class T>
243 inline
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 }
252 template<class T>
253 inline
256 {
257  return _options.get_order();
258 }
259 template<class T>
260 inline
263 {
264  return _options.get_family();
265 }
266 template<class T>
267 inline
268 std::string
270 {
271  return _options.get_family_name();
272 }
273 template<class T>
274 inline
275 const quadrature_option&
277 {
278  return _options;
279 }
280 #ifdef TO_CLEAN
281 template<class T>
282 inline
283 void
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 }
292 template<class T>
293 inline
294 void
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
rheolef::quadrature_on_geo::quadrature_on_geo
quadrature_on_geo(const quadrature_on_geo &q)
Definition: quadrature.h:95
rheolef::quadrature_rep::quadrature_rep
quadrature_rep(quadrature_option opt=quadrature_option())
Definition: quadrature_rep.cc:135
rheolef::quadrature_on_geo::init_triangle
void init_triangle(quadrature_option opt)
Definition: quadrature_t.cc:50
rheolef::quadrature_on_geo::init_edge
void init_edge(quadrature_option opt)
Definition: quadrature_e.cc:38
rheolef::weighted_point
Definition: quadrature.h:71
rheolef::quadrature::reset
void reset(const std::string &name)
Definition: quadrature_rep.cc:254
rheolef::point_basic
Definition: point.h:87
rheolef::quadrature::size
size_type size(reference_element hat_K) const
Definition: quadrature.h:217
rheolef::quadrature::operator()
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
Definition: quadrature.h:220
rheolef::quadrature_on_geo
Definition: quadrature.h:82
rheolef::quadrature_rep
Definition: quadrature.h:128
rheolef::quadrature_on_geo::operator=
quadrature_on_geo & operator=(const quadrature_on_geo &q)
Definition: quadrature.h:96
rheolef::quadrature_on_geo::operator<<
friend std::ostream & operator<<(std::ostream &, const quadrature_on_geo< U > &)
rheolef::quadrature_rep::operator()
const weighted_point< T > & operator()(reference_element hat_K, size_type q) const
Definition: quadrature_rep.cc:210
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::base
smart_pointer< rep > base
Definition: quadrature.h:192
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::set_order
void set_order(size_type order)
Definition: quadrature_rep.cc:284
rheolef::quadrature::const_iterator
rep::const_iterator const_iterator
Definition: quadrature.h:195
rheolef::quadrature_rep::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition: quadrature.h:136
rheolef::quadrature_rep::operator=
const quadrature_rep & operator=(const quadrature_rep< T > &q)
Definition: quadrature_rep.cc:157
rheolef::quadrature_on_geo::initialize
void initialize(reference_element hat_K, quadrature_option opt)
Definition: quadrature_rep.cc:105
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
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_rep::get_options
const quadrature_option & get_options() const
Definition: quadrature.h:276
rheolef::weighted_point::weighted_point
weighted_point(const point_basic< T > &x1, const T &w1)
Definition: quadrature.h:73
rheolef::quadrature_on_geo::init_point
void init_point(quadrature_option opt)
Definition: quadrature_p.cc:27
rheolef::quadrature::quadrature
quadrature(const std::string &name="")
Definition: quadrature_rep.cc:276
rheolef::quadrature::get_family
family_type get_family() const
Definition: quadrature.h:215
rheolef::quadrature_rep::family_type
quadrature_option::family_type family_type
Definition: quadrature.h:134
rheolef::quadrature_rep::get_order
size_type get_order() const
Definition: quadrature.h:255
rheolef::quadrature::end
const_iterator end(reference_element hat_K) const
Definition: quadrature.h:219
rheolef::quadrature::family_type
rep::family_type family_type
Definition: quadrature.h:194
rheolef::quadrature::get_options
const quadrature_option & get_options() const
Definition: quadrature_rep.cc:264
rheolef::persistent_table
see the persistent_table page for the full documentation
Definition: persistent_table.h:84
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::smart_pointer
see the smart_pointer page for the full documentation
Definition: smart_pointer.h:351
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::quadrature_rep::operator<<
friend std::ostream & operator<<(std::ostream &, const quadrature_rep< U > &)
rheolef::quadrature_on_geo::x
point_basic< T > x
Definition: quadrature.h:89
rheolef::quadrature_on_geo::init_square
void init_square(quadrature_option opt)
Definition: quadrature_q.cc:46
rheolef::weighted_point::w
T w
Definition: quadrature.h:75
rheolef::quadrature::get_order
size_type get_order() const
Definition: quadrature.h:214
rheolef::quadrature_rep::set_order
void set_order(size_type order)
Definition: quadrature.h:284
rheolef::quadrature::name
std::string name() const
Definition: quadrature.h:213
rheolef::weighted_point::x
point_basic< T > x
Definition: quadrature.h:74
rheolef::geo_element_indirect::orientation_type
short int orientation_type
Definition: geo_element_indirect.h:38
rheolef::quadrature_on_geo::quadrature_on_geo
quadrature_on_geo()
Definition: quadrature.h:238
rheolef::quadrature::orientation_type
rep::orientation_type orientation_type
Definition: quadrature.h:196
rheolef::quadrature::get_family_name
std::string get_family_name() const
Definition: quadrature.h:216
rheolef::quadrature_rep::_quad
std::array< quadrature_on_geo< T >, reference_element::max_variant > _quad
Definition: quadrature.h:175
rheolef::quadrature_rep::set_family
void set_family(family_type ft)
Definition: quadrature.h:295
rheolef::quadrature_rep::get_family
family_type get_family() const
Definition: quadrature.h:262
rheolef::quadrature::set_family
void set_family(family_type ft)
Definition: quadrature_rep.cc:292
rheolef::weighted_point::weighted_point
weighted_point()
Definition: quadrature.h:72
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::smart_pointer_base< T, details::constructor_copy< T > >::data
const T & data() const
Definition: smart_pointer.h:266
rheolef::quadrature_on_geo::size_type
base::size_type size_type
Definition: quadrature.h:88
rheolef::quadrature_rep::get_family_name
std::string get_family_name() const
Definition: quadrature.h:269
rheolef::quadrature_rep::size
size_type size(reference_element hat_K) const
Definition: quadrature_rep.cc:202
rheolef::quadrature_rep::_initialized
std::vector< bool > _initialized
Definition: quadrature.h:176
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::quadrature_rep::_initialize
void _initialize(reference_element hat_K) const
Definition: quadrature_rep.cc:171
rheolef::quadrature_rep::begin
const_iterator begin(reference_element hat_K) const
Definition: quadrature_rep.cc:186
rheolef::quadrature_on_geo::init_tetrahedron
void init_tetrahedron(quadrature_option opt)
Definition: quadrature_Te.cc:28
rheolef::quadrature::rep
quadrature_rep< T > rep
Definition: quadrature.h:191
rheolef::quadrature_on_geo::init_hexahedron
void init_hexahedron(quadrature_option opt)
Definition: quadrature_H.cc:28
rheolef::space_constant::vector
Definition: space_constant.h:137
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::quadrature_rep::_options
quadrature_option _options
Definition: quadrature.h:173
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::quadrature_rep::end
const_iterator end(reference_element hat_K) const
Definition: quadrature_rep.cc:194
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
rheolef::quadrature_on_geo::base
std::vector< weighted_point< T > > base
Definition: quadrature.h:87
rheolef::quadrature_option
integrate_option quadrature_option
Definition: integrate_option.h:186
rheolef::std
Definition: vec_expr_v2.h:391
rheolef::quadrature
Definition: quadrature.h:185
rheolef::quadrature_rep::name
std::string name() const
Definition: quadrature_rep.cc:165
rheolef::quadrature_rep::size_type
quadrature_on_geo< T >::size_type size_type
Definition: quadrature.h:133
rheolef::quadrature_rep::const_iterator
std::vector< weighted_point< T > >::const_iterator const_iterator
Definition: quadrature.h:135
rheolef::quadrature::begin
const_iterator begin(reference_element hat_K) const
Definition: quadrature.h:218
rheolef::quadrature_on_geo::n_node_gauss
static size_type n_node_gauss(size_type r)
Definition: quadrature.h:245
rheolef::quadrature_rep::~quadrature_rep
~quadrature_rep()
Definition: quadrature_rep.cc:246
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::quadrature_on_geo::init_prism
void init_prism(quadrature_option opt)
Definition: quadrature_Pr.cc:28
rheolef::quadrature_on_geo::wx
void wx(const point_basic< T > &x, const T &w)
Definition: quadrature.h:116
rheolef::quadrature::size_type
rep::size_type size_type
Definition: quadrature.h:193
rheolef::quadrature_rep::make_ptr
static quadrature_rep * make_ptr(const std::string &name)
Definition: quadrature.h:178