← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty-proposed

« back to all changes in this revision

Viewing changes to nfem/basis/quadrature.h

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/basis/quadrature.h
 
1
#ifndef _RHEO_QUADRATURE_H
 
2
#define _RHEO_QUADRATURE_H
 
3
///
 
4
/// This file is part of Rheolef.
 
5
///
 
6
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
7
///
 
8
/// Rheolef is free software; you can redistribute it and/or modify
 
9
/// it under the terms of the GNU General Public License as published by
 
10
/// the Free Software Foundation; either version 2 of the License, or
 
11
/// (at your option) any later version.
 
12
///
 
13
/// Rheolef is distributed in the hope that it will be useful,
 
14
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
15
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
16
/// GNU General Public License for more details.
 
17
///
 
18
/// You should have received a copy of the GNU General Public License
 
19
/// along with Rheolef; if not, write to the Free Software
 
20
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
21
/// 
 
22
/// =========================================================================
 
23
 
 
24
/*Class:quadrature
 
25
NAME: @code{quadrature} - quadrature formulae on the reference lement
 
26
@cindex  quadrature formulae
 
27
@clindex quadrature
 
28
@clindex reference_element
 
29
SYNOPSYS:
 
30
@noindent
 
31
The @code{quadrature} class defines a container for a quadrature
 
32
formulae on the reference element (@pxref{reference_element internal}).
 
33
This container stores the nodes coordinates and the weights.
 
34
The constructor takes two arguments: the reference element 
 
35
@math{K}
 
36
and the order @math{r} of the quadrature formulae.
 
37
The formulae is exact when computing the integral
 
38
of a polynom @math{p} that degree is less or equal to order @math{r}.
 
39
@example
 
40
                  n
 
41
    /            ___
 
42
    | p(x) dx =  \    p(x_q) w_q
 
43
    / K          /__
 
44
                 q=1
 
45
@end example
 
46
LIMITATIONS:
 
47
@noindent
 
48
The formulae is optimal when it uses a minimal number
 
49
of nodes @math{n}.
 
50
Optimal quadrature formula are hard-coded in this class.
 
51
Not all reference elements and orders are yet 
 
52
implemented. This class will be completed in the future.
 
53
 
 
54
AUTHORS:
 
55
    LMC-IMAG, 38041 Grenoble cedex 9, France
 
56
   | Pierre.Saramito@imag.fr
 
57
DATE:   30 november 2003
 
58
End:
 
59
*/
 
60
 
 
61
#include "rheolef/reference_element.h"
 
62
#include "rheolef/point.h"
 
63
 
 
64
struct weighted_point {
 
65
    weighted_point ();
 
66
    weighted_point (const point& x, Float w);
 
67
    point x;
 
68
    Float w;
 
69
};
 
70
 
 
71
class quadrature_option_type {
 
72
public:
 
73
    typedef enum {
 
74
        gauss           = 0,
 
75
        gauss_lobatto   = 1,
 
76
        gauss_radau     = 2,
 
77
        middle_edge     = 3,
 
78
        superconvergent = 4,
 
79
        max_family      = 5
 
80
    } family_type; // update also family_name[] in quatrature.cc
 
81
 
 
82
    typedef size_t size_type;
 
83
    quadrature_option_type(
 
84
        family_type ft = quadrature_option_type::gauss,
 
85
        size_type k = std::numeric_limits<size_type>::max());
 
86
    size_t         get_order() const;
 
87
    family_type    get_family() const;
 
88
    std::string    get_family_name() const;
 
89
    void set_order (size_t r);
 
90
    void set_family (family_type ft);
 
91
protected:
 
92
    family_type   _family;
 
93
    size_t        _order;
 
94
};
 
95
class quadrature_on_geo : public std::vector<weighted_point> {
 
96
public:
 
97
 
 
98
// typedefs:
 
99
 
 
100
    typedef std::vector<weighted_point>:: size_type size_type;
 
101
    typedef point                                   x;
 
102
 
 
103
// alocators/deallocators:
 
104
 
 
105
    quadrature_on_geo ();
 
106
 
 
107
// modifier:
 
108
 
 
109
    void initialize (reference_element hat_K, quadrature_option_type opt);
 
110
    void init_point (quadrature_option_type opt);
 
111
    void init_edge (quadrature_option_type opt);
 
112
    void init_triangle (quadrature_option_type opt);
 
113
    void init_square (quadrature_option_type opt);
 
114
    void init_tetrahedron (quadrature_option_type opt);
 
115
    void init_prism (quadrature_option_type opt);
 
116
    void init_hexahedron (quadrature_option_type opt);
 
117
 
 
118
    void wx (const point& x, const Float& w) {
 
119
            push_back (weighted_point(x,w)); }
 
120
 
 
121
    static size_type n_node_gauss (size_t r);
 
122
 
 
123
    friend std::ostream& operator<< (std::ostream&, const quadrature_on_geo&);
 
124
private:
 
125
    quadrature_on_geo (const quadrature_on_geo&);
 
126
    quadrature_on_geo operator= (const quadrature_on_geo&);
 
127
};
 
