Rheolef  7.1
an efficient C++ finite element environment
quadrature_rep.cc
Go to the documentation of this file.
1 #include "rheolef/quadrature.h"
22 #include "rheolef/reference_element_face_transformation.h"
23 #include "rheolef/rheostream.h"
24 #include <cctype> // for std::isdigit
25 namespace rheolef {
26 using namespace std;
27 
28 // ----------------------------------------------------------------------------
29 // quadrature options
30 // ----------------------------------------------------------------------------
31 static 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 };
40 static
41 quadrature_option::family_type
42 family_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 }
52 void
53 quadrature_option::set_family (string name)
54 {
55  set_family (family_name2type (name));
56 }
57 string
58 quadrature_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 }
64 string
66 {
67  if (get_family() == default_family && get_order() == default_order) {
68  return "";
69  } else {
70  return get_family_name() + "(" + itos(get_order()) + ")";
71  }
72 #ifdef TO_CLEAN
73  return get_family_name() + "(" + itos(long(get_order())) + ")";
74 #endif // TO_CLEAN
75 }
76 void
77 quadrature_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 // ----------------------------------------------------------------------------
103 template<class T>
104 void
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 }
121 template<class T>
122 void
123 quadrature_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 // ----------------------------------------------------------------------------
134 template<class T>
136  : _options(opt),
137  _quad(),
138  _initialized (reference_element::max_variant, false)
139 {
140 }
141 template<class T>
143  : _options(name),
144  _quad(),
145  _initialized (reference_element::max_variant, false)
146 {
147 }
148 template<class T>
150  : _options (q._options),
151  _quad (q._quad),
152  _initialized (q._initialized)
153 {
154 }
155 template<class T>
156 const quadrature_rep<T>&
158  _options = q._options;
159  _quad = q._quad;
160  _initialized = q._initialized;
161  return *this;
162 }
163 template<class T>
164 string
166 {
167  return _options.name();
168 }
169 template<class T>
170 void
172 {
173  reference_element::variant_type K_type = hat_K.variant();
174  _quad[K_type].initialize (hat_K, _options);
175  _initialized [K_type] = true;
176 }
177 template<class T>
178 void
179 quadrature_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 }
184 template<class T>
187 {
188  reference_element::variant_type K_type = hat_K.variant();
189  if (!_initialized [K_type]) _initialize (hat_K);
190  return _quad[K_type].begin();
191 }
192 template<class T>
195 {
196  reference_element::variant_type K_type = hat_K.variant();
197  if (!_initialized [K_type]) _initialize (hat_K);
198  return _quad[K_type].end();
199 }
200 template<class T>
203 {
204  reference_element::variant_type K_type = hat_K.variant();
205  if (!_initialized [K_type]) _initialize (hat_K);
206  return _quad[K_type].size();
207 }
208 template<class T>
209 const weighted_point<T>&
211 {
212  reference_element::variant_type K_type = hat_K.variant();
213  if (!_initialized [K_type]) _initialize (hat_K);
214  return _quad[K_type][q];
215 }
216 template<class T>
217 ostream&
218 operator<< (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 }
226 template<class T>
227 ostream&
228 operator<< (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 // ----------------------------------------------------------------------------
245 template <class T>
247 {
248  if (_options.name() != "") {
250  }
251 }
252 template<class T>
253 void
254 quadrature<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 }
262 template<class T>
263 const quadrature_option&
265 {
266  return base::data().get_options();
267 }
268 template<class T>
270  : base(),
272 {
273  reset (opt.name());
274 }
275 template<class T>
276 quadrature<T>::quadrature (const std::string& name)
277  : base(),
279 {
280  reset (name);
281 }
282 template<class T>
283 void
285 {
286  quadrature_option qopt = get_options();
287  qopt.set_order (order);
288  reset (qopt.name());
289 }
290 template<class T>
291 void
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) \
302 template class quadrature_on_geo<Float>; \
303 template class quadrature_rep<Float>; \
304 template class quadrature<T>; \
305 template ostream& operator<< (ostream&, const quadrature_on_geo<T>&); \
306 template ostream& operator<< (ostream&, const quadrature_rep<T>&);
307 
309 
310 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:425
see the Float 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
integrate_option quadrature_option
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:261
void reset(T &x)
This file is part of Rheolef.
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
void load(idiststream &in, Float &p, field &uh)