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