GetFEM  5.4.2
getfem_integration.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_integration.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date December 17, 2000.
35  @brief Integration methods (exact and approximated) on convexes
36 
37  @section im_list List of integration methods for getfem::int_method_descriptor
38 
39  - "IM_EXACT_SIMPLEX(n)" : exact integration on simplexes.
40  - "IM_PRODUCT(IM1, IM2)" : product of two integration methods
41  - "IM_EXACT_PARALLELEPIPED(n)" : exact integration on parallelepipeds
42  - "IM_EXACT_PRISM(n)" : exact integration on prisms
43  - "IM_GAUSS1D(K)" : Gauss method on the segment,
44  order K=1, 3, 5, 7, ..., 99.
45  - "IM_GAUSSLOBATTO1D(K)" : Gauss-Lobatto method on the segment,
46  order K=1, 3, 5, 7, ..., 99.
47  - "IM_NC(N,K)" : Newton-Cotes approximative
48  integration on simplexes, order K
49  - "IM_NC_PARALLELEPIPED(N,K)" : product of Newton-Cotes integration
50  on parallelepipeds
51  - "IM_NC_PRISM(N,K)" : product of Newton-Cotes integration
52  on prisms
53  - "IM_GAUSS_PARALLELEPIPED(N,K)" : product of Gauss1D integration
54  on parallelepipeds (K=1,3,..99)
55  - "IM_TRIANGLE(K)" : Gauss methods on triangles
56  (K=1..10,13,17,19)
57  - "IM_QUAD(K)" : Gauss methods on quadrilaterons
58  (K=2, 3, 5, 7, 9, 17)
59  - "IM_HEXAHEDRON(K)" : Gauss methods on hexahedrons
60  (K=5,9,11)
61  - "IM_TETRAHEDRON(K)" : Gauss methods on tetrahedrons
62  (K=1, 2, 3, 5, 6, or 8)
63  - "IM_SIMPLEX4D(3)" : Gauss method on a 4-dimensional
64  simplex.
65  - "IM_CUBE4D(K)" : Gauss method on a 4-dimensional cube
66  (K=5,9).
67  - "IM_STRUCTURED_COMPOSITE(IM1, K)" : Composite method on a grid with
68  K divisions
69  - "IM_HCT_COMPOSITE(IM1)" : Composite integration suited to the
70  HCT composite finite element.
71  - "IM_QUAQC1_COMPOSITE(IM1)" : Composite integration suited to the
72  QUADC1 composite finite element.
73 
74  - "IM_QUASI_POLAR(IM1, IP1)" : if IM1 is an integration method on a
75  square, gives an integration method on a triangle which is
76  close to a polar integration with respect to vertex IP1.
77  if IM1 is an integration method on a tetrahedron, gives an
78  integration method on a tetrahedron which is close to a
79  cylindrical integration with respect to vertex IP1
80  (does not work very well).
81  if IM1 is an integration method on a prism. Gives an integration
82  method on a tetrahedron which is close to a
83  cylindrical integration with respect to vertex IP1.
84 
85  - "IM_QUASI_POLAR(IM1, IP1, IP2)" : IM1 should be an integration method
86  on a prism. Gives an integration method on a tetrahedron which
87  is close to a cylindrical integration with respect to IP1-IP2
88  axis.
89 
90  - "IM_PYRAMID_COMPOSITE(IM1)" : Composite integration for a pyramid
91  decomposed into two tetrahedrons
92 */
93 #ifndef GETFEM_INTEGRATION_H__
94 #define GETFEM_INTEGRATION_H__
95 
96 #include "getfem_config.h"
97 #include "bgeot_convex_ref.h"
98 #include "bgeot_geometric_trans.h"
99 #include "bgeot_node_tab.h"
100 #include "bgeot_poly_composite.h"
102 
103 
104 
105 
106 namespace getfem
107 {
108  /** Description of an exact integration of polynomials.
109  * This class is not to be manipulate by itself. Use ppoly_integration and
110  * the functions written to produce the basic descriptions.
111  */
113  protected :
114 
116 
117  mutable std::vector<long_scalar_type> int_monomials;
118  mutable std::vector< std::vector<long_scalar_type> > int_face_monomials;
119 
120  public :
121 
122  /// Dimension of convex of reference.
123  dim_type dim(void) const { return cvs->dim(); }
124  /// {Structure of convex of reference.
125  bgeot::pconvex_structure structure(void) const { return cvs; }
126 
127  virtual long_scalar_type int_monomial(const bgeot::power_index &power)
128  const = 0;
129  virtual long_scalar_type int_monomial_on_face(const bgeot::power_index
130  &power, short_type f) const = 0;
131 
132  /// Evaluate the integral of the polynomial P on the reference element.
133  long_scalar_type int_poly(const base_poly &P) const;
134  /** Evaluate the integral of the polynomial P on the face f of the
135  * reference element.
136  */
137  long_scalar_type int_poly_on_face(const base_poly &P,
138  short_type f) const;
139 
140  virtual ~poly_integration() {}
141 
142  };
143 
144  typedef std::shared_ptr<const poly_integration> ppoly_integration;
145 
146  /** Description of an approximate integration of polynomials of
147  * several variables on a reference element.
148  * This class is not to be manipulate by itself. Use papprox\_integration
149  * and the functions written to produce the basic descriptions.
150  */
151 
152  class integration_method;
153  typedef std::shared_ptr<const integration_method> pintegration_method;
154 
155 
156  class approx_integration {
157  protected :
158 
159  typedef bgeot::node_tab PT_TAB;
160  bgeot::pconvex_ref cvr;
161  bgeot::pstored_point_tab pint_points;
162  std::vector<scalar_type> int_coeffs;
163  std::vector<size_type> repartition;
164 
165  // index 0 : points for volumic integration, index > 0 : points for faces
166  std::vector<PT_TAB> pt_to_store;
167  bool valid;
168  bool built_on_the_fly;
169 
170  public :
171 
172  /// Dimension of reference convex.
173  dim_type dim(void) const { return cvr->structure()->dim(); }
174  size_type nb_points(void) const { return int_coeffs.size(); }
175  /// Number of integration nodes on the reference element.
176  size_type nb_points_on_convex(void) const { return repartition[0]; }
177  /// Number of integration nodes on the face f of the reference element.
178  size_type nb_points_on_face(short_type f) const
179  { return repartition[f+1] - repartition[f]; }
180  size_type ind_first_point_on_face(short_type f) const
181  { return repartition[f]; }
182  /// Convenience method returning the number of faces of the reference convex
183  short_type nb_convex_faces() const
184  { return cvr->structure()->nb_faces(); /* == repartition.size() - 1;*/ }
185  bool is_built_on_the_fly(void) const { return built_on_the_fly; }
186  void set_built_on_the_fly(void) { built_on_the_fly = true; }
187  /// Structure of the reference element.
188  bgeot::pconvex_structure structure(void) const
189  { return basic_structure(cvr->structure()); }
190  bgeot::pconvex_ref ref_convex(void) const { return cvr; }
191 
192  const std::vector<size_type> &repart(void) const { return repartition; }
193 
194  /// Gives an array of integration nodes.
195  // const bgeot::stored_point_tab &integration_points(void) const
196  // { return *(pint_points); }
197  bgeot::pstored_point_tab pintegration_points(void) const
198  { return pint_points; }
199  /// Gives the integration node i on the reference element.
200  const base_node &point(size_type i) const
201  { return (*pint_points)[i]; }
202  /// Gives the integration node i of the face f.
203  const base_node &
204  point_on_face(short_type f, size_type i) const
205  { return (*pint_points)[repartition[f] + i]; }
206  /// Gives an array of the integration coefficients.
207  const std::vector<scalar_type> &integration_coefficients(void) const
208  { return int_coeffs; }
209  /// Gives the integration coefficient corresponding to node i.
210  scalar_type coeff(size_type i) const { return int_coeffs[i]; }
211  /// Gives the integration coefficient corresponding to node i of face f.
212  scalar_type coeff_on_face(short_type f, size_type i) const
213  { return int_coeffs[repartition[f] + i]; }
214 
215  void add_point(const base_node &pt, scalar_type w,
216  short_type f = short_type(-1),
217  bool include_empty = false);
218  void add_point_norepeat(const base_node &pt, scalar_type w,
219  short_type f=short_type(-1));
220  void add_point_full_symmetric(base_node pt, scalar_type w);
221  void add_method_on_face(pintegration_method ppi, short_type f);
222  void valid_method(void);
223 
224  approx_integration(void) : valid(false), built_on_the_fly(false) { }
225  approx_integration(bgeot::pconvex_ref cr)
226  : cvr(cr), repartition(cr->structure()->nb_faces()+1),
227  pt_to_store(cr->structure()->nb_faces()+1), valid(false),
228  built_on_the_fly(false)
229  { std::fill(repartition.begin(), repartition.end(), 0); }
230  virtual ~approx_integration() {}
231  };
232 
233  typedef std::shared_ptr<const approx_integration> papprox_integration;
234 
235  /**
236  the list of main integration method types
237  */
238  typedef enum { IM_APPROX, IM_EXACT, IM_NONE } integration_method_type;
239 
240  /**
241  this structure is not intended to be used directly. It is built via
242  the int_method_descriptor() function
243  */
245  ppoly_integration ppi; /* for exact integrations */
246  papprox_integration pai; /* for approximate integrations (i.e. cubatures) */
247  integration_method_type im_type;
248  void remove() { pai = papprox_integration(); ppi = ppoly_integration(); }
249 
250  public:
251  integration_method_type type(void) const { return im_type; }
252  const papprox_integration &approx_method(void) const { return pai; }
253  const ppoly_integration &exact_method(void) const { return ppi; }
254 
255  void set_approx_method(papprox_integration paii)
256  { remove(); pai = paii; im_type = IM_APPROX; }
257  void set_exact_method(ppoly_integration ppii)
258  { remove(); ppi = ppii; im_type = IM_EXACT; }
259 
260  bgeot::pstored_point_tab pintegration_points(void) const {
261  if (type() == IM_EXACT) {
262  size_type n = ppi->structure()->dim();
263  std::vector<base_node> spt(1); spt[0] = base_node(n);
264  return store_point_tab(spt);
265  }
266  else if (type() == IM_APPROX)
267  return pai->pintegration_points();
268  else GMM_ASSERT1(false, "IM_NONE has no points");
269  }
270 
271  // const bgeot::stored_point_tab &integration_points(void) const
272  // { return *(pintegration_points()); }
273 
274  bgeot::pconvex_structure structure(void) const {
275  switch (type()) {
276  case IM_EXACT: return ppi->structure();
277  case IM_APPROX: return pai->structure();
278  case IM_NONE: GMM_ASSERT1(false, "IM_NONE has no structure");
279  default: GMM_ASSERT3(false, "");
280  }
281  return 0;
282  }
283 
284  integration_method(ppoly_integration p) {
285  DAL_STORED_OBJECT_DEBUG_CREATED(this, "Exact integration method");
286  ppi = p; im_type = IM_EXACT;
287  }
288 
289  integration_method(papprox_integration p) {
290  DAL_STORED_OBJECT_DEBUG_CREATED(this, "Approximate integration method");
291  pai = p; im_type = IM_APPROX;
292  }
293 
295  DAL_STORED_OBJECT_DEBUG_CREATED(this, "Integration method");
296  im_type = IM_NONE;
297  }
298 
299  virtual ~integration_method()
300  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Integration method"); }
301  };
302 
303 
304  /** Get an integration method from its name .
305  @see @ref im_list
306  @param name the integration method name, for example "IM_TRIANGLE(6)"
307  @param throw_if_not_found choose if an exception must be thrown
308  when the integration method does not exist (if no exception, a
309  null pointer is returned).
310  */
311  pintegration_method int_method_descriptor(std::string name,
312  bool throw_if_not_found = true);
313 
314  /** Get the string name of an integration method .
315  */
316  std::string name_of_int_method(pintegration_method p);
317 
318  /**
319  return an exact integration method for convex type handled by pgt.
320  If pgt is not linear, classical_exact_im will fail.
321  */
322  pintegration_method classical_exact_im(bgeot::pgeometric_trans pgt);
323  /**
324  try to find an approximate integration method for the geometric
325  transformation pgt which is able to integrate exactly polynomials
326  of degree <= "degree". It may return a higher order integration
327  method if no method match the exact degree.
328  */
329  pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree);
330 
331  /// return IM_EXACT_SIMPLEX(n)
332  pintegration_method exact_simplex_im(size_type n);
333  /// return IM_EXACT_PARALLELEPIPED(n)
334  pintegration_method exact_parallelepiped_im(size_type n);
335  /// return IM_EXACT_PRISM(n)
336  pintegration_method exact_prism_im(size_type n);
337  /// use classical_exact_im instead.
338  pintegration_method exact_classical_im(bgeot::pgeometric_trans pgt) IS_DEPRECATED;
339  /// return IM_NONE
340  pintegration_method im_none(void);
341 
342 
343  class mesh_im;
344  papprox_integration composite_approx_int_method(const bgeot::mesh_precomposite &mp,
345  const mesh_im &mf,
346  bgeot::pconvex_ref cr);
347 
348  /* try to integrate all monomials up to order 'order' and return the
349  max. error */
350  scalar_type test_integration_error(papprox_integration pim, dim_type order);
351 
352  papprox_integration get_approx_im_or_fail(pintegration_method pim);
353 
354  /* Function allowing the add of an integration method outwards
355  of getfem_integration.cc */
356 
357  typedef dal::naming_system<integration_method>::param_list im_param_list;
358 
359  void add_integration_name(std::string name,
361 
362 } /* end of namespace getfem. */
363 
364 
365 #endif /* BGEOT_INTEGRATION_H__ */
bgeot_convex_ref.h
Reference convexes.
getfem_config.h
defines and typedefs for namespace getfem
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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
bgeot_node_tab.h
Structure which dynamically collects points identifying points that are nearer than a certain very sm...
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
dal::naming_system
Associate a name to a method descriptor and store method descriptors.
Definition: dal_naming_system.h:56
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::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
bgeot::basic_structure
pconvex_structure basic_structure(pconvex_structure cv)
Original structure (if concerned)
Definition: bgeot_convex_structure.h:172
dal_naming_system.h
Naming system.
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
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
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::exact_parallelepiped_im
pintegration_method exact_parallelepiped_im(size_type n)
return IM_EXACT_PARALLELEPIPED(n)
Definition: getfem_integration.cc:1163
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::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::integration_method
this structure is not intended to be used directly.
Definition: getfem_integration.h:244
bgeot::power_index
Vector of integer (16 bits type) which represent the powers of a monomial.
Definition: bgeot_poly.h:64
bgeot_geometric_trans.h
Geometric transformations on convexes.
bgeot::node_tab
Store a set of points, identifying points that are nearer than a certain very small distance.
Definition: bgeot_node_tab.h:52
getfem::im_none
pintegration_method im_none(void)
return IM_NONE
Definition: getfem_integration.cc:1293
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
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
dal::static_stored_object
base class for static stored objects
Definition: dal_static_stored_objects.h:206
getfem::poly_integration::dim
dim_type dim(void) const
Dimension of convex of reference.
Definition: getfem_integration.h:123
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
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
bgeot_poly_composite.h
Handle composite polynomials.
getfem::poly_integration
Description of an exact integration of polynomials.
Definition: getfem_integration.h:112
getfem::integration_method_type
integration_method_type
the list of main integration method types
Definition: getfem_integration.h:238