Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
quadrature_rep.cc
Go to the documentation of this file.
1
21#include "rheolef/quadrature.h"
22#include "rheolef/reference_element_face_transformation.h"
23#include "rheolef/rheostream.h"
24#include <cctype> // for std::isdigit
25namespace rheolef {
26using namespace std;
27
28// ----------------------------------------------------------------------------
29// quadrature options
30// ----------------------------------------------------------------------------
31static const char* static_family_name[quadrature_option::max_family+1] = {
32 "gauss" ,
33 "gauss_lobatto" ,
34 "gauss_radau" ,
35 "middle_edge" ,
36 "superconvergent" ,
37 "equispaced" ,
38 "undefined"
39};
40static
41quadrature_option::family_type
42family_name2type (string name)
43{
44 for (size_t i = 0; i < quadrature_option::max_family; ++i) {
45 if (static_family_name[i] == name) {
46 return (quadrature_option::family_type) (i);
47 }
48 }
49 error_macro ("unexpected quadrature family name `" << name << "'");
50 return quadrature_option::max_family;
51}
52void
53quadrature_option::set_family (string name)
54{
55 set_family (family_name2type (name));
56}
57string
58quadrature_option::get_family_name() const
59{
60 check_macro (_family >= 0 && _family <= max_family,
61 "unexpected quadrature family number = " << _family);
62 return static_family_name[_family];
63}
64string
65quadrature_option::name() const
66{
67 if (get_family() == default_family && get_order() == default_order) {
68 return "";
69 } else {
70 return get_family_name() + "(" + std::to_string(get_order()) + ")";
71 }
72#ifdef TO_CLEAN
73 return get_family_name() + "(" + std::to_string(long(get_order())) + ")";
74#endif // TO_CLEAN
75}
76void
77quadrature_option::reset (const std::string& name)
78{
79 if (name == "") {
80 // set options to default:
81 set_family (default_family);
82 set_order (default_order);
83 return;
84 }
85 size_type i = name.find ('(');
86 string family = name.substr (0, i);
87 size_type i1 = i+1, l1 = name.size()-i1-1;
88 check_macro (name.size() > i1+1,
89 "invalid quadrature name \""<<name<<"\"; expect e.g. \"gauss(3)\"");
90 string order_str = name.substr (i1, l1);
91 for (size_type i = 0; i < order_str.size(); ++i) {
92 // accepts also e.g. "gauss(-1)" as default initialization
93 check_macro (std::isdigit(order_str[i]) || order_str[i] == '-',
94 "invalid quadrature name \""<<name<<"\"; expect e.g. \"gauss(3)\"");
95 }
96 size_type order = atoi(order_str.c_str());
97 set_family (family);
98 set_order (order);
99}
100// ----------------------------------------------------------------------------
101// quadrature on a specific reference element
102// ----------------------------------------------------------------------------
103template<class T>
104void
106{
107 base::resize (0);
108 switch (hat_K.variant()) {
109 case reference_element::p: init_point (opt); break;
110 case reference_element::e: init_edge (opt); break;
111 case reference_element::t: init_triangle (opt); break;
112 case reference_element::q: init_square (opt); break;
113 case reference_element::T: init_tetrahedron (opt); break;
114 case reference_element::P: init_prism (opt); break;
115 case reference_element::H: init_hexahedron (opt); break;
116 default:
117 fatal_macro ("quadrature formula on element type `"
118 << hat_K.name() << "' are not yet implemented.");
119 }
120}
121template<class T>
122void
123quadrature_on_geo<T>::get_nodes (Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const
124{
125 size_type nq = base::size();
126 node.resize (nq);
127 for (size_type q = 0; q < nq; ++q) {
128 node[q] = base::operator[](q).x;
129 }
130}
131// ----------------------------------------------------------------------------
132// quadrature on the full set of reference elements
133// ----------------------------------------------------------------------------
134template<class T>
136 : _options(opt),
137 _quad(),
138 _initialized (reference_element::max_variant, false)
139{
140}
141template<class T>
142quadrature_rep<T>::quadrature_rep (const std::string& name)
143 : _options(name),
144 _quad(),
145 _initialized (reference_element::max_variant, false)
146{
147}
148template<class T>
150 : _options (q._options),
151 _quad (q._quad),
152 _initialized (q._initialized)
153{
154}
155template<class T>
158 _options = q._options;
159 _quad = q._quad;
160 _initialized = q._initialized;
161 return *this;
162}
163template<class T>
164string
166{
167 return _options.name();
168}
169template<class T>
170void
172{
174 _quad[K_type].initialize (hat_K, _options);
175 _initialized [K_type] = true;
176}
177template<class T>
178void
179quadrature_rep<T>::get_nodes (reference_element hat_K, Eigen::Matrix<point_basic<T>, Eigen::Dynamic, 1>& node) const
180{
181 if (!_initialized [hat_K.variant()]) _initialize (hat_K);
182 _quad [hat_K.variant()].get_nodes (node);
183}
184template<class T>
187{
189 if (!_initialized [K_type]) _initialize (hat_K);
190 return _quad[K_type].begin();
191}
192template<class T>
195{
197 if (!_initialized [K_type]) _initialize (hat_K);
198 return _quad[K_type].end();
199}
200template<class T>
203{
205 if (!_initialized [K_type]) _initialize (hat_K);
206 return _quad[K_type].size();
207}
208template<class T>
211{
213 if (!_initialized [K_type]) _initialize (hat_K);
214 return _quad[K_type][q];
215}
216template<class T>
217ostream&
218operator<< (ostream& out, const quadrature_on_geo<T>& x)
219{
220 out << setprecision (numeric_limits<T>::digits10)
221 << x.size() << endl;
222 for (size_t r = 0; r < x.size(); r++)
223 out << x[r].x << "\t" << x[r].w << endl;
224 return out;
225}
226template<class T>
227ostream&
228operator<< (ostream& out, const quadrature_rep<T>& x)
229{
230 out << "quadrature" << endl
231 << x._options.get_family_name() << " " << x._options.get_order() << endl;
232 for (size_t i = 0; i != size_t(reference_element::max_variant); i++) {
234 if (! x._initialized [K_type]) continue;
235 reference_element hat_K (K_type);
236 out << "reference_element " << hat_K.name() << endl
237 << x._quad[K_type];
238 }
239 out << "end_quadrature" << endl;
240 return out;
241}
242// ----------------------------------------------------------------------------
243// quadrature : pointer interface
244// ----------------------------------------------------------------------------
245template <class T>
247{
248 if (_options.name() != "") {
249 persistent_table<quadrature<T>>::unload (name());
250 }
251}
252template<class T>
253void
254quadrature<T>::reset (const std::string& name)
255{
256 if (name == "") {
257 base::operator= (new_macro(quadrature_rep<T>()));
258 } else {
259 base::operator= (persistent_table<quadrature<T>>::load (name));
260 }
261}
262template<class T>
265{
266 return base::data().get_options();
267}
268template<class T>
269quadrature<T>::quadrature (quadrature_option opt)
270 : base(),
272{
273 reset (opt.name());
274}
275template<class T>
276quadrature<T>::quadrature (const std::string& name)
277 : base(),
279{
280 reset (name);
281}
282template<class T>
283void
285{
286 quadrature_option qopt = get_options();
287 qopt.set_order (order);
288 reset (qopt.name());
289}
290template<class T>
291void
293{
294 quadrature_option qopt = get_options();
295 qopt.set_family (ft);
296 reset (qopt.name());
297}
298// ----------------------------------------------------------------------------
299// instanciation in library
300// ----------------------------------------------------------------------------
301#define _RHEOLEF_instanciation(T) \
302template class quadrature_on_geo<Float>; \
303template class quadrature_rep<Float>; \
304template class quadrature<T>; \
305template ostream& operator<< (ostream&, const quadrature_on_geo<T>&); \
306template ostream& operator<< (ostream&, const quadrature_rep<T>&);
307
309
310} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the integrate_option page for the full documentation
see the persistent_table page for the full documentation
void initialize(reference_element hat_K, quadrature_option opt)
base::size_type size_type
Definition quadrature.h:88
void get_nodes(Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &node) const
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
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
quadrature_rep(quadrature_option opt=quadrature_option())
std::string name() const
std::string name() const
Definition quadrature.h:213
void reset(const std::string &name)
void set_family(family_type ft)
void set_order(size_type order)
rep::family_type family_type
Definition quadrature.h:194
quadrature(const std::string &name="")
rep::size_type size_type
Definition quadrature.h:193
const quadrature_option & get_options() const
see the reference_element page for the full documentation
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type max_variant
static const variant_type p
variant_type variant() const
static const variant_type T
static const variant_type P
static const variant_type t
see the smart_pointer page for the full documentation
static const char * static_family_name[quadrature_option::max_family+1]
#define error_macro(message)
Definition dis_macros.h:49
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition catchmark.h:99
STL namespace.
void load(idiststream &in, Float &p, field &uh)