GetFEM  5.4.3
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  }
570  else {
571 
572  std::stringstream name;
573  name << "IM_EXACT_SIMPLEX(" << int(nc) << ")";
574  ppoly_integration ppi = int_method_descriptor(name.str())->exact_method();
575 
576  size_type sum = 0, l;
577  c.fill(scalar_type(0.0));
578  if (k == 0) c.fill(1.0 / scalar_type(nc+1));
579 
580  gmm::dense_matrix<long_scalar_type> M(R, R);
581  std::vector<long_scalar_type> F(R), U(R);
582  std::vector<bgeot::power_index> base(R);
583  std::vector<base_node> nodes(R);
584 
585  bgeot::power_index pi(nc);
586 
587  for (size_type r = 0; r < R; ++r, ++pi) {
588  base[r] = pi; nodes[r] = c;
589  if (k != 0 && nc > 0) {
590  l = 0; c[l] += 1.0 / scalar_type(k); sum++;
591  while (sum > k) {
592  sum -= int(floor(0.5+(c[l] * k)));
593  c[l] = 0.0; l++; if (l == nc) break;
594  c[l] += 1.0 / scalar_type(k); sum++;
595  }
596  }
597  }
598 
599 // if (nc == 1) {
600 // M = bgeot::vsmatrix<long_scalar_type>((R+1)/2, (R+1)/2);
601 // U = F = bgeot::vsvector<long_scalar_type>((R+1)/2);
602 // gmm::clear(M);
603 // }
604 
605  for (size_type r = 0; r < R; ++r) {
606 // if (nc == 1) {
607 // if (r < (R+1)/2) {
608 // F[r] = ppi->int_monomial(base[R-1-r]);
609 // cout << "F[" << r << "] = " << F[r] << endl;
610 // }
611 // }
612 // else {
613  F[r] = ppi->int_monomial(base[r]);
614  //cout << "F[" << r << "] = " << F[r] << endl;
615 // }
616  for (size_type q = 0; q < R; ++q) {
617 // if (nc == 1) {
618 // if (r < (R+1)/2) {
619 // if (q < (R+1)/2)
620 // M(r, q) += bgeot::eval_monomial(base[R-1-r], nodes[q].begin());
621 // else
622 // M(r, R-1-q) += bgeot::eval_monomial(base[R-1-r],
623 // nodes[q].begin());
624 // }
625 // }
626 // else
627  M(r, q) = bgeot::eval_monomial(base[r], nodes[q].begin());
628  }
629  }
630 
631  gmm::lu_solve(M, U, F);
632  for (size_type r = 0; r < R; ++r)
633  add_point(nodes[r], bgeot::to_scalar(U[r]));
634 
635  std::stringstream name2;
636  name2 << "IM_NC(" << int(nc-1) << "," << int(k) << ")";
637  for (short_type f = 0; f < structure()->nb_faces(); ++f)
638  add_method_on_face(int_method_descriptor(name2.str()), f);
639  }
640  valid_method();
641  }
642 
643  static pintegration_method
644  Newton_Cotes(im_param_list &params,
645  std::vector<dal::pstatic_stored_object> &dependencies) {
646  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
647  << params.size() << " should be 2.");
648  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
649  "Bad type of parameters");
650  int n = int(::floor(params[0].num() + 0.01));
651  int k = int(::floor(params[1].num() + 0.01));
652  GMM_ASSERT1(n >= 0 && n < 100 && k >= 0 && k <= 150 &&
653  double(n) == params[0].num() && double(k) == params[1].num(),
654  "Bad parameters");
655  papprox_integration
656  pai = std::make_shared<Newton_Cotes_approx_integration_>(dim_type(n),
657  short_type(k));
658  pintegration_method p = std::make_shared<integration_method>(pai);
659  dependencies.push_back(p->approx_method()->ref_convex());
660  dependencies.push_back(p->approx_method()->pintegration_points());
661  return p;
662  }
663 
664  /* ********************************************************************* */
665  /* integration on direct product of convex structures */
666  /* ********************************************************************* */
667 
668  struct a_int_pro_integration : public approx_integration
669  {
670  a_int_pro_integration(papprox_integration a, papprox_integration b);
671  };
672 
673 
674  a_int_pro_integration::a_int_pro_integration(papprox_integration a,
675  papprox_integration b) {
676  cvr = bgeot::convex_ref_product(a->ref_convex(), b->ref_convex());
677  size_type n1 = a->nb_points_on_convex();
678  size_type n2 = b->nb_points_on_convex();
679  std::vector<base_node> int_points;
680  int_points.resize(n1 * n2);
681  int_coeffs.resize(n1 * n2);
682  repartition.resize(cvr->structure()->nb_faces()+1);
683  repartition[0] = n1 * n2;
684  for (size_type i1 = 0; i1 < n1; ++i1)
685  for (size_type i2 = 0; i2 < n2; ++i2) {
686  int_coeffs[i1 + i2 * n1] = a->coeff(i1) * b->coeff(i2);
687  int_points[i1 + i2 * n1].resize(dim());
688  std::copy(a->point(i1).begin(), a->point(i1).end(),
689  int_points[i1 + i2 * n1].begin());
690  std::copy(b->point(i2).begin(), b->point(i2).end(),
691  int_points[i1 + i2 * n1].begin() + a->dim());
692  }
693  short_type f = 0;
694  for (short_type f1 = 0; f1 < a->structure()->nb_faces(); ++f1, ++f) {
695  n1 = a->nb_points_on_face(f1);
696  n2 = b->nb_points_on_convex();
697  size_type w = repartition[f];
698  repartition[f+1] = w + n1 * n2;
699  int_points.resize(repartition[f+1]);
700  int_coeffs.resize(repartition[f+1]);
701  for (size_type i1 = 0; i1 < n1; ++i1)
702  for (size_type i2 = 0; i2 < n2; ++i2) {
703  int_coeffs[w + i1 + i2 * n1] = a->coeff_on_face(f1, i1)
704  * b->coeff(i2);
705  int_points[w + i1 + i2 * n1].resize(dim());
706  std::copy(a->point_on_face(f1, i1).begin(),
707  a->point_on_face(f1, i1).end(),
708  int_points[w + i1 + i2 * n1].begin());
709  std::copy(b->point(i2).begin(), b->point(i2).end(),
710  int_points[w + i1 + i2 * n1].begin() + a->dim());
711  }
712  }
713  for (short_type f2 = 0; f2 < b->structure()->nb_faces(); ++f2, ++f) {
714  n1 = a->nb_points_on_convex();
715  n2 = b->nb_points_on_face(f2);
716  size_type w = repartition[f];
717  repartition[f+1] = w + n1 * n2;
718  int_points.resize(repartition[f+1]);
719  int_coeffs.resize(repartition[f+1]);
720  for (size_type i1 = 0; i1 < n1; ++i1)
721  for (size_type i2 = 0; i2 < n2; ++i2) {
722  int_coeffs[w + i1 + i2 * n1] = a->coeff(i1)
723  * b->coeff_on_face(f2, i2);
724  int_points[w + i1 + i2 * n1].resize(dim());
725  std::copy(a->point(i1).begin(), a->point(i1).end(),
726  int_points[w + i1 + i2 * n1].begin());
727  std::copy(b->point_on_face(f2, i2).begin(),
728  b->point_on_face(f2, i2).end(),
729  int_points[w + i1 + i2 * n1].begin() + a->dim());
730  }
731  }
732  pint_points = bgeot::store_point_tab(int_points);
733  valid = true;
734  }
735 
736  static pintegration_method
737  product_approx(im_param_list &params,
738  std::vector<dal::pstatic_stored_object> &dependencies) {
739  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
740  << params.size() << " should be 2.");
741  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
742  "Bad type of parameters");
743  pintegration_method a = params[0].method();
744  pintegration_method b = params[1].method();
745  GMM_ASSERT1(a->type() == IM_APPROX && b->type() == IM_APPROX,
746  "Bad parameters");
747  papprox_integration
748  pai = std::make_shared<a_int_pro_integration>(a->approx_method(),
749  b->approx_method());
750  pintegration_method p = std::make_shared<integration_method>(pai);
751  dependencies.push_back(p->approx_method()->ref_convex());
752  dependencies.push_back(p->approx_method()->pintegration_points());
753  return p;
754  }
755 
756  static pintegration_method
757  product_which(im_param_list &params,
758  std::vector<dal::pstatic_stored_object> &dependencies) {
759  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
760  << params.size() << " should be 2.");
761  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
762  "Bad type of parameters");
763  pintegration_method a = params[0].method();
764  pintegration_method b = params[1].method();
765  if (a->type() == IM_EXACT || b->type() == IM_EXACT)
766  return product_exact(params, dependencies);
767  else return product_approx(params, dependencies);
768  }
769 
770 
771  /* ********************************************************************* */
772  /* integration on parallelepiped with Newton Cotes formulae */
773  /* ********************************************************************* */
774 
775  static pintegration_method
776  Newton_Cotes_para(im_param_list &params,
777  std::vector<dal::pstatic_stored_object> &) {
778  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
779  << params.size() << " should be 2.");
780  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
781  "Bad type of parameters");
782  int n = int(::floor(params[0].num() + 0.01));
783  int k = int(::floor(params[1].num() + 0.01));
784  GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
785  double(n) == params[0].num() && double(k) == params[1].num(),
786  "Bad parameters");
787 
788  std::stringstream name;
789  if (n == 1)
790  name << "IM_NC(1," << k << ")";
791  else
792  name << "IM_PRODUCT(IM_NC_PARALLELEPIPED(" << n-1 << "," << k
793  << "),IM_NC(1," << k << "))";
794  return int_method_descriptor(name.str());
795  }
796 
797  static pintegration_method
798  Newton_Cotes_prism(im_param_list &params,
799  std::vector<dal::pstatic_stored_object> &) {
800  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
801  << params.size() << " should be 2.");
802  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
803  "Bad type of parameters");
804  int n = int(::floor(params[0].num() + 0.01));
805  int k = int(::floor(params[1].num() + 0.01));
806  GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
807  double(n) == params[0].num() && double(k) == params[1].num(),
808  "Bad parameters");
809 
810  std::stringstream name;
811  name << "IM_PRODUCT(IM_NC(" << n-1 << "," << k << "),IM_NC(1,"
812  << k << "))";
813  return int_method_descriptor(name.str());
814  }
815 
816  /* ********************************************************************* */
817  /* integration on parallelepiped with Gauss formulae */
818  /* ********************************************************************* */
819 
820  static pintegration_method
821  Gauss_paramul(im_param_list &params,
822  std::vector<dal::pstatic_stored_object> &) {
823  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
824  << params.size() << " should be 2.");
825  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
826  "Bad type of parameters");
827  int n = int(::floor(params[0].num() + 0.01));
828  int k = int(::floor(params[1].num() + 0.01));
829  GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
830  double(n) == params[0].num() && double(k) == params[1].num(),
831  "Bad parameters");
832 
833  std::stringstream name;
834  if (n == 1)
835  name << "IM_GAUSS1D(" << k << ")";
836  else
837  name << "IM_PRODUCT(IM_GAUSS_PARALLELEPIPED(" << n-1 << "," << k
838  << "),IM_GAUSS1D(" << k << "))";
839  return int_method_descriptor(name.str());
840  }
841 
842  /* ******************************************************************** */
843  /* Quasi-polar integration */
844  /* ******************************************************************** */
845 
846  struct quasi_polar_integration : public approx_integration {
847  quasi_polar_integration(papprox_integration base_im,
848  size_type ip1, size_type ip2=size_type(-1)) :
849  approx_integration
850  ((base_im->structure() == bgeot::parallelepiped_structure(3)) ?
851  bgeot::pyramid_QK_of_reference(1)
852  : bgeot::simplex_of_reference(base_im->dim())) {
853  size_type N = base_im->dim();
854 
855  enum { SQUARE, PRISM, TETRA_CYL, PRISM2, PYRAMID } what;
856  if (N == 2) what = SQUARE;
857  else if (base_im->structure() == bgeot::prism_P1_structure(3))
858  what = (ip2 == size_type(-1) || ip1 == ip2) ? PRISM2 : PRISM;
859  else if (base_im->structure() == bgeot::simplex_structure(3))
860  what = TETRA_CYL;
861  else if (base_im->structure() == bgeot::parallelepiped_structure(3))
862  what = PYRAMID;
863  else GMM_ASSERT1(false, "Incoherent integration method");
864 
865  // The first geometric transformation collapse a face of
866  // a parallelepiped or collapse a parrallelepiped on a pyramid.
867  // The second geometric transformation chooses the orientation.
868  // The third is used for the TETRA_CYL case only.
869  bgeot::pgeometric_trans pgt1 = bgeot::parallelepiped_geotrans(N,1);
870  std::vector<base_node> nodes1 = pgt1->convex_ref()->points();
871  bgeot::pgeometric_trans pgt2 = bgeot::simplex_geotrans(N, 1);
872  std::vector<base_node> nodes2(N+1);
873  if (what == PYRAMID) {
874  pgt2 = bgeot::pyramid_QK_geotrans(1);
875  nodes2.resize(5);
876  }
877  std::vector<size_type> other_nodes; // for the construction of node2
878  bgeot::pgeometric_trans pgt3 = bgeot::simplex_geotrans(N, 1);
879  std::vector<base_node> nodes3(N+1);
880 
881  switch (what) {
882  case SQUARE :
883  nodes1[3] = nodes1[1];
884  nodes2[ip1] = nodes1[1]; ip2 = ip1;
885  other_nodes.push_back(0);
886  other_nodes.push_back(2);
887  break;
888  case PRISM :
889  nodes1[4] = nodes1[0]; nodes1[5] = nodes1[1];
890  nodes2[ip1] = nodes1[0];
891  nodes2[ip2] = nodes1[1];
892  other_nodes.push_back(2);
893  other_nodes.push_back(6);
894  break;
895  case TETRA_CYL :
896  nodes3[0] = nodes1[0]; nodes3[1] = nodes1[1];
897  nodes3[2] = nodes1[2]; nodes3[3] = nodes1[4];
898  // nodes1[4] = nodes1[0]; nodes1[7] = base_node(1.0, 1.0, 2.0);
899  nodes2[ip1] = nodes1[1]; ip2 = ip1;
900  other_nodes.push_back(0);
901  other_nodes.push_back(2);
902  other_nodes.push_back(4);
903  break;
904  case PRISM2 :
905  nodes2[ip1] = nodes1[4];
906  other_nodes.push_back(0);
907  other_nodes.push_back(1);
908  other_nodes.push_back(2);
909  break;
910  case PYRAMID :
911  ip2 = ip1 = 0;
912  nodes1[0] = base_node(-1.,-1., 0.);
913  nodes1[1] = base_node( 1.,-1., 0.);
914  nodes1[2] = base_node(-1., 1., 0.);
915  nodes1[3] = base_node( 1., 1., 0.);
916  nodes1[4] = base_node( 0., 0., 1.);
917  nodes1[5] = nodes1[6] = nodes1[7] = nodes1[4];
918  nodes2[ip1] = nodes1[0];
919  other_nodes.push_back(4);
920  other_nodes.push_back(3);
921  other_nodes.push_back(2);
922  other_nodes.push_back(1);
923  }
924 
925  for (size_type i = 0; i <= nodes2.size()-1; ++i)
926  if (i != ip1 && i != ip2) {
927  GMM_ASSERT3(!other_nodes.empty(), "");
928  nodes2[i] = nodes1[other_nodes.back()];
929  other_nodes.pop_back();
930  }
931 
932  base_matrix G1, G2, G3;
933  base_matrix K(N, N), grad(N, N), K3(N, N), K4(N, N);
934  base_node normal1(N), normal2(N);
935  bgeot::geotrans_inv_convex gic(nodes2, pgt2);
936  scalar_type J1, J2, J3(1), J4(1);
937 
938  bgeot::vectors_to_base_matrix(G1, nodes1);
939  bgeot::vectors_to_base_matrix(G2, nodes2);
940 
941  for (size_type nc = 0; nc < 2; ++nc) {
942 
943  if (what == TETRA_CYL) {
944  if (nc == 1) nodes3[0] = nodes1[3];
945  bgeot::vectors_to_base_matrix(G3, nodes3);
946  pgt3->poly_vector_grad(nodes1[0], grad);
947  gmm::mult(gmm::transposed(grad), gmm::transposed(G3), K3);
948  J3 = gmm::abs(gmm::lu_inverse(K3)); /* = 1 */
949  }
950 
951  for (size_type i=0; i < base_im->nb_points(); ++i) {
952 
953  gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
954 
955  size_type fp = size_type(-1); // Search the face number in the
956  if (i >= base_im->nb_points_on_convex()) { // original element
957  size_type ii = i - base_im->nb_points_on_convex();
958  for (short_type ff=0; ff < base_im->structure()->nb_faces(); ++ff) {
959  if (ii < base_im->nb_points_on_face(ff)) { fp = ff; break; }
960  else ii -= base_im->nb_points_on_face(ff);
961  }
962  normal1 = base_im->ref_convex()->normals()[fp];
963  }
964 
965  base_node P = base_im->point(i);
966  if (what == TETRA_CYL) {
967  P = pgt3->transform(P, nodes3);
968  scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
969  K4(0, 1) = - y / one_minus_z;
970  K4(1, 1) = 1.0 - x / one_minus_z;
971  K4(2, 1) = - x * y / gmm::sqr(one_minus_z);
972  J4 = gmm::abs(gmm::lu_det(K4));
973  P[1] *= (1.0 - x / one_minus_z);
974  }
975  if (what == PRISM2) {
976  scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
977  K4(0,0) = one_minus_z; K4(2,0) = -x;
978  K4(1,1) = one_minus_z; K4(2,1) = -y;
979  J4 = gmm::abs(gmm::lu_det(K4));
980  P[0] *= one_minus_z;
981  P[1] *= one_minus_z;
982  }
983 
984  base_node P1 = pgt1->transform(P, nodes1), P2(N);
985  pgt1->poly_vector_grad(P, grad);
986  gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
987  J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
988 
989  if (fp != size_type(-1) && J1 > 1E-10 && J4 > 1E-10) {
990  if (what == TETRA_CYL) {
991  gmm::mult(K3, normal1, normal2);
992  normal1 = normal2;
993  }
994  gmm::lu_inverse(K4);
995  gmm::lu_inverse(K);
996  gmm::mult(K4, normal1, normal2);
997  gmm::mult(K, normal2, normal1);
998  normal2 = normal1;
999  J1 *= gmm::vect_norm2(normal2);
1000  normal2 /= gmm::vect_norm2(normal2);
1001  }
1002  gic.invert(P1, P2);
1003  GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
1004  "Point not in the convex ref : " << P2);
1005 
1006  pgt2->poly_vector_grad(P2, grad);
1007  gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
1008  J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
1009 
1010  if (i < base_im->nb_points_on_convex())
1011  add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
1012  else if (J1 > 1E-10) {
1013  short_type f = short_type(-1);
1014  for (short_type ff = 0; ff <= N; ++ff)
1015  if (gmm::abs(pgt2->convex_ref()->is_in_face(ff, P2)) < 1E-8) {
1016  GMM_ASSERT1(f == short_type(-1),
1017  "An integration point is common to two faces");
1018  f = ff;
1019  }
1020  if (f != short_type(-1)) {
1021  gmm::mult(K, normal2, normal1);
1022  add_point(P2,base_im->coeff(i)*J1*gmm::vect_norm2(normal1)/J2,f);
1023  }
1024  // else { cout << "Point " << P2 << " eliminated" << endl; }
1025  }
1026  }
1027  if (what != TETRA_CYL) break;
1028  }
1029  valid_method();
1030  }
1031  };
1032 
1033 
1034  static pintegration_method
1035  quasi_polar(im_param_list &params,
1036  std::vector<dal::pstatic_stored_object> &dependencies) {
1037  GMM_ASSERT1(params.size() >= 1 && params[0].type() == 1,
1038  "Bad parameters for quasi polar integration: the first "
1039  "parameter should be an integration method");
1040  pintegration_method a = params[0].method();
1041  GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
1042  int ip1 = 0, ip2 = 0;
1043  if (a->approx_method()->structure() == bgeot::parallelepiped_structure(3)) {
1044  GMM_ASSERT1(params.size() == 1, "Bad number of parameters");
1045  } else {
1046  GMM_ASSERT1(params.size() == 2 || params.size() == 3,
1047  "Bad number of parameters : " << params.size()
1048  << " should be 2 or 3.");
1049  GMM_ASSERT1(params[1].type() == 0
1050  && params.back().type() == 0, "Bad type of parameters");
1051  ip1 = int(::floor(params[1].num() + 0.01));
1052  ip2 = int(::floor(params.back().num() + 0.01));
1053  }
1054  int N = a->approx_method()->dim();
1055  GMM_ASSERT1(N >= 2 && N <= 3 && ip1 >= 0 && ip2 >= 0 && ip1 <= N
1056  && ip2 <= N, "Bad parameters");
1057 
1058  papprox_integration
1059  pai = std::make_shared<quasi_polar_integration>(a->approx_method(),
1060  ip1, ip2);
1061  pintegration_method p = std::make_shared<integration_method>(pai);
1062  dependencies.push_back(p->approx_method()->ref_convex());
1063  dependencies.push_back(p->approx_method()->pintegration_points());
1064  return p;
1065  }
1066 
1067  static pintegration_method
1068  pyramid(im_param_list &params,
1069  std::vector<dal::pstatic_stored_object> &dependencies) {
1070  GMM_ASSERT1(params.size() == 1 && params[0].type() == 1,
1071  "Bad parameters for pyramid integration: the first "
1072  "parameter should be an integration method");
1073  pintegration_method a = params[0].method();
1074  GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
1075  int N = a->approx_method()->dim();
1076  GMM_ASSERT1(N == 3, "Bad parameters");
1077 
1078  papprox_integration
1079  pai = std::make_shared<quasi_polar_integration>(a->approx_method(), 0, 0);
1080  pintegration_method p = std::make_shared<integration_method>(pai);
1081  dependencies.push_back(p->approx_method()->ref_convex());
1082  dependencies.push_back(p->approx_method()->pintegration_points());
1083  return p;
1084  }
1085 
1086 
1087 
1088  /* ******************************************************************** */
1089  /* Naming system */
1090  /* ******************************************************************** */
1091 
1092  pintegration_method
1093  structured_composite_int_method(im_param_list &,
1094  std::vector<dal::pstatic_stored_object> &);
1095  pintegration_method HCT_composite_int_method(im_param_list &params,
1096  std::vector<dal::pstatic_stored_object> &dependencies);
1097 
1098  pintegration_method QUADC1_composite_int_method(im_param_list &params,
1099  std::vector<dal::pstatic_stored_object> &dependencies);
1100 
1101  pintegration_method pyramid_composite_int_method(im_param_list &params,
1102  std::vector<dal::pstatic_stored_object> &dependencies);
1103 
1104  struct im_naming_system : public dal::naming_system<integration_method> {
1105  im_naming_system() : dal::naming_system<integration_method>("IM") {
1106  add_suffix("NONE",im_none);
1107  add_suffix("EXACT_SIMPLEX", exact_simplex);
1108  add_suffix("PRODUCT", product_which);
1109  add_suffix("EXACT_PARALLELEPIPED",exact_parallelepiped);
1110  add_suffix("EXACT_PRISM", exact_prism);
1111  add_suffix("GAUSS1D", gauss1d);
1112  add_suffix("NC", Newton_Cotes);
1113  add_suffix("NC_PARALLELEPIPED", Newton_Cotes_para);
1114  add_suffix("NC_PRISM", Newton_Cotes_prism);
1115  add_suffix("GAUSS_PARALLELEPIPED", Gauss_paramul);
1116  add_suffix("QUASI_POLAR", quasi_polar);
1117  add_suffix("PYRAMID", pyramid);
1118  add_suffix("STRUCTURED_COMPOSITE",
1119  structured_composite_int_method);
1120  add_suffix("HCT_COMPOSITE",
1121  HCT_composite_int_method);
1122  add_suffix("QUADC1_COMPOSITE",
1123  QUADC1_composite_int_method);
1124  add_suffix("PYRAMID_COMPOSITE",
1125  pyramid_composite_int_method);
1126  add_generic_function(im_list_integration);
1127  }
1128  };
1129 
1130  pintegration_method int_method_descriptor(std::string name,
1131  bool throw_if_not_found) {
1132  size_type i = 0;
1133  return dal::singleton<im_naming_system>::instance().method
1134  (name, i, throw_if_not_found);
1135  }
1136 
1137  std::string name_of_int_method(pintegration_method p) {
1138  if (!(p.get())) return "IM_NONE";
1139  return dal::singleton<im_naming_system>::instance().shorter_name_of_method(p);
1140  }
1141 
1142  // allows the add of an integration method.
1143  void add_integration_name(std::string name,
1144  dal::naming_system<integration_method>::pfunction f) {
1145  dal::singleton<im_naming_system>::instance().add_suffix(name, f);
1146  }
1147 
1148 
1149  /* Fonctions pour la ref. directe. */
1150 
1151  pintegration_method exact_simplex_im(size_type n) {
1152  THREAD_SAFE_STATIC pintegration_method pim = nullptr;
1153  THREAD_SAFE_STATIC size_type d = -2;
1154  if (d != n) {
1155  std::stringstream name;
1156  name << "IM_EXACT_SIMPLEX(" << n << ")";
1157  pim = int_method_descriptor(name.str());
1158  d = n;
1159  }
1160  return pim;
1161  }
1162 
1163  pintegration_method exact_parallelepiped_im(size_type n) {
1164  THREAD_SAFE_STATIC pintegration_method pim = nullptr;
1165  THREAD_SAFE_STATIC size_type d = -2;
1166  if (d != n) {
1167  std::stringstream name;
1168  name << "IM_EXACT_PARALLELEPIPED(" << n << ")";
1169  pim = int_method_descriptor(name.str());
1170  d = n;
1171  }
1172  return pim;
1173  }
1174 
1175  pintegration_method exact_prism_im(size_type n) {
1176  THREAD_SAFE_STATIC pintegration_method pim = nullptr;
1177  THREAD_SAFE_STATIC size_type d = -2;
1178  if (d != n) {
1179  std::stringstream name;
1180  name << "IM_EXACT_PRISM(" << n << ")";
1181  pim = int_method_descriptor(name.str());
1182  d = n;
1183  }
1184  return pim;
1185  }
1186 
1187  pintegration_method exact_classical_im(bgeot::pgeometric_trans pgt) {
1188  return classical_exact_im(pgt);
1189  }
1190 
1191  static pintegration_method classical_exact_im(bgeot::pconvex_structure cvs) {
1192  cvs = bgeot::basic_structure(cvs);
1193  THREAD_SAFE_STATIC bgeot::pconvex_structure cvs_last = nullptr;
1194  THREAD_SAFE_STATIC pintegration_method im_last = nullptr;
1195  bool found = false;
1196 
1197  if (cvs_last == cvs)
1198  return im_last;
1199 
1200  size_type n = cvs->dim();
1201  size_type nbp = cvs->nb_points();
1202  std::stringstream name;
1203 
1204  /* Identifying P1-simplexes. */
1205 
1206  if (nbp == n+1)
1207  if (cvs == bgeot::simplex_structure(dim_type(n)))
1208  { name << "IM_EXACT_SIMPLEX("; found = true; }
1209 
1210  /* Identifying Q1-parallelepiped. */
1211 
1212  if (!found && nbp == (size_type(1) << n))
1213  if (cvs == bgeot::parallelepiped_structure(dim_type(n)))
1214  { name << "IM_EXACT_PARALLELEPIPED("; found = true; }
1215 
1216  /* Identifying Q1-prisms. */
1217 
1218  if (!found && nbp == 2 * n)
1219  if (cvs == bgeot::prism_P1_structure(dim_type(n)))
1220  { name << "IM_EXACT_PRISM("; found = true; }
1221 
1222  // To be completed
1223 
1224  if (found) {
1225  name << int(n) << ')';
1226  im_last = int_method_descriptor(name.str());
1227  cvs_last = cvs;
1228  return im_last;
1229  }
1230 
1231  GMM_ASSERT1(false, "This element is not taken into account. Contact us");
1232  }
1233 
1234 
1235  pintegration_method classical_exact_im(bgeot::pgeometric_trans pgt) {
1236  return classical_exact_im(pgt->structure());
1237  }
1238 
1239  static pintegration_method
1240  classical_approx_im_(bgeot::pconvex_structure cvs, dim_type degree) {
1241  size_type n = cvs->dim();
1242  std::stringstream name;
1243 
1244  if(bgeot::is_torus_structure(cvs) && n == 3) n = 2;
1245 
1246  degree = std::max<dim_type>(degree, 1);
1247  bgeot::pconvex_structure a, b;
1248  if (bgeot::basic_structure(cvs) == bgeot::simplex_structure(dim_type(n))) {
1249  /* Identifying P1-simplexes. */
1250  switch (n) {
1251  case 0: return int_method_descriptor("IM_NC(0,0)");
1252  case 1: name << "IM_GAUSS1D"; break;
1253  case 2: name << "IM_TRIANGLE"; break;
1254  case 3: name << "IM_TETRAHEDRON"; break;
1255  case 4: name << "IM_SIMPLEX4D"; break;
1256  default: GMM_ASSERT1(false, "no approximate integration method "
1257  "for simplexes of dimension " << n);
1258  }
1259  for (size_type k = degree; k < size_type(degree+10); ++k) {
1260  std::stringstream name2;
1261  name2 << name.str() << "(" << k << ")";
1262  pintegration_method im = int_method_descriptor(name2.str(), false);
1263  if (im) return im;
1264  }
1265  GMM_ASSERT1(false, "could not find an " << name.str()
1266  << " of degree >= " << int(degree));
1267  } else if (bgeot::basic_structure(cvs) == bgeot::pyramid_QK_structure(1)) {
1268  GMM_ASSERT1(n == 3, "Wrong dimension");
1269  name << "IM_PYRAMID(IM_GAUSS_PARALLELEPIPED(3," << degree << "))";
1270  } else if (cvs->is_product(&a,&b) ||
1271  (bgeot::basic_structure(cvs).get() &&
1272  bgeot::basic_structure(cvs)->is_product(&a,&b))) {
1273  name << "IM_PRODUCT("
1274  << name_of_int_method(classical_approx_im_(a,degree)) << ","
1275  << name_of_int_method(classical_approx_im_(b,degree)) << ")";
1276  } else GMM_ASSERT1(false, "unknown convex structure!");
1277  return int_method_descriptor(name.str());
1278  }
1279 
1280  pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt,
1281  dim_type degree) {
1282  THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_last = nullptr;
1283  THREAD_SAFE_STATIC dim_type degree_last;
1284  THREAD_SAFE_STATIC pintegration_method im_last = nullptr;
1285  if (pgt_last == pgt && degree == degree_last)
1286  return im_last;
1287  im_last = classical_approx_im_(pgt->structure(),degree);
1288  degree_last = degree;
1289  pgt_last = pgt;
1290  return im_last;
1291  }
1292 
1293  pintegration_method im_none() {
1294  THREAD_SAFE_STATIC pintegration_method im_last = nullptr;
1295  if (!im_last) im_last = int_method_descriptor("IM_NONE");
1296  return im_last;
1297  }
1298 
1299  /* try to integrate all monomials up to order 'order' and return the
1300  max. error */
1301  scalar_type test_integration_error(papprox_integration pim, dim_type order) {
1302  short_type dim = pim->dim();
1303  pintegration_method exact = classical_exact_im(pim->structure());
1304  opt_long_scalar_type error(0);
1305  for (bgeot::power_index idx(dim); idx.degree() <= order; ++idx) {
1306  opt_long_scalar_type sum(0), realsum;
1307  for (size_type i=0; i < pim->nb_points_on_convex(); ++i) {
1308  opt_long_scalar_type prod = pim->coeff(i);
1309  for (size_type d=0; d < dim; ++d)
1310  prod *= pow(opt_long_scalar_type(pim->point(i)[d]), idx[d]);
1311  sum += prod;
1312  }
1313  realsum = exact->exact_method()->int_monomial(idx);
1314  error = std::max(error, gmm::abs(realsum-sum));
1315  }
1316  return bgeot::to_scalar(error);
1317  }
1318 
1319  papprox_integration get_approx_im_or_fail(pintegration_method pim) {
1320  GMM_ASSERT1(pim->type() == IM_APPROX, "error estimate work only with "
1321  "approximate integration methods");
1322  return pim->approx_method();
1323  }
1324 
1325 } /* 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:978
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:558
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:508
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:397
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:368
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:1163
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:1137
getfem::exact_simplex_im
pintegration_method exact_simplex_im(size_type n)
return IM_EXACT_SIMPLEX(n)
Definition: getfem_integration.cc:1151
getfem::exact_prism_im
pintegration_method exact_prism_im(size_type n)
return IM_EXACT_PRISM(n)
Definition: getfem_integration.cc:1175
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:1280
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:1235
getfem::exact_classical_im
pintegration_method exact_classical_im(bgeot::pgeometric_trans pgt) IS_DEPRECATED
use classical_exact_im instead.
Definition: getfem_integration.cc:1187
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:1130
getfem::im_none
pintegration_method im_none(void)
return IM_NONE
Definition: getfem_integration.cc:1293