← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to nfem/plib/basis_on_pointset.h

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito, Pierre Saramito, Sylvestre Ledru
  • Date: 2012-05-14 14:02:09 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20120514140209-dzbdlidkotyflf9e
Tags: 6.1-1
http://bugs.debian.org/671996
[ Pierre Saramito ]
* New upstream release 6.1 (minor changes):
  - support arbitrarily polynomial order Pk
  - source code supports g++-4.7 (closes: #671996)

[ Sylvestre Ledru ]
* update of the watch file

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/plib/basis_on_pointset.h
 
1
#ifndef _RHEO_BASIS_ON_POINTSET_H
 
2
#define _RHEO_BASIS_ON_POINTSET_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
//! Pre-compute a piola transformation polynomial basis
 
25
//! on a lattice definied on the reference element.
 
26
//! The computation is performed one time for all
 
27
//! the first time the reference element appears.
 
28
//! Thus, when looping on the mesh for the Lagrange interpolation
 
29
//! on a Lagrange nodal basis : the lattice is the set of nodes
 
30
//! the evaluation of the fem polynomial basis.
 
31
//! The evaluation is performed only one time.
 
32
 
 
33
#include "rheolef/basis.h"
 
34
#include "rheolef/quadrature.h"
 
35
namespace rheolef {
 
36
 
 
37
// -----------------------------------------------------------------------
 
38
// basis evaluated on lattice of quadrature formulae
 
39
// -----------------------------------------------------------------------
 
40
template<class T>
 
41
class basis_on_pointset {
 
42
public:
 
43
 
 
44
    typedef typename std::vector<T>::size_type                     size_type;
 
45
    typedef typename std::vector<T>::const_iterator                const_iterator;
 
46
    typedef typename std::vector<point_basic<T> >::const_iterator  const_iterator_grad;
 
47
 
 
48
// allocators:
 
49
 
 
50
    basis_on_pointset ();
 
51
    basis_on_pointset (const quadrature<T>& quad, const basis_basic<T>& b);
 
52
    basis_on_pointset (const basis_basic<T>& nb,  const basis_basic<T>& b);
 
53
 
 
54
// modifiers:
 
55
 
 
56
    void set (const quadrature<T>& quad, const basis_basic<T>& b);
 
57
    void set (const basis_basic<T>& nb,  const basis_basic<T>& b);
 
58
 
 
59
// accessors:
 
60
 
 
61
    const basis_basic<T>& get_basis() const { return _b; }
 
62
    size_type size (reference_element hat_K) const;
 
63
 
 
64
    // old interface:
 
65
    const_iterator begin (reference_element hat_K, size_type q) const;
 
66
    const_iterator end   (reference_element hat_K, size_type q) const;
 
67
 
 
68
    const_iterator_grad begin_grad (reference_element hat_K, size_type q) const;
 
69
    const_iterator_grad end_grad   (reference_element hat_K, size_type q) const;
 
70
 
 
71
    // new interface:
 
72
    void evaluate      (reference_element hat_K, size_type q) const;
 
73
    void evaluate_grad (reference_element hat_K, size_type q) const;
 
74
    void evaluate      (reference_element hat_K, const point_basic<T>& hat_xq) const;
 
75
    void evaluate_grad (reference_element hat_K, const point_basic<T>& hat_xq) const;
 
76
    const_iterator      begin() const;
 
77
    const_iterator      end() const;
 
78
    const_iterator_grad begin_grad() const;
 
79
    const_iterator_grad end_grad() const;
 
80
 
 
81
// data:
 
82
    typedef enum {
 
83
      quad_mode     = 0,
 
84
      nodal_mode    = 1,
 
85
      specific_mode = 2,
 
86
      max_mode      = 3
 
87
    } mode_type;
 
88
protected:
 
89
    basis_basic<T>                                      _b;
 
90
    mutable mode_type                                   _mode;
 
91
    const quadrature<T>*                                _p_quad; // when mode: on quadrature pointset
 
92
    basis_basic<T>                                      _nb;     // when mode: on nodal basis pointset
 
93
    mutable std::vector<std::vector<T> >                _val         [reference_element::max_variant];
 
94
    mutable std::vector<std::vector<point_basic<T> > >  _grad_val    [reference_element::max_variant];
 
95
    mutable std::vector<bool>                           _initialized;
 
96
    mutable std::vector<bool>                           _grad_initialized;
 
97
    mutable std::vector<T>                              _val_specific;
 
98
    mutable std::vector<point_basic<T> >                _grad_val_specific;
 
99
    mutable size_type                                   _curr_K_variant;
 
100
    mutable size_type                                   _curr_q;
 
101
    void _initialize          (reference_element hat_K) const;
 
102
    void _grad_initialize     (reference_element hat_K) const;
 
103
};
 
104
// ----------------------------------------------------------
 
105
// inlined
 
106
// ----------------------------------------------------------
 
107
template<class T>
 
108
inline
 
109
basis_on_pointset<T>::basis_on_pointset ()
 
110
 : _b(),
 
111
   _mode(max_mode),
 
112
   _p_quad(0),
 
113
   _nb(),
 
114
   _val(),
 
115
   _grad_val(),
 
116
   _initialized(reference_element::max_variant, false),
 
117
   _grad_initialized(reference_element::max_variant, false),
 
118
   _val_specific(),
 
119
   _curr_K_variant (std::numeric_limits<size_type>::max()),
 
120
   _curr_q (std::numeric_limits<size_type>::max())
 
121
{
 
122
}
 
123
template<class T>
 
124
inline
 
125
basis_on_pointset<T>::basis_on_pointset (const quadrature<T>& quad, const basis_basic<T>& b)
 
126
 : _b(b),
 
127
   _mode(quad_mode),
 
128
   _p_quad(&quad),
 
129
   _nb(),
 
130
   _val(),
 
131
   _grad_val(),
 
132
   _initialized(reference_element::max_variant, false),
 
133
   _grad_initialized(reference_element::max_variant, false),
 
134
   _val_specific(),
 
135
   _curr_K_variant (std::numeric_limits<size_type>::max()),
 
136
   _curr_q (std::numeric_limits<size_type>::max())
 
137
{
 
138
}
 
139
template<class T>
 
140
inline
 
141
basis_on_pointset<T>::basis_on_pointset (const basis_basic<T>& nb, const basis_basic<T>& b)
 
142
 : _b(b),
 
143
   _mode(nodal_mode),
 
144
   _p_quad(0),
 
145
   _nb(nb),
 
146
   _val(),
 
147
   _grad_val(),
 
148
   _initialized(reference_element::max_variant, false),
 
149
   _grad_initialized(reference_element::max_variant, false),
 
150
   _val_specific(),
 
151
   _curr_K_variant (std::numeric_limits<size_type>::max()),
 
152
   _curr_q (std::numeric_limits<size_type>::max())
 
153
{
 
154
}
 
155
template<class T>
 
156
inline
 
157
void
 
158
basis_on_pointset<T>::set (const quadrature<T>& quad, const basis_basic<T>& b)
 
159
{
 
160
    // data was not defined on this quadrature formulae : reset
 
161
    _b  = b;
 
162
    _mode = quad_mode;
 
163
    _p_quad = &quad;
 
164
    std::fill (_initialized.begin(), _initialized.end(), false);
 
165
    std::fill (_grad_initialized.begin(), _grad_initialized.end(), false);
 
166
}
 
167
template<class T>
 
168
inline
 
169
void
 
170
basis_on_pointset<T>::set (const basis_basic<T>& nb, const basis_basic<T>& b)
 
171
{
 
172
    // data was not defined on this quadrature formulae : reset
 
173
    _b  = b;
 
174
    _mode = nodal_mode;
 
175
    _nb = nb;
 
176
    std::fill (_initialized.begin(), _initialized.end(), false);
 
177
    std::fill (_grad_initialized.begin(), _grad_initialized.end(), false);
 
178
}
 
179
// ==================================================================================
 
180
// old interface
 
181
// ==================================================================================
 
182
template<class T>
 
183
inline
 
184
typename basis_on_pointset<T>::const_iterator
 
185
basis_on_pointset<T>::begin (reference_element hat_K, size_type q) const
 
186
{
 
187
    reference_element::variant_type K_variant = hat_K.variant();
 
188
    if (!_initialized [K_variant]) _initialize (hat_K);
 
189
    return _val[K_variant][q].begin();
 
190
}
 
191
template<class T>
 
192
inline
 
193
typename basis_on_pointset<T>::const_iterator
 
194
basis_on_pointset<T>::end (reference_element hat_K, size_type q) const
 
195
{
 
196
    reference_element::variant_type K_variant = hat_K.variant();
 
197
    if (!_initialized [K_variant]) _initialize (hat_K);
 
198
    return _val[K_variant][q].end();
 
199
}
 
200
template<class T>
 
201
inline
 
202
typename basis_on_pointset<T>::const_iterator_grad
 
203
basis_on_pointset<T>::begin_grad (reference_element hat_K, size_type q) const
 
204
{
 
205
    reference_element::variant_type K_variant = hat_K.variant();
 
206
    if (!_grad_initialized [K_variant]) _grad_initialize (hat_K);
 
207
    return _grad_val[K_variant][q].begin();
 
208
}
 
209
template<class T>
 
210
inline
 
211
typename basis_on_pointset<T>::const_iterator_grad
 
212
basis_on_pointset<T>::end_grad (reference_element hat_K, size_type q) const
 
213
{
 
214
    reference_element::variant_type K_variant = hat_K.variant();
 
215
    if (!_grad_initialized [K_variant]) _grad_initialize (hat_K);
 
216
    return _grad_val[K_variant][q].end();
 
217
}
 
218
// ==================================================================================
 
219
// new interface
 
220
// ==================================================================================
 
221
// -------------------
 
222
// scalar
 
223
// -------------------
 
224
template<class T>
 
225
inline
 
226
void
 
227
basis_on_pointset<T>::evaluate (reference_element hat_K, size_type q) const
 
228
{
 
229
    reference_element::variant_type K_variant = hat_K.variant();
 
230
    if (!_initialized [K_variant]) _initialize (hat_K);
 
231
    _curr_K_variant = K_variant;
 
232
    _curr_q     = q;
 
233
}
 
234
template<class T>
 
235
inline
 
236
void
 
237
basis_on_pointset<T>::evaluate (reference_element hat_K, const point_basic<T>& hat_xq) const
 
238
{
 
239
    size_type nb = _b.size(hat_K);
 
240
    if (nb != _val_specific.size ()) { _val_specific.resize (nb); }
 
241
    _b.eval (hat_K, hat_xq, _val_specific);
 
242
    _mode = specific_mode;
 
243
}
 
244
template<class T>
 
245
inline
 
246
typename basis_on_pointset<T>::const_iterator
 
247
basis_on_pointset<T>::begin() const
 
248
{
 
249
    if (_mode != specific_mode) return _val[_curr_K_variant][_curr_q].begin();
 
250
    else                        return _val_specific.begin();
 
251
}
 
252
template<class T>
 
253
inline
 
254
typename basis_on_pointset<T>::const_iterator
 
255
basis_on_pointset<T>::end() const
 
256
{
 
257
    if (_mode != specific_mode) return _val[_curr_K_variant][_curr_q].end();
 
258
    else                        return _val_specific.end();
 
259
}
 
260
// -------------------
 
261
// grad
 
262
// -------------------
 
263
template<class T>
 
264
inline
 
265
void
 
266
basis_on_pointset<T>::evaluate_grad (reference_element hat_K, size_type q) const
 
267
{
 
268
    reference_element::variant_type K_variant = hat_K.variant();
 
269
    if (!_grad_initialized [K_variant]) _grad_initialize (hat_K);
 
270
    _curr_K_variant = K_variant;
 
271
    _curr_q     = q;
 
272
}
 
273
template<class T>
 
274
inline
 
275
void
 
276
basis_on_pointset<T>::evaluate_grad (reference_element hat_K, const point_basic<T>& hat_xq) const
 
277
{
 
278
    size_type nb = _b.size(hat_K);
 
279
    if (nb != _grad_val_specific.size ()) { _grad_val_specific.resize (nb); }
 
280
    _b.grad_eval (hat_K, hat_xq, _grad_val_specific);
 
281
    _mode = specific_mode;
 
282
}
 
283
template<class T>
 
284
inline
 
285
typename basis_on_pointset<T>::const_iterator_grad
 
286
basis_on_pointset<T>::begin_grad() const
 
287
{
 
288
    if (_mode != specific_mode) return _grad_val[_curr_K_variant][_curr_q].begin();
 
289
    else                        return _grad_val_specific.begin();
 
290
}
 
291
template<class T>
 
292
inline
 
293
typename basis_on_pointset<T>::const_iterator_grad
 
294
basis_on_pointset<T>::end_grad() const
 
295
{
 
296
    if (_mode != specific_mode) return _grad_val[_curr_K_variant][_curr_q].end();
 
297
    else                        return _grad_val_specific.end();
 
298
}
 
299
 
 
300
}// namespace rheolef
 
301
#endif // _RHEO_BASIS_ON_POINTSET_H