GetFEM  5.5
getfem_integration.cc
1 /*===========================================================================
2 
3  Copyright (C) 2000-2026 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program. If not, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
21 #include "getfem/bgeot_torus.h"
22 #include "getfem/dal_singleton.h"
23 #include "getfem/getfem_integration.h"
24 #include "gmm/gmm_dense_lu.h"
25 #include "getfem/bgeot_permutations.h"
26 #include "getfem/bgeot_geotrans_inv.h"
27 #include "getfem/getfem_im_list.h"
28 #include "getfem/dal_naming_system.h"
29 #include "getfem/getfem_omp.h"
30 #include "getfem/getfem_torus.h"
31 
32 namespace getfem {
33 
34  /*
35  * dummy integration method
36  */
37  static pintegration_method im_none(im_param_list &params,
38  std::vector<dal::pstatic_stored_object> &) {
39  GMM_ASSERT1(params.size() == 0, "IM_NONE does not accept any parameter");
40  return std::make_shared<integration_method>();
41  }
42 
43  long_scalar_type poly_integration::int_poly(const base_poly &P) const {
44  long_scalar_type res = 0.0;
45  if (P.size() > int_monomials.size()) {
46  std::vector<long_scalar_type> *hum = &int_monomials;
47  size_type i = P.size(), j = int_monomials.size();
48  hum->resize(i);
49  bgeot::power_index mi(P.dim()); mi[P.dim()-1] = P.degree();
50  for (size_type k = i; k > j; --k, --mi)
51  (*hum)[k-1] = int_monomial(mi);
52  }
53  base_poly::const_iterator it = P.begin(), ite = P.end();
54  std::vector<long_scalar_type>::const_iterator itb = int_monomials.begin();
55  for ( ; it != ite; ++it, ++itb) {
56  res += long_scalar_type(*it) * long_scalar_type(*itb);
57  }
58  return res;
59  }
60 
61  long_scalar_type
62  poly_integration::int_poly_on_face(const base_poly &P,short_type f) const {
63  long_scalar_type res = 0.0;
64  std::vector<long_scalar_type> *hum = &(int_face_monomials[f]);
65  if (P.size() > hum->size()) {
66  size_type i = P.size(), j = hum->size();
67  hum->resize(i);
68  bgeot::power_index mi(P.dim()); mi[P.dim()-1] = P.degree();
69  for (size_type k = i; k > j; --k, --mi)
70  (*hum)[k-1] = int_monomial_on_face(mi, f);
71  }
72  base_poly::const_iterator it = P.begin(), ite = P.end();
73  std::vector<long_scalar_type>::const_iterator itb = hum->begin();
74  for ( ; it != ite; ++it, ++itb)
75  res += long_scalar_type(*it) * long_scalar_type(*itb);
76  return res;
77  }
78 
79  /* ******************************************************************** */
80  /* integration on simplex */
81  /* ******************************************************************** */
82 
83  struct simplex_poly_integration_ : public poly_integration {
84 
85  long_scalar_type int_monomial(const bgeot::power_index &power) const;
86 
87  long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
88  short_type f) const;
89 
90  simplex_poly_integration_(bgeot::pconvex_structure c)
91  { cvs = c; int_face_monomials.resize(c->nb_faces()); }
92  };
93 
94 
95  long_scalar_type
96  simplex_poly_integration_::int_monomial(const bgeot::power_index &power) const{
97  long_scalar_type res = LONG_SCAL(1);
98  short_type fa = 1;
99  bgeot::power_index::const_iterator itm = power.begin(),
100  itme = power.end();
101  for ( ; itm != itme; ++itm)
102  for (int k = 1; k <= *itm; ++k, ++fa)
103  res *= long_scalar_type(k) / long_scalar_type(fa);
104 
105  for (int k = 0; k < cvs->dim(); k++) { res /= long_scalar_type(fa); fa++; }
106  return res;
107  }
108 
109  long_scalar_type simplex_poly_integration_::int_monomial_on_face
110  (const bgeot::power_index &power, short_type f) const {
111  long_scalar_type res = LONG_SCAL(0);
112 
113  if (f == 0 || power[f-1] == 0.0) {
114  res = (f == 0) ? sqrt(long_scalar_type(cvs->dim())) : LONG_SCAL(1);
115  short_type fa = 1;
116  bgeot::power_index::const_iterator itm = power.begin(),
117  itme = power.end();
118  for ( ; itm != itme; ++itm)
119  for (int k = 1; k <= *itm; ++k, ++fa)
120  res *= long_scalar_type(k) / long_scalar_type(fa);
121 
122  for (int k = 1; k < cvs->dim(); k++) { res/=long_scalar_type(fa); fa++; }
123  }
124  return res;
125  }
126 
127  static pintegration_method
128  exact_simplex(im_param_list &params,
129  std::vector<dal::pstatic_stored_object> &dependencies) {
130  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
131  << params.size() << " should be 1.");
132  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
133  int n = int(::floor(params[0].num() + 0.01));
134  GMM_ASSERT1(n > 0 && n < 100 && double(n) == params[0].num(),
135  "Bad parameters");
136  dependencies.push_back(bgeot::simplex_structure(dim_type(n)));
137  ppoly_integration ppi = std::make_shared<simplex_poly_integration_>
138  (bgeot::simplex_structure(dim_type(n)));
139  return std::make_shared<integration_method>(ppi);
140  }
141 
142  /* ******************************************************************** */
143  /* integration on direct product of convex structures */
144  /* ******************************************************************** */
145 
146  struct plyint_mul_structure_ : public poly_integration {
147  ppoly_integration cv1, cv2;
148 
149  long_scalar_type int_monomial(const bgeot::power_index &power) const;
150 
151  long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
152  short_type f) const;
153 
154  plyint_mul_structure_(ppoly_integration a, ppoly_integration b);
155  };
156 
157  long_scalar_type plyint_mul_structure_::int_monomial
158  (const bgeot::power_index &power) const {
159  bgeot::power_index mi1(cv1->dim()), mi2(cv2->dim());
160  std::copy(power.begin(), power.begin() + cv1->dim(), mi1.begin());
161  std::copy(power.begin() + cv1->dim(), power.end(), mi2.begin());
162  return cv1->int_monomial(mi1) * cv2->int_monomial(mi2);
163  }
164 
165  long_scalar_type plyint_mul_structure_::int_monomial_on_face
166  (const bgeot::power_index &power, short_type f) const {
167  bgeot::power_index mi1(cv1->dim()), mi2(cv2->dim());
168  std::copy(power.begin(), power.begin() + cv1->dim(), mi1.begin());
169  std::copy(power.begin() + cv1->dim(), power.end(), mi2.begin());
170  short_type nfx = cv1->structure()->nb_faces();
171  if (f < nfx)
172  return cv1->int_monomial_on_face(mi1,f) * cv2->int_monomial(mi2);
173  else
174  return cv1->int_monomial(mi1)
175  * cv2->int_monomial_on_face(mi2, short_type(f-nfx));
176  }
177 
178  plyint_mul_structure_::plyint_mul_structure_(ppoly_integration a,
179  ppoly_integration b) {
180  cv1 = a; cv2 = b;
181  cvs = bgeot::convex_product_structure(cv1->structure(),
182  cv2->structure());
183  int_face_monomials.resize(cvs->nb_faces());
184  }
185 
186  static pintegration_method
187  product_exact(im_param_list &params,
188  std::vector<dal::pstatic_stored_object> &dependencies) {
189  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
190  << params.size() << " should be 2.");
191  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
192  "Bad type of parameters");
193  pintegration_method a = params[0].method();
194  pintegration_method b = params[1].method();
195  GMM_ASSERT1(a->type() == IM_EXACT && b->type() == IM_EXACT,
196  "Bad parameters");
197  dependencies.push_back(a); dependencies.push_back(b);
198  dependencies.push_back(bgeot::convex_product_structure(a->structure(),
199  b->structure()));
200  ppoly_integration ppi = std::make_shared<plyint_mul_structure_>
201  (a->exact_method(), b->exact_method());
202  return std::make_shared<integration_method>(ppi);
203  }
204 
205  /* ******************************************************************** */
206  /* integration on parallelepiped. */
207  /* ******************************************************************** */
208 
209  static pintegration_method
210  exact_parallelepiped(im_param_list &params,
211  std::vector<dal::pstatic_stored_object> &) {
212  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
213  << params.size() << " should be 1.");
214  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
215  int n = int(::floor(params[0].num() + 0.01));
216  GMM_ASSERT1(n > 0 && n < 100 && double(n) == params[0].num(),
217  "Bad parameters");
218 
219  std::stringstream name;
220  if (n == 1)
221  name << "IM_EXACT_SIMPLEX(1)";
222  else
223  name << "IM_PRODUCT(IM_EXACT_PARALLELEPIPED(" << n-1
224  << "),IM_EXACT_SIMPLEX(1)))";
225  return int_method_descriptor(name.str());
226  }
227 
228  static pintegration_method
229  exact_prism(im_param_list &params,
230  std::vector<dal::pstatic_stored_object> &) {
231  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
232  << params.size() << " should be 1.");
233  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
234  int n = int(::floor(params[0].num() + 0.01));
235  GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
236  "Bad parameters");
237 
238  std::stringstream name;
239  name << "IM_PRODUCT(IM_EXACT_SIMPLEX(" << n-1
240  << "),IM_EXACT_SIMPLEX(1))";
241  return int_method_descriptor(name.str());
242  }
243 
244  /* ********************************************************************* */
245  /* Approximated integration */
246  /* ********************************************************************* */
247 
248  void approx_integration::add_point(const base_node &pt,
249  scalar_type w,
250  short_type f,
251  bool include_empty) {
252  GMM_ASSERT1(!valid, "Impossible to modify a valid integration method.");
253  if (gmm::abs(w) > 1.0E-15 || include_empty) {
254  ++f;
255  if (gmm::abs(w) <= 1.0E-15) w = scalar_type(0);
256  GMM_ASSERT1(f <= cvr->structure()->nb_faces(), "Wrong argument.");
257  size_type i = pt_to_store[f].search_node(pt);
258  if (i == size_type(-1)) {
259  i = pt_to_store[f].add_node(pt);
260  int_coeffs.resize(int_coeffs.size() + 1);
261  for (size_type j = f; j <= cvr->structure()->nb_faces(); ++j)
262  repartition[j] += 1;
263  for (size_type j = int_coeffs.size(); j > repartition[f]; --j)
264  int_coeffs[j-1] = int_coeffs[j-2];
265  int_coeffs[repartition[f]-1] = 0.0;
266  }
267  int_coeffs[((f == 0) ? 0 : repartition[f-1]) + i] += w;
268  }
269  }
270 
271  void approx_integration::add_point_norepeat(const base_node &pt,
272  scalar_type w,
273  short_type f) {
274  short_type f2 = f; ++f2;
275  if (pt_to_store[f2].search_node(pt) == size_type(-1)) add_point(pt,w,f);
276  }
277 
278  void approx_integration::add_point_full_symmetric(base_node pt,
279  scalar_type w) {
280  dim_type n = cvr->structure()->dim();
281  dim_type k;
282  base_node pt2(n);
283  if (n+1 == cvr->structure()->nb_faces()) {
284  // simplices
285  // for a point at (x,y) you have to consider the other points
286  // at (y,x) (x,1-x-y) (1-x-y,x) (y,1-x-y) (1-x-y,y)
287  base_node pt3(n+1);
288  std::copy(pt.begin(), pt.end(), pt3.begin());
289  pt3[n] = 1; for (k = 0; k < n; ++k) pt3[n] -= pt[k];
290  std::vector<int> ind(n); std::fill(ind.begin(), ind.end(), 0);
291  std::vector<bool> ind2(n+1);
292  for (;;) {
293 
294  std::fill(ind2.begin(), ind2.end(), false);
295  bool good = true;
296  for (k = 0; k < n; ++k)
297  if (ind2[ind[k]]) { good = false; break; } else ind2[ind[k]] = true;
298  if (good) {
299  for (k = 0; k < n; ++k) pt2[k] = pt3[ind[k]];
300  add_point_norepeat(pt2, w);
301  }
302  ind[0]++; k = 0;
303  while(ind[k] == n+1) { ind[k++] = 0; if (k == n) return; ind[k]++; }
304  }
305  }
306 
307  else if (cvr->structure()->nb_points() == (size_type(1) << n)) {
308  // parallelepidedic
309  for (size_type i = 0; i < (size_type(1) << n); ++i) {
310  for (k = 0; k < n; ++k)
311  if (i & (size_type(1) << k)) pt2[k]=pt[k]; else pt2[k] = 1.0-pt[k];
312  add_point_norepeat(pt2, w);
313  }
314  }
315  else
316  GMM_ASSERT1(false, "Fully symmetric option is only valid for"
317  "simplices and parallelepipedic elements");
318  }
319 
320  void approx_integration::add_method_on_face(pintegration_method ppi,
321  short_type f) {
322  papprox_integration pai = ppi->approx_method();
323  GMM_ASSERT1(!valid, "Impossible to modify a valid integration method.");
324  GMM_ASSERT1(*key_of_stored_object(pai->structure())
325  == *key_of_stored_object(structure()->faces_structure()[f]),
326  "structures missmatch");
327  GMM_ASSERT1(ppi->type() == IM_APPROX, "Impossible with an exact method.");
328 
329  dim_type N = pai->structure()->dim();
330  scalar_type det = 1.0;
331  base_node pt(N+1);
332  std::vector<base_node> pts(N);
333  for (size_type i = 0; i < N; ++i)
334  pts[i] = (cvr->dir_points_of_face(f))[i+1]
335  - (cvr->dir_points_of_face(f))[0];
336  if (N) {
337  base_matrix a(N+1, N), b(N, N), tmp(N, N);
338  for (dim_type i = 0; i < N+1; ++i)
339  for (dim_type j = 0; j < N; ++j)
340  a(i, j) = pts[j][i];
341 
342  gmm::mult(gmm::transposed(a), a, b);
343  det = ::sqrt(gmm::abs(gmm::lu_det(b)));
344  }
345  for (size_type i = 0; i < pai->nb_points_on_convex(); ++i) {
346  pt = (cvr->dir_points_of_face(f))[0];
347  for (dim_type j = 0; j < N; ++j)
348  pt += pts[j] * ((*(pai->pintegration_points()))[i])[j];
349  add_point(pt, pai->coeff(i) * det, f);
350  }
351  }
352 
353  void approx_integration::valid_method() {
354  std::vector<base_node> ptab(int_coeffs.size());
355  // std::vector<scalar_type> int_coeffs2(int_coeffs);
356  size_type i = 0;
357  for (short_type f = 0; f <= cvr->structure()->nb_faces(); ++f) {
358  // size_type j = i;
359  for (PT_TAB::const_iterator it = pt_to_store[f].begin();
360  it != pt_to_store[f].end(); ++it /* , ++j */) {
361  // int_coeffs[i] = int_coeffs2[j];
362  ptab[i++] = *it;
363  }
364  }
365  GMM_ASSERT1(i == int_coeffs.size(), "internal error.");
366  pint_points = bgeot::store_point_tab(ptab);
367  pt_to_store = std::vector<PT_TAB>();
368  valid = true;
369  }
370 
371 
372  /* ********************************************************************* */
373  /* methods stored in getfem_im_list.h */
374  /* ********************************************************************* */
375 
376  /// search a method in getfem_im_list.h
377  static pintegration_method
378  im_list_integration(std::string name,
379  std::vector<dal::pstatic_stored_object> &dependencies) {
380  // cerr << "searching " << name << endl;
381  for (int i = 0; i < NB_IM; ++i)
382  if (!name.compare(im_desc_tab[i].method_name)) {
383  bgeot::pgeometric_trans pgt
384  = bgeot::geometric_trans_descriptor(im_desc_tab[i].geotrans_name);
385  dim_type N = pgt->structure()->dim();
386  base_node pt(N);
387  auto pai = std::make_shared<approx_integration>(pgt->convex_ref());
388  size_type fr = im_desc_tab[i].firstreal;
389  for (size_type j = 0; j < im_desc_tab[i].nb_points; ++j) {
390  for (dim_type k = 0; k < N; ++k)
391  pt[k] = atof(im_desc_real[fr+j*(N+1)+k]);
392  // pt[k] = LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+k]);
393 
394  switch (im_desc_node_type[im_desc_tab[i].firsttype + j]) {
395  case 2: {
396  base_node pt2(pt.size());
397  for (bgeot::permutation p(pt.size()); !p.finished(); ++p) {
398  p.apply_to(pt,pt2);
399  pai->add_point_full_symmetric(pt2,
400  atof(im_desc_real[fr+j*(N+1)+N]));
401  // pai->add_point_full_symmetric(pt2,
402  // LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
403  }
404  } break;
405  case 1: {
406  pai->add_point_full_symmetric(pt,atof(im_desc_real[fr+j*(N+1)+N]));
407  // pai->add_point_full_symmetric(pt,
408  // LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
409  } break;
410  case 0: {
411  pai->add_point(pt, atof(im_desc_real[fr+j*(N+1)+N]));
412  // pai->add_point(pt,LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
413  } break;
414  default: GMM_ASSERT1(false, "");
415  }
416  }
417 
418  for (short_type f = 0; N > 0 && f < pgt->structure()->nb_faces(); ++f)
419  pai->add_method_on_face
420  (int_method_descriptor
421  (im_desc_face_meth[im_desc_tab[i].firstface + f]), f);
422 
423  pai->valid_method();
424  // cerr << "finding " << name << endl;
425 
426  pintegration_method p(std::make_shared<integration_method>(pai));
427  dependencies.push_back(p->approx_method()->ref_convex());
428  dependencies.push_back(p->approx_method()->pintegration_points());
429  return p;
430  }
431  return pintegration_method();
432  }
433 
434 
435  /* ********************************************************************* */
436  /* Gauss method. */
437  /* ********************************************************************* */
438 
439  struct Legendre_polynomials {
440  std::vector<base_poly> polynomials;
441  std::vector<std::vector<long_scalar_type>> roots;
442  int nb_lp;
443  Legendre_polynomials() : nb_lp(-1) {}
444  void init(short_type de) {
445  polynomials.resize(de+2);
446  roots.resize(de+2);
447  if (nb_lp < 0) {
448  polynomials[0] = base_poly(1,0);
449  polynomials[0].one();
450  polynomials[1] = base_poly(1,1,0);
451  roots[1].resize(1);
452  roots[1][0] = 0.0;
453  nb_lp = 1;
454  }
455  while (nb_lp < de) {
456  ++nb_lp;
457  polynomials[nb_lp] =
458  (base_poly(1,1,0) * polynomials[nb_lp-1]
459  * ((2.0 * long_scalar_type(nb_lp) - 1.0) / long_scalar_type(nb_lp)))
460  + (polynomials[nb_lp-2]
461  * ((1.0 - long_scalar_type(nb_lp)) / long_scalar_type(nb_lp)));
462  roots[nb_lp].resize(nb_lp);
463  roots[nb_lp][nb_lp/2] = 0.0;
464  long_scalar_type a = -1.0, b, c, d, e, cv, ev, ecart, ecart2;
465  for (int k = 0; k < nb_lp / 2; ++k) { // + symetry ...
466  b = roots[nb_lp-1][k];
467 
468  c = a, d = b;
469  cv = polynomials[nb_lp].eval(&c);
470  int imax = 10000;
471  ecart = 2.0 * (d - c);
472  while(c != d) {
473  --imax; if (imax == 0) break;
474  e = (c + d) / 2.0;
475  ecart2 = d - c; if (ecart2 >= ecart) break;
476  ecart = ecart2;
477  ev = polynomials[nb_lp].eval(&e);
478  if (ev * cv < 0.) { d = e; } else { c = e; cv = ev; }
479  }
480 
481  roots[nb_lp][k] = c;
482  roots[nb_lp][nb_lp-k-1] = -c;
483  a = b;
484  }
485  }
486  }
487  };
488 
489  struct gauss_approx_integration_ : public approx_integration {
490  gauss_approx_integration_(short_type nbpt);
491  };
492 
493  gauss_approx_integration_::gauss_approx_integration_(short_type nbpt) {
494  GMM_ASSERT1(nbpt <= 32000, "too much points");
495 
496  cvr = bgeot::simplex_of_reference(1);
497  std::vector<base_node> int_points(nbpt+2);
498  int_coeffs.resize(nbpt+2);
499  repartition.resize(3);
500  repartition[0] = nbpt;
501  repartition[1] = nbpt + 1;
502  repartition[2] = nbpt + 2;
503 
504  Legendre_polynomials Lp;
505  Lp.init(nbpt);
506 
507  for (short_type i = 0; i < nbpt; ++i) {
508  int_points[i].resize(1);
509  long_scalar_type lr = Lp.roots[nbpt][i];
510  int_points[i][0] = 0.5 + 0.5 * bgeot::to_scalar(lr);
511  int_coeffs[i] = bgeot::to_scalar((1.0 - gmm::sqr(lr))
512  / gmm::sqr( long_scalar_type(nbpt)
513  * (Lp.polynomials[nbpt-1].eval(&lr))));
514  }
515 
516  int_points[nbpt].resize(1);
517  int_points[nbpt][0] = 1.0; int_coeffs[nbpt] = 1.0;
518 
519  int_points[nbpt+1].resize(1);
520  int_points[nbpt+1][0] = 0.0; int_coeffs[nbpt+1] = 1.0;
521  pint_points = bgeot::store_point_tab(int_points);
522  valid = true;
523  }
524 
525 
526  static pintegration_method
527  gauss1d(im_param_list &params,
528  std::vector<dal::pstatic_stored_object> &dependencies) {
529  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
530  << params.size() << " should be 1.");
531  GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
532  int n = int(::floor(params[0].num() + 0.01));
533  GMM_ASSERT1(n >= 0 && n < 32000 && double(n) == params[0].num(),
534  "Bad parameters");
535  if (n & 1) {
536  std::stringstream name;
537  name << "IM_GAUSS1D(" << n-1 << ")";
538  return int_method_descriptor(name.str());
539  }
540  else {
541  papprox_integration
542  pai = std::make_shared<gauss_approx_integration_>(short_type(n/2 + 1));
543  pintegration_method p = std::make_shared<integration_method>(pai);
544  dependencies.push_back(p->approx_method()->ref_convex());
545  dependencies.push_back(p->approx_method()->pintegration_points());
546  return p;
547  }
548  }
549 
550  /* ********************************************************************* */
551  /* integration on simplexes */
552  /* ********************************************************************* */
553 
554  struct Newton_Cotes_approx_integration_ : public approx_integration
555  {
556  // void calc_base_func(base_poly &p, short_type K, base_node &c) const;
557  Newton_Cotes_approx_integration_(dim_type nc, short_type k);
558  };
559 
560  Newton_Cotes_approx_integration_::Newton_Cotes_approx_integration_
561  (dim_type nc, short_type k)
562  : approx_integration(bgeot::simplex_of_reference(nc)) {
563  size_type R = bgeot::alpha(nc,k);
564 
565  base_node c(nc);
566  if (nc == 0)
567  add_point(c, scalar_type(1));
568  else {
569  std::stringstream name;
570  name << "IM_EXACT_SIMPLEX(" << int(nc) << ")";
571  ppoly_integration ppi = int_method_descriptor(name.str())->exact_method();
572 
573  c.fill(scalar_type(0.0));
574  if (k == 0)
575  c.fill(1.0 / scalar_type(nc+1));
576 
577  gmm::dense_matrix<long_scalar_type> M(R, R);
578  std::vector<long_scalar_type> F(R), U(R);
579  std::vector<bgeot::power_index> base(R);
580  std::vector<base_node> nodes(R);
581 
582  bgeot::power_index pi(nc);
583 
584  size_type sum = 0;
585  if (k != 0 && nc > 0)
586  for (size_type r = 0; r < R; ++r, ++pi) {
587  base[r] = pi;
588  nodes[r] = c;
589  size_type l = 0;
590  c[l] += 1.0 / scalar_type(k);
591  sum++;
592  while (sum > k) {
593  sum -= int(floor(0.5+(c[l] * k)));
594  c[l++] = 0.0;
595  if (l == nc)
596  break;
597  c[l] += 1.0 / scalar_type(k);
598  sum++;
599  }
600  }
601  else // not sure if the following loop is really necessary
602  for (size_type r = 0; r < R; ++r, ++pi) {
603  base[r] = pi;
604  nodes[r] = c;
605  }
606 // if (nc == 1) {
607 // M = bgeot::vsmatrix<long_scalar_type>((R+1)/2, (R+1)/2);
608 // U = F = bgeot::vsvector<long_scalar_type>((R+1)/2);
609 // gmm::clear(M);
610 // }
611 
612  for (size_type r = 0; r < R; ++r) {
613 // if (nc == 1) {
614 // if (r < (R+1)/2) {
615 // F[r] = ppi->int_monomial(base[R-1-r]);
616 // cout << "F[" << r << "] = " << F[r] << endl;
617 // }
618 // }
619 // else {
620  F[r] = ppi->int_monomial(base[r]);
621  //cout << "F[" << r << "] = " << F[r] << endl;
622 // }
623  for (size_type q = 0; q < R; ++q) {
624 // if (nc == 1) {
625 // if (r < (R+1)/2) {
626 // if (q < (R+1)/2)
627 // M(r, q) += bgeot::eval_monomial(base[R-1-r], nodes[q].begin());
628 // else
629 // M(r, R-1-q) += bgeot::eval_monomial(base[R-1-r],
630 // nodes[q].begin());
631 // }
632 // }
633 // else
634  M(r, q) = bgeot::eval_monomial(base[r], nodes[q].begin());
635  }
636  }
637 
638  gmm::lu_solve(M, U, F);
639  for (size_type r = 0; r < R; ++r)
640  add_point(nodes[r], bgeot::to_scalar(U[r]));
641 
642  std::stringstream name2;
643  name2 << "IM_NC(" << int(nc-1) << "," << int(k) << ")";
644  for (short_type f = 0; f < structure()->nb_faces(); ++f)
645  add_method_on_face(int_method_descriptor(name2.str()), f);
646  }
647  valid_method();
648  }
649 
650  static pintegration_method
651  Newton_Cotes(im_param_list &params,
652  std::vector<dal::pstatic_stored_object> &dependencies) {
653  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
654  << params.size() << " should be 2.");
655  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
656  "Bad type of parameters");
657  int n = int(::floor(params[0].num() + 0.01));
658  int k = int(::floor(params[1].num() + 0.01));
659  GMM_ASSERT1(n >= 0 && n < 100 && k >= 0 && k <= 150 &&
660  double(n) == params[0].num() && double(k) == params[1].num(),
661  "Bad parameters");
662  papprox_integration
663  pai = std::make_shared<Newton_Cotes_approx_integration_>(dim_type(n),
664  short_type(k));
665  pintegration_method p = std::make_shared<integration_method>(pai);
666  dependencies.push_back(p->approx_method()->ref_convex());
667  dependencies.push_back(p->approx_method()->pintegration_points());
668  return p;
669  }
670 
671  /* ********************************************************************* */
672  /* integration on direct product of convex structures */
673  /* ********************************************************************* */
674 
675  struct a_int_pro_integration : public approx_integration
676  {
677  a_int_pro_integration(papprox_integration a, papprox_integration b);
678  };
679 
680 
681  a_int_pro_integration::a_int_pro_integration(papprox_integration a,
682  papprox_integration b) {
683  cvr = bgeot::convex_ref_product(a->ref_convex(), b->ref_convex());
684  size_type n1 = a->nb_points_on_convex();
685  size_type n2 = b->nb_points_on_convex();
686  std::vector<base_node> int_points;
687  int_points.resize(n1 * n2);
688  int_coeffs.resize(n1 * n2);
689  repartition.resize(cvr->structure()->nb_faces()+1);
690  repartition[0] = n1 * n2;
691  for (size_type i1 = 0; i1 < n1; ++i1)
692  for (size_type i2 = 0; i2 < n2; ++i2) {
693  int_coeffs[i1 + i2 * n1] = a->coeff(i1) * b->coeff(i2);
694  int_points[i1 + i2 * n1].resize(dim());
695  std::copy(a->point(i1).begin(), a->point(i1).end(),
696  int_points[i1 + i2 * n1].begin());
697  std::copy(b->point(i2).begin(), b->point(i2).end(),
698  int_points[i1 + i2 * n1].begin() + a->dim());
699  }
700  short_type f = 0;
701  for (short_type f1 = 0; f1 < a->structure()->nb_faces(); ++f1, ++f) {
702  n1 = a->nb_points_on_face(f1);
703  n2 = b->nb_points_on_convex();
704  size_type w = repartition[f];
705  repartition[f+1] = w + n1 * n2;
706  int_points.resize(repartition[f+1]);
707  int_coeffs.resize(repartition[f+1]);
708  for (size_type i1 = 0; i1 < n1; ++i1)
709  for (size_type i2 = 0; i2 < n2; ++i2) {
710  int_coeffs[w + i1 + i2 * n1] = a->coeff_on_face(f1, i1)
711  * b->coeff(i2);
712  int_points[w + i1 + i2 * n1].resize(dim());
713  std::copy(a->point_on_face(f1, i1).begin(),
714  a->point_on_face(f1, i1).end(),
715  int_points[w + i1 + i2 * n1].begin());
716  std::copy(b->point(i2).begin(), b->point(i2).end(),
717  int_points[w + i1 + i2 * n1].begin() + a->dim());
718  }
719  }
720  for (short_type f2 = 0; f2 < b->structure()->nb_faces(); ++f2, ++f) {
721  n1 = a->nb_points_on_convex();
722  n2 = b->nb_points_on_face(f2);
723  size_type w = repartition[f];
724  repartition[f+1] = w + n1 * n2;
725  int_points.resize(repartition[f+1]);
726  int_coeffs.resize(repartition[f+1]);
727  for (size_type i1 = 0; i1 < n1; ++i1)
728  for (size_type i2 = 0; i2 < n2; ++i2) {
729  int_coeffs[w + i1 + i2 * n1] = a->coeff(i1)
730  * b->coeff_on_face(f2, i2);
731  int_points[w + i1 + i2 * n1].resize(dim());
732  std::copy(a->point(i1).begin(), a->point(i1).end(),
733  int_points[w + i1 + i2 * n1].begin());
734  std::copy(b->point_on_face(f2, i2).begin(),
735  b->point_on_face(f2, i2).end(),
736  int_points[w + i1 + i2 * n1].begin() + a->dim());
737  }
738  }
739  pint_points = bgeot::store_point_tab(int_points);
740  valid = true;
741  }
742 
743  static pintegration_method
744  product_approx(im_param_list &params,
745  std::vector<dal::pstatic_stored_object> &dependencies) {
746  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
747  << params.size() << " should be 2.");
748  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
749  "Bad type of parameters");
750  pintegration_method a = params[0].method();
751  pintegration_method b = params[1].method();
752  GMM_ASSERT1(a->type() == IM_APPROX && b->type() == IM_APPROX,
753  "Bad parameters");
754  papprox_integration
755  pai = std::make_shared<a_int_pro_integration>(a->approx_method(),
756  b->approx_method());
757  pintegration_method p = std::make_shared<integration_method>(pai);
758  dependencies.push_back(p->approx_method()->ref_convex());
759  dependencies.push_back(p->approx_method()->pintegration_points());
760  return p;
761  }
762 
763  static pintegration_method
764  product_which(im_param_list &params,
765  std::vector<dal::pstatic_stored_object> &dependencies) {
766  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
767  << params.size() << " should be 2.");
768  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
769  "Bad type of parameters");
770  pintegration_method a = params[0].method();
771  pintegration_method b = params[1].method();
772  if (a->type() == IM_EXACT || b->type() == IM_EXACT)
773  return product_exact(params, dependencies);
774  else return product_approx(params, dependencies);
775  }
776 
777 
778  /* ********************************************************************* */
779  /* integration on parallelepiped with Newton Cotes formulae */
780  /* ********************************************************************* */
781 
782  static pintegration_method
783  Newton_Cotes_para(im_param_list &params,
784  std::vector<dal::pstatic_stored_object> &) {
785  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
786  << params.size() << " should be 2.");
787  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
788  "Bad type of parameters");
789  int n = int(::floor(params[0].num() + 0.01));
790  int k = int(::floor(params[1].num() + 0.01));
791  GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
792  double(n) == params[0].num() && double(k) == params[1].num(),
793  "Bad parameters");
794 
795  std::stringstream name;
796  if (n == 1)
797  name << "IM_NC(1," << k << ")";
798  else
799  name << "IM_PRODUCT(IM_NC_PARALLELEPIPED(" << n-1 << "," << k
800  << "),IM_NC(1," << k << "))";
801  return int_method_descriptor(name.str());
802  }
803 
804  static pintegration_method
805  Newton_Cotes_prism(im_param_list &params,
806  std::vector<dal::pstatic_stored_object> &) {
807  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
808  << params.size() << " should be 2.");
809  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
810  "Bad type of parameters");
811  int n = int(::floor(params[0].num() + 0.01));
812  int k = int(::floor(params[1].num() + 0.01));
813  GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
814  double(n) == params[0].num() && double(k) == params[1].num(),
815  "Bad parameters");
816 
817  std::stringstream name;
818  name << "IM_PRODUCT(IM_NC(" << n-1 << "," << k << "),IM_NC(1,"
819  << k << "))";
820  return int_method_descriptor(name.str());
821  }
822 
823  /* ********************************************************************* */
824  /* integration on parallelepiped with Gauss formulae */
825  /* ********************************************************************* */
826 
827  static pintegration_method
828  Gauss_paramul(im_param_list &params,
829  std::vector<dal::pstatic_stored_object> &) {
830  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
831  << params.size() << " should be 2.");
832  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
833  "Bad type of parameters");
834  int n = int(::floor(params[0].num() + 0.01));
835  int k = int(::floor(params[1].num() + 0.01));
836  GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
837  double(n) == params[0].num() && double(k) == params[1].num(),
838  "Bad parameters");
839 
840  std::stringstream name;
841  if (n == 1)
842  name << "IM_GAUSS1D(" << k << ")";
843  else
844  name << "IM_PRODUCT(IM_GAUSS_PARALLELEPIPED(" << n-1 << "," << k
845  << "),IM_GAUSS1D(" << k << "))";
846  return int_method_descriptor(name.str());
847  }
848 
849  /* ******************************************************************** */
850  /* Quasi-polar integration */
851  /* ******************************************************************** */
852 
853  struct quasi_polar_integration : public approx_integration {
854  quasi_polar_integration(papprox_integration base_im,
855  size_type ip1, size_type ip2=size_type(-1)) :
856  approx_integration
857  ((base_im->structure() == bgeot::parallelepiped_structure(3)) ?
858  bgeot::pyramid_QK_of_reference(1)
859  : bgeot::simplex_of_reference(base_im->dim())) {
860  size_type N = base_im->dim();
861 
862  enum { SQUARE, PRISM, TETRA_CYL, PRISM2, PYRAMID } what;
863  if (N == 2) what = SQUARE;
864  else if (base_im->structure() == bgeot::prism_P1_structure(3))
865  what = (ip2 == size_type(-1) || ip1 == ip2) ? PRISM2 : PRISM;
866  else if (base_im->structure() == bgeot::simplex_structure(3))
867  what = TETRA_CYL;
868  else if (base_im->structure() == bgeot::parallelepiped_structure(3))
869  what = PYRAMID;
870  else GMM_ASSERT1(false, "Incoherent integration method");
871 
872  // The first geometric transformation collapse a face of
873  // a parallelepiped or collapse a parrallelepiped on a pyramid.
874  // The second geometric transformation chooses the orientation.
875  // The third is used for the TETRA_CYL case only.
876  bgeot::pgeometric_trans pgt1 = bgeot::parallelepiped_geotrans(N,1);
877  std::vector<base_node> nodes1 = pgt1->convex_ref()->points();
878  bgeot::pgeometric_trans pgt2 = bgeot::simplex_geotrans(N, 1);
879  std::vector<base_node> nodes2(N+1);
880  if (what == PYRAMID) {
881  pgt2 = bgeot::pyramid_QK_geotrans(1);
882  nodes2.resize(5);
883  }
884  std::vector<size_type> other_nodes; // for the construction of node2
885  bgeot::pgeometric_trans pgt3 = bgeot::simplex_geotrans(N, 1);
886  std::vector<base_node> nodes3(N+1);
887 
888  switch (what) {
889  case SQUARE :
890  nodes1[3] = nodes1[1];
891  nodes2[ip1] = nodes1[1]; ip2 = ip1;
892  other_nodes.push_back(0);
893  other_nodes.push_back(2);
894  break;
895  case PRISM :
896  nodes1[4] = nodes1[0]; nodes1[5] = nodes1[1];
897  nodes2[ip1] = nodes1[0];
898  nodes2[ip2] = nodes1[1];
899  other_nodes.push_back(2);
900  other_nodes.push_back(6);
901  break;
902  case TETRA_CYL :
903  nodes3[0] = nodes1[0]; nodes3[1] = nodes1[1];
904  nodes3[2] = nodes1[2]; nodes3[3] = nodes1[4];
905  // nodes1[4] = nodes1[0]; nodes1[7] = base_node(1.0, 1.0, 2.0);
906  nodes2[ip1] = nodes1[1]; ip2 = ip1;
907  other_nodes.push_back(0);
908  other_nodes.push_back(2);
909  other_nodes.push_back(4);
910  break;
911  case PRISM2 :
912  nodes2[ip1] = nodes1[4];
913  other_nodes.push_back(0);
914  other_nodes.push_back(1);
915  other_nodes.push_back(2);
916  break;
917  case PYRAMID :
918  ip2 = ip1 = 0;
919  nodes1[0] = base_node(-1.,-1., 0.);
920  nodes1[1] = base_node( 1.,-1., 0.);
921  nodes1[2] = base_node(-1., 1., 0.);
922  nodes1[3] = base_node( 1., 1., 0.);
923  nodes1[4] = base_node( 0., 0., 1.);
924  nodes1[5] = nodes1[6] = nodes1[7] = nodes1[4];
925  nodes2[ip1] = nodes1[0];
926  other_nodes.push_back(4);
927  other_nodes.push_back(3);
928  other_nodes.push_back(2);
929  other_nodes.push_back(1);
930  }
931 
932  for (size_type i = 0; i <= nodes2.size()-1; ++i)
933  if (i != ip1 && i != ip2) {
934  GMM_ASSERT3(!other_nodes.empty(), "");
935  nodes2[i] = nodes1[other_nodes.back()];
936  other_nodes.pop_back();
937  }
938 
939  base_matrix G1, G2, G3;
940  base_matrix K(N, N), grad(N, N), K3(N, N), K4(N, N);
941  base_node normal1(N), normal2(N);
942  bgeot::geotrans_inv_convex gic(nodes2, pgt2);
943  scalar_type J1, J2, J3(1), J4(1);
944 
945  bgeot::vectors_to_base_matrix(G1, nodes1);
946  bgeot::vectors_to_base_matrix(G2, nodes2);
947 
948  for (size_type nc = 0; nc < 2; ++nc) {
949 
950  if (what == TETRA_CYL) {
951  if (nc == 1) nodes3[0] = nodes1[3];
952  bgeot::vectors_to_base_matrix(G3, nodes3);
953  pgt3->poly_vector_grad(nodes1[0], grad);
954  gmm::mult(gmm::transposed(grad), gmm::transposed(G3), K3);
955  J3 = gmm::abs(gmm::lu_inverse(K3)); /* = 1 */
956  }
957 
958  for (size_type i=0; i < base_im->nb_points(); ++i) {
959 
960  gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
961 
962  size_type fp = size_type(-1); // Search the face number in the
963  if (i >= base_im->nb_points_on_convex()) { // original element
964  size_type ii = i - base_im->nb_points_on_convex();
965  for (short_type ff=0; ff < base_im->structure()->nb_faces(); ++ff) {
966  if (ii < base_im->nb_points_on_face(ff)) { fp = ff; break; }
967  else ii -= base_im->nb_points_on_face(ff);
968  }
969  normal1 = base_im->ref_convex()->normals()[fp];
970  }
971 
972  base_node P = base_im->point(i);
973  if (what == TETRA_CYL) {
974  P = pgt3->transform(P, nodes3);
975  scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
976  K4(0, 1) = - y / one_minus_z;
977  K4(1, 1) = 1.0 - x / one_minus_z;
978  K4(2, 1) = - x * y / gmm::sqr(one_minus_z);
979  J4 = gmm::abs(gmm::lu_det(K4));
980  P[1] *= (1.0 - x / one_minus_z);
981  }
982  if (what == PRISM2) {
983  scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
984  K4(0,0) = one_minus_z; K4(2,0) = -x;
985  K4(1,1) = one_minus_z; K4(2,1) = -y;
986  J4 = gmm::abs(gmm::lu_det(K4));
987  P[0] *= one_minus_z;
988  P[1] *= one_minus_z;
989  }
990 
991  base_node P1 = pgt1->transform(P, nodes1), P2(N);
992  pgt1->poly_vector_grad(P, grad);
993  gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
994  J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
995 
996  if (fp != size_type(-1) && J1 > 1E-10 && J4 > 1E-10) {
997  if (what == TETRA_CYL) {
998  gmm::mult(K3, normal1, normal2);
999  normal1 = normal2;
1000  }
1001  gmm::lu_inverse(K4);
1002  gmm::lu_inverse(K);
1003  gmm::mult(K4, normal1, normal2);
1004  gmm::mult(K, normal2, normal1);
1005  normal2 = normal1;
1006  J1 *= gmm::vect_norm2(normal2);
1007  normal2 /= gmm::vect_norm2(normal2);
1008  }
1009  gic.invert(P1, P2);
1010  GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
1011  "Point not in the convex ref : " << P2);
1012 
1013  pgt2->poly_vector_grad(P2, grad);
1014  gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
1015  J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
1016 
1017  if (i < base_im->nb_points_on_convex())
1018  add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
1019  else if (J1 > 1E-10) {
1020  short_type f = short_type(-1);
1021  for (short_type ff = 0; ff <= N; ++ff)
1022  if (gmm::abs(pgt2->convex_ref()->is_in_face(ff, P2)) < 1E-8) {
1023  GMM_ASSERT1(f == short_type(-1),
1024  "An integration point is common to two faces");
1025  f = ff;
1026  }
1027  if (f != short_type(-1)) {
1028  gmm::mult(K, normal2, normal1);
1029  add_point(P2,base_im->coeff(i)*J1*gmm::vect_norm2(normal1)/J2,f);
1030  }
1031  // else { cout << "Point " << P2 << " eliminated" << endl; }
1032  }
1033  }
1034  if (what != TETRA_CYL) break;
1035  }
1036  valid_method();
1037  }
1038  };
1039 
1040 
1041  static pintegration_method
1042  quasi_polar(im_param_list &params,
1043  std::vector<dal::pstatic_stored_object> &dependencies) {
1044  GMM_ASSERT1(params.size() >= 1 && params[0].type() == 1,
1045  "Bad parameters for quasi polar integration: the first "
1046  "parameter should be an integration method");
1047  pintegration_method a = params[0].method();
1048  GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
1049  int ip1 = 0, ip2 = 0;
1050  if (a->approx_method()->structure() == bgeot::parallelepiped_structure(3)) {
1051  GMM_ASSERT1(params.size() == 1, "Bad number of parameters");
1052  } else {
1053  GMM_ASSERT1(params.size() == 2 || params.size() == 3,
1054  "Bad number of parameters : " << params.size()
1055  << " should be 2 or 3.");
1056  GMM_ASSERT1(params[1].type() == 0
1057  && params.back().type() == 0, "Bad type of parameters");
1058  ip1 = int(::floor(params[1].num() + 0.01));
1059  ip2 = int(::floor(params.back().num() + 0.01));
1060  }
1061  int N = a->approx_method()->dim();
1062  GMM_ASSERT1(N >= 2 && N <= 3 && ip1 >= 0 && ip2 >= 0 && ip1 <= N
1063  && ip2 <= N, "Bad parameters");
1064 
1065  papprox_integration
1066  pai = std::make_shared<quasi_polar_integration>(a->approx_method(),
1067  ip1, ip2);
1068  pintegration_method p = std::make_shared<integration_method>(pai);
1069  dependencies.push_back(p->approx_method()->ref_convex());
1070  dependencies.push_back(p->approx_method()->pintegration_points());
1071  return p;
1072  }
1073 
1074  static pintegration_method
1075  pyramid(im_param_list &params,
1076  std::vector<dal::pstatic_stored_object> &dependencies) {
1077  GMM_ASSERT1(params.size() == 1 && params[0].type() == 1,
1078  "Bad parameters for pyramid integration: the first "
1079  "parameter should be an integration method");
1080  pintegration_method a = params[0].method();
1081  GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
1082  int N = a->approx_method()->dim();
1083  GMM_ASSERT1(N == 3, "Bad parameters");
1084 
1085  papprox_integration
1086  pai = std::make_shared<quasi_polar_integration>(a->approx_method(), 0, 0);
1087  pintegration_method p = std::make_shared<integration_method>(pai);
1088  dependencies.push_back(p->approx_method()->ref_convex());
1089  dependencies.push_back(p->approx_method()->pintegration_points());
1090  return p;
1091  }
1092 
1093 
1094 
1095  /* ******************************************************************** */
1096  /* Naming system */
1097  /* ******************************************************************** */
1098 
1099  pintegration_method
1100  structured_composite_int_method(im_param_list &,
1101  std::vector<dal::pstatic_stored_object> &);
1102  pintegration_method HCT_composite_int_method(im_param_list &params,
1103  std::vector<dal::pstatic_stored_object> &dependencies);
1104 
1105  pintegration_method QUADC1_composite_int_method(im_param_list &params,
1106  std::vector<dal::pstatic_stored_object> &dependencies);
1107 
1108  pintegration_method pyramid_composite_int_method(im_param_list &params,
1109  std::vector<dal::pstatic_stored_object> &dependencies);
1110 
1111  struct im_naming_system : public dal::naming_system<integration_method> {
1112  im_naming_system() : dal::naming_system<integration_method>("IM") {
1113  add_suffix("NONE",im_none);
1114  add_suffix("EXACT_SIMPLEX", exact_simplex);
1115  add_suffix("PRODUCT", product_which);
1116  add_suffix("EXACT_PARALLELEPIPED",exact_parallelepiped);
1117  add_suffix("EXACT_PRISM", exact_prism);
1118  add_suffix("GAUSS1D", gauss1d);
1119  add_suffix("NC", Newton_Cotes);
1120  add_suffix("NC_PARALLELEPIPED", Newton_Cotes_para);
1121  add_suffix("NC_PRISM", Newton_Cotes_prism);
1122  add_suffix("GAUSS_PARALLELEPIPED", Gauss_paramul);
1123  add_suffix("QUASI_POLAR", quasi_polar);
1124  add_suffix("PYRAMID", pyramid);
1125  add_suffix("STRUCTURED_COMPOSITE",
1126  structured_composite_int_method);
1127  add_suffix("HCT_COMPOSITE",
1128  HCT_composite_int_method);
1129  add_suffix("QUADC1_COMPOSITE",
1130  QUADC1_composite_int_method);
1131  add_suffix("PYRAMID_COMPOSITE",
1132  pyramid_composite_int_method);
1133  add_generic_function(im_list_integration);
1134  }
1135  };
1136 
1137  pintegration_method int_method_descriptor(std::string name,
1138  bool throw_if_not_found) {
1139  size_type i = 0;
1140  return dal::singleton<im_naming_system>::instance().method
1141  (name, i, throw_if_not_found);
1142  }
1143 
1144  std::string name_of_int_method(pintegration_method p) {
1145  if (!(p.get())) return "IM_NONE";
1146  return dal::singleton<im_naming_system>::instance().shorter_name_of_method(p);
1147  }
1148 
1149  // allows the add of an integration method.
1150  void add_integration_name(std::string name,
1151  dal::naming_system<integration_method>::pfunction f) {
1152  dal::singleton<im_naming_system>::instance().add_suffix(name, f);
1153  }
1154 
1155 
1156  /* Fonctions pour la ref. directe. */
1157 
1158  pintegration_method exact_simplex_im(size_type n) {
1159  THREAD_SAFE_STATIC pintegration_method pim = nullptr;
1160  THREAD_SAFE_STATIC size_type d = -2;
1161  if (d != n) {
1162  std::stringstream name;
1163  name << "IM_EXACT_SIMPLEX(" << n << ")";
1164  pim = int_method_descriptor(name.str());
1165  d = n;
1166  }
1167  return pim;
1168  }
1169 
1170  pintegration_method exact_parallelepiped_im(size_type n) {
1171  THREAD_SAFE_STATIC pintegration_method pim = nullptr;
1172  THREAD_SAFE_STATIC size_type d = -2;
1173  if (d != n) {
1174  std::stringstream name;
1175  name << "IM_EXACT_PARALLELEPIPED(" << n << ")";
1176  pim = int_method_descriptor(name.str());
1177  d = n;
1178  }
1179  return pim;
1180  }
1181 
1182  pintegration_method exact_prism_im(size_type n) {
1183  THREAD_SAFE_STATIC pintegration_method pim = nullptr;
1184  THREAD_SAFE_STATIC size_type d = -2;
1185  if (d != n) {
1186  std::stringstream name;
1187  name << "IM_EXACT_PRISM(" << n << ")";
1188  pim = int_method_descriptor(name.str());
1189  d = n;
1190  }
1191  return pim;
1192  }
1193 
1194  pintegration_method exact_classical_im(bgeot::pgeometric_trans pgt) {
1195  return classical_exact_im(pgt);
1196  }
1197 
1198  static pintegration_method classical_exact_im(bgeot::pconvex_structure cvs) {
1199  cvs = bgeot::basic_structure(cvs);
1200  THREAD_SAFE_STATIC bgeot::pconvex_structure cvs_last = nullptr;
1201  THREAD_SAFE_STATIC pintegration_method im_last = nullptr;
1202  bool found = false;
1203 
1204  if (cvs_last == cvs)
1205  return im_last;
1206 
1207  size_type n = cvs->dim();
1208  size_type nbp = cvs->nb_points();
1209  std::stringstream name;
1210 
1211  /* Identifying P1-simplexes. */
1212 
1213  if (nbp == n+1)
1214  if (cvs == bgeot::simplex_structure(dim_type(n)))
1215  { name << "IM_EXACT_SIMPLEX("; found = true; }
1216 
1217  /* Identifying Q1-parallelepiped. */
1218 
1219  if (!found && nbp == (size_type(1) << n))
1220  if (cvs == bgeot::parallelepiped_structure(dim_type(n)))
1221  { name << "IM_EXACT_PARALLELEPIPED("; found = true; }
1222 
1223  /* Identifying Q1-prisms. */
1224 
1225  if (!found && nbp == 2 * n)
1226  if (cvs == bgeot::prism_P1_structure(dim_type(n)))
1227  { name << "IM_EXACT_PRISM("; found = true; }
1228 
1229  // To be completed
1230 
1231  if (found) {
1232  name << int(n) << ')';
1233  im_last = int_method_descriptor(name.str());
1234  cvs_last = cvs;
1235  return im_last;
1236  }
1237 
1238  GMM_ASSERT1(false, "This element is not taken into account. Contact us");
1239  }
1240 
1241 
1242  pintegration_method classical_exact_im(bgeot::pgeometric_trans pgt) {
1243  return classical_exact_im(pgt->structure());
1244  }
1245 
1246  static pintegration_method
1247  classical_approx_im_(bgeot::pconvex_structure cvs, dim_type degree) {
1248  size_type n = cvs->dim();
1249  std::stringstream name;
1250 
1251  if(bgeot::is_torus_structure(cvs) && n == 3) n = 2;
1252 
1253  degree = std::max<dim_type>(degree, 1);
1254  bgeot::pconvex_structure a, b;
1255  if (bgeot::basic_structure(cvs) == bgeot::simplex_structure(dim_type(n))) {
1256  /* Identifying P1-simplexes. */
1257  switch (n) {
1258  case 0: return int_method_descriptor("IM_NC(0,0)");
1259  case 1: name << "IM_GAUSS1D"; break;
1260  case 2: name << "IM_TRIANGLE"; break;
1261  case 3: name << "IM_TETRAHEDRON"; break;
1262  case 4: name << "IM_SIMPLEX4D"; break;
1263  default: GMM_ASSERT1(false, "no approximate integration method "
1264  "for simplexes of dimension " << n);
1265  }
1266  for (size_type k = degree; k < size_type(degree+10); ++k) {
1267  std::stringstream name2;
1268  name2 << name.str() << "(" << k << ")";
1269  pintegration_method im = int_method_descriptor(name2.str(), false);
1270  if (im) return im;
1271  }
1272  GMM_ASSERT1(false, "could not find an " << name.str()
1273  << " of degree >= " << int(degree));
1274  } else if (bgeot::basic_structure(cvs) == bgeot::pyramid_QK_structure(1)) {
1275  GMM_ASSERT1(n == 3, "Wrong dimension");
1276  name << "IM_PYRAMID(IM_GAUSS_PARALLELEPIPED(3," << degree << "))";
1277  } else if (cvs->is_product(&a,&b) ||
1278  (bgeot::basic_structure(cvs).get() &&
1279  bgeot::basic_structure(cvs)->is_product(&a,&b))) {
1280  name << "IM_PRODUCT("
1281  << name_of_int_method(classical_approx_im_(a,degree)) << ","
1282  << name_of_int_method(classical_approx_im_(b,degree)) << ")";
1283  } else GMM_ASSERT1(false, "unknown convex structure!");
1284  return int_method_descriptor(name.str());
1285  }
1286 
1287  pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt,
1288  dim_type degree) {
1289  THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_last = nullptr;
1290  THREAD_SAFE_STATIC dim_type degree_last;
1291  THREAD_SAFE_STATIC pintegration_method im_last = nullptr;
1292  if (pgt_last == pgt && degree == degree_last)
1293  return im_last;
1294  im_last = classical_approx_im_(pgt->structure(),degree);
1295  degree_last = degree;
1296  pgt_last = pgt;
1297  return im_last;
1298  }
1299 
1300  pintegration_method im_none() {
1301  THREAD_SAFE_STATIC pintegration_method im_last = nullptr;
1302  if (!im_last) im_last = int_method_descriptor("IM_NONE");
1303  return im_last;
1304  }
1305 
1306  /* try to integrate all monomials up to order 'order' and return the
1307  max. error */
1308  scalar_type test_integration_error(papprox_integration pim, dim_type order) {
1309  short_type dim = pim->dim();
1310  pintegration_method exact = classical_exact_im(pim->structure());
1311  opt_long_scalar_type error(0);
1312  for (bgeot::power_index idx(dim); idx.degree() <= order; ++idx) {
1313  opt_long_scalar_type sum(0), realsum;
1314  for (size_type i=0; i < pim->nb_points_on_convex(); ++i) {
1315  opt_long_scalar_type prod = pim->coeff(i);
1316  for (size_type d=0; d < dim; ++d)
1317  prod *= pow(opt_long_scalar_type(pim->point(i)[d]), idx[d]);
1318  sum += prod;
1319  }
1320  realsum = exact->exact_method()->int_monomial(idx);
1321  error = std::max(error, gmm::abs(realsum-sum));
1322  }
1323  return bgeot::to_scalar(error);
1324  }
1325 
1326  papprox_integration get_approx_im_or_fail(pintegration_method pim) {
1327  GMM_ASSERT1(pim->type() == IM_APPROX, "error estimate work only with "
1328  "approximate integration methods");
1329  return pim->approx_method();
1330  }
1331 
1332 } /* end of namespace getfem. */
bgeot_geotrans_inv.h
Inversion of geometric transformations.
bgeot_torus.h
Provides mesh of torus.
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:63
bgeot::permutation
generation of permutations, and ranking/unranking of these.
Definition: bgeot_permutations.h:43
bgeot::power_index
Vector of integer (16 bits type) which represent the powers of a monomial.
Definition: bgeot_poly.h:63
bgeot::power_index::degree
short_type degree() const
Gives the degree.
Definition: bgeot_poly.cc:89
dal::naming_system
Associate a name to a method descriptor and store method descriptors.
Definition: dal_naming_system.h:55
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:164
getfem::poly_integration
Description of an exact integration of polynomials.
Definition: getfem_integration.h:111
getfem::poly_integration::int_poly
long_scalar_type int_poly(const base_poly &P) const
Evaluate the integral of the polynomial P on the reference element.
Definition: getfem_integration.cc:43
getfem::poly_integration::structure
bgeot::pconvex_structure structure(void) const
{Structure of convex of reference.
Definition: getfem_integration.h:124
getfem::poly_integration::int_poly_on_face
long_scalar_type int_poly_on_face(const base_poly &P, short_type f) const
Evaluate the integral of the polynomial P on the face f of the reference element.
Definition: getfem_integration.cc:62
dal_naming_system.h
Naming system.
dal_singleton.h
A simple singleton implementation.
getfem_im_list.h
This file is generated by cubature/make_getfem_list.
getfem_integration.h
Integration methods (exact and approximated) on convexes.
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
getfem_torus.h
Provides mesh and mesh fem of torus.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:556
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm_dense_lu.h
LU factorizations and determinant computation for dense matrices.
gmm::lu_inverse
void lu_inverse(const DenseMatrixLU &LU, const Pvector &pvector, const DenseMatrix &AInv_)
Given an LU factored matrix, build the inverse of the matrix.
Definition: gmm_dense_lu.h:211
gmm::lu_solve
void lu_solve(const DenseMatrix &LU, const Pvector &pvector, VectorX &x, const VectorB &b)
LU Solve : Solve equation Ax=b, given an LU factored matrix.
Definition: gmm_dense_lu.h:129
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:26
bgeot::convex_ref_product
pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b)
tensorial product of two convex ref.
Definition: bgeot_convex_ref.cc:719
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::pyramid_QK_structure
pconvex_structure pyramid_QK_structure(dim_type k)
Give a pointer on the 3D pyramid structure for a degree k = 1 or 2.
Definition: bgeot_convex_structure.cc:513
bgeot::pyramid_QK_of_reference
pconvex_ref pyramid_QK_of_reference(dim_type k)
pyramidal element of reference of degree k (k = 1 or 2 only)
Definition: bgeot_convex_ref.cc:512
bgeot::simplex_of_reference
pconvex_ref simplex_of_reference(dim_type nc, short_type K)
returns a simplex of reference of dimension nc and degree k
Definition: bgeot_convex_ref.cc:330
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:53
bgeot::parallelepiped_structure
pconvex_structure parallelepiped_structure(dim_type nc, dim_type k)
Give a pointer on the structures of a parallelepiped of dimension d.
Definition: bgeot_convex_structure.cc:402
bgeot::prism_P1_structure
pconvex_structure prism_P1_structure(dim_type nc)
Give a pointer on the structures of a prism of dimension d.
Definition: bgeot_convex_structure.h:205
bgeot::simplex_structure
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
Definition: bgeot_convex_structure.cc:147
bgeot::geometric_trans_descriptor
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Definition: bgeot_geometric_trans.cc:1162
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::convex_product_structure
pconvex_structure convex_product_structure(pconvex_structure a, pconvex_structure b)
Give a pointer on the structures of a convex which is the direct product of the convexes represented ...
Definition: bgeot_convex_structure.cc:373
bgeot::basic_structure
pconvex_structure basic_structure(pconvex_structure cv)
Original structure (if concerned)
Definition: bgeot_convex_structure.h:171
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:185
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:46
dal
Dynamic Array Library.
Definition: dal_backtrace.cc:28
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::exact_parallelepiped_im
pintegration_method exact_parallelepiped_im(size_type n)
return IM_EXACT_PARALLELEPIPED(n)
Definition: getfem_integration.cc:1170
getfem::name_of_int_method
std::string name_of_int_method(pintegration_method p)
Get the string name of an integration method .
Definition: getfem_integration.cc:1144
getfem::exact_simplex_im
pintegration_method exact_simplex_im(size_type n)
return IM_EXACT_SIMPLEX(n)
Definition: getfem_integration.cc:1158
getfem::exact_prism_im
pintegration_method exact_prism_im(size_type n)
return IM_EXACT_PRISM(n)
Definition: getfem_integration.cc:1182
getfem::classical_approx_im
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree)
try to find an approximate integration method for the geometric transformation pgt which is able to i...
Definition: getfem_integration.cc:1287
getfem::classical_exact_im
pintegration_method classical_exact_im(bgeot::pgeometric_trans pgt)
return an exact integration method for convex type handled by pgt.
Definition: getfem_integration.cc:1242
getfem::exact_classical_im
pintegration_method exact_classical_im(bgeot::pgeometric_trans pgt) IS_DEPRECATED
use classical_exact_im instead.
Definition: getfem_integration.cc:1194
getfem::int_method_descriptor
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Definition: getfem_integration.cc:1137
getfem::im_none
pintegration_method im_none(void)
return IM_NONE
Definition: getfem_integration.cc:1300