128
//<quadrature:
 
129
class quadrature {
 
130
public:
 
131
 
 
132
// typedefs:
 
133
 
 
134
    typedef quadrature_on_geo::size_type size_type;
 
135
    typedef quadrature_option_type::family_type family_type;
 
136
    typedef std::vector<weighted_point>::const_iterator const_iterator;
 
137
 
 
138
// allocators:
 
139
 
 
140
    quadrature (quadrature_option_type opt = quadrature_option_type());
 
141
 
 
142
// modifiers:
 
143
 
 
144
    void set_order (size_type order);
 
145
    void set_family (family_type ft);
 
146
 
 
147
// accessors:
 
148
 
 
149
    size_type      get_order() const;
 
150
    family_type    get_family() const;
 
151
    std::string    get_family_name() const;
 
152
    size_type      size  (reference_element hat_K) const;
 
153
    const_iterator begin (reference_element hat_K) const;
 
154
    const_iterator end   (reference_element hat_K) const;
 
155
    friend std::ostream& operator<< (std::ostream&, const quadrature&);
 
156
protected:
 
157
    quadrature_option_type     _options;
 
158
    mutable quadrature_on_geo  _quad        [reference_element::max_size];
 
159
    mutable std::vector<bool>  _initialized;
 
160
    void _initialize (reference_element hat_K) const;
 
161
private:
 
162
    quadrature (const quadrature&);
 
163
    quadrature operator= (const quadrature&);
 
164
};
 
165
//>quadrature:
 
166
// ------------------------------------------------------------
 
167
// inlined
 
168
// ------------------------------------------------------------
 
169
inline 
 
170
quadrature_option_type::quadrature_option_type (family_type ft, size_type k)
 
171
 : _family(ft),
 
172
   _order(k)
 
173
{
 
174
}
 
175
inline
 
176
quadrature_option_type::size_type
 
177
quadrature_option_type::get_order () const
 
178
{
 
179
    return _order;
 
180
}
 
181
inline
 
182
quadrature_option_type::family_type
 
183
quadrature_option_type::get_family () const
 
184
{
 
185
    return _family;
 
186
}
 
187
inline
 
188
void
 
189
quadrature_option_type::set_order (size_t r)
 
190
{
 
191
    _order = r;
 
192
}
 
193
inline
 
194
void
 
195
quadrature_option_type::set_family (family_type ft)
 
196
{
 
197
    _family = ft;
 
198
}
 
199
inline
 
200
weighted_point::weighted_point ()
 
201
 : x(), w(0) {}
 
202
 
 
203
inline
 
204
weighted_point::weighted_point (const point& x1, Float w1)
 
205
 : x(x1), w(w1) {}
 
206
 
 
207
inline 
 
208
quadrature_on_geo::quadrature_on_geo ()
 
209
 : std::vector<weighted_point>()
 
210
{
 
211
}
 
212
inline 
 
213
quadrature::size_type
 
214
quadrature::get_order () const
 
215
{
 
216
    return _options.get_order();
 
217
}
 
218
inline 
 
219
quadrature::family_type
 
220
quadrature::get_family () const
 
221
{
 
222
    return _options.get_family();
 
223
}
 
224
inline 
 
225
std::string
 
226
quadrature::get_family_name () const
 
227
{
 
228
    return _options.get_family_name();
 
229
}
 
230
inline 
 
231
quadrature::quadrature (quadrature_option_type opt)
 
232
 : _options(opt),
 
233
   _quad(),
 
234
   _initialized(reference_element::max_size, false)
 
235
{
 
236
}
 
237
 
 
238
inline
 
239
void
 
240
quadrature::set_order (size_type r)
 
241
{
 
242
  if (get_order() != r) {
 
243
    // do not re-initialize nodes-weights if unchanged
 
244
    _options.set_order(r);
 
245
    std::fill (_initialized.begin(), _initialized.end(), false);
 
246
  }
 
247
}
 
248
inline
 
249
void
 
250
quadrature::set_family (family_type ft)
 
251
{
 
252
  if (get_family() != ft) {
 
253
    // do not re-initialize nodes-weights if unchanged
 
254
    _options.set_family(ft);
 
255
    std::fill (_initialized.begin(), _initialized.end(), false);
 
256
  }
 
257
}
 
258
inline
 
259
quadrature_on_geo::size_type
 
260
quadrature_on_geo::n_node_gauss (size_type r)
 
261
{
 
262
    // when using n nodes : gauss quadrature formulae order is r=2*n-1
 
263
    if (r == 0) return 1;
 
264
    size_type n = (r % 2 == 0) ? r/2+1 : (r+1)/2;
 
265
    return std::max(size_t(1), n);
 
266
}
 
267
#endif // _RHEO_QUADRATURE_H
 
268