GetFEM  5.5
getfem_integration_composite.cc
1 /*===========================================================================
2 
3  Copyright (C) 2002-2026 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program. If not, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
21 
22 #include "getfem/bgeot_poly_composite.h"
23 #include "getfem/getfem_integration.h"
24 #include "getfem/getfem_mesh_im.h"
25 #include "getfem/dal_naming_system.h"
26 
27 namespace getfem {
28 
29  papprox_integration
30  composite_approx_int_method(const bgeot::mesh_precomposite &mp,
31  const mesh_im &mi,
32  bgeot::pconvex_ref cr) {
33  auto p = std::make_shared<approx_integration>(cr);
34  base_vector w;
35  for (dal::bv_visitor cv(mi.convex_index()); !cv.finished(); ++cv) {
36  pintegration_method pim = mi.int_method_of_element(cv);
37  bgeot::pgeometric_trans pgt = mi.linked_mesh().trans_of_convex(cv);
38  if (pim->type() != IM_APPROX || !(pgt->is_linear())) {
39  GMM_ASSERT1(false, "Approx integration and linear transformation "
40  "are required.");
41  }
42  papprox_integration pai = pim->approx_method();
43 
44  for (size_type j = 0; j < pai->nb_points_on_convex(); ++j) {
45  base_node pt = pgt->transform((*(pim->pintegration_points()))[j],
46  mi.linked_mesh().points_of_convex(cv));
47  p->add_point(pt, pai->coeff(j) * gmm::abs(mp.det[cv]));
48  }
49  for (short_type f = 0; f < pgt->structure()->nb_faces(); ++f) {
50 
51  base_node barycentre = gmm::mean_value(mi.linked_mesh().points_of_face_of_convex(cv, f).begin(),
52  mi.linked_mesh().points_of_face_of_convex(cv, f).end());
53  short_type f2 = short_type(-1);
54  for (short_type f3 = 0; f3 < cr->structure()->nb_faces(); ++f3) {
55  if (gmm::abs(cr->is_in_face(f3, barycentre)) < 1.0E-7)
56  { f2 = f3; break;}
57  }
58  if (f2 != short_type(-1)) {
59  w.resize(gmm::mat_nrows(mp.gtransinv[cv]));
60  gmm::mult(mp.gtransinv[cv], pgt->normals()[f], w);
61  scalar_type coeff_mul = gmm::abs(gmm::vect_norm2(w) * mp.det[cv]);
62  for (size_type j = 0; j < pai->nb_points_on_face(f); ++j) {
63  base_node pt = pgt->transform
64  (pai->point_on_face(f, j),
65  mi.linked_mesh().points_of_convex(cv));
66  p->add_point(pt, pai->coeff_on_face(f, j) * coeff_mul, f2);
67  }
68  }
69  }
70  }
71  p->valid_method();
72 
73  return p;
74  }
75 
76  typedef dal::naming_system<integration_method>::param_list im_param_list;
77 
78  pintegration_method structured_composite_int_method(im_param_list &params,
79  std::vector<dal::pstatic_stored_object> &dependencies) {
80  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
81  << params.size() << " should be 2.");
82  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 0,
83  "Bad type of parameters");
84  pintegration_method pim = params[0].method();
85  int k = int(::floor(params[1].num() + 0.01));
86  GMM_ASSERT1(pim->type() == IM_APPROX && k > 0 && k <= 500 &&
87  double(k) == params[1].num(), "Bad parameters");
88 
89  bgeot::pbasic_mesh pm;
90  bgeot::pmesh_precomposite pmp;
91 
92  structured_mesh_for_convex(pim->approx_method()->ref_convex(),
93  short_type(k), pm, pmp);
94  mesh m(*pm);
95  mesh_im mi(m);
96  mi.set_integration_method(pm->convex_index(), pim);
97 
98  pintegration_method
99  p = std::make_shared<integration_method>
100  (composite_approx_int_method(*pmp, mi,
101  pim->approx_method()->ref_convex()));
102  dependencies.push_back(p->approx_method()->ref_convex());
103  dependencies.push_back(p->approx_method()->pintegration_points());
104  return p;
105  }
106 
107 
108 
109  struct just_for_singleton_HCT__ { mesh m; bgeot::mesh_precomposite mp; };
110 
111  pintegration_method HCT_composite_int_method(im_param_list &params,
112  std::vector<dal::pstatic_stored_object> &dependencies) {
113 
114  just_for_singleton_HCT__ &jfs
115  = dal::singleton<just_for_singleton_HCT__>::instance();
116 
117  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
118  << params.size() << " should be 1.");
119  GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
120  pintegration_method pim = params[0].method();
121  GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
122 
123 
124  jfs.m.clear();
125  size_type i0 = jfs.m.add_point(base_node(1.0/3.0, 1.0/3.0));
126  size_type i1 = jfs.m.add_point(base_node(0.0, 0.0));
127  size_type i2 = jfs.m.add_point(base_node(1.0, 0.0));
128  size_type i3 = jfs.m.add_point(base_node(0.0, 1.0));
129  jfs.m.add_triangle(i0, i2, i3);
130  jfs.m.add_triangle(i0, i3, i1);
131  jfs.m.add_triangle(i0, i1, i2);
132  jfs.mp.initialise(jfs.m);
133 
134  mesh_im mi(jfs.m);
135  mi.set_integration_method(jfs.m.convex_index(), pim);
136 
137  pintegration_method
138  p = std::make_shared<integration_method>
139  (composite_approx_int_method(jfs.mp, mi,
140  pim->approx_method()->ref_convex()));
141  dependencies.push_back(p->approx_method()->ref_convex());
142  dependencies.push_back(p->approx_method()->pintegration_points());
143  return p;
144  }
145 
146 
147  struct just_for_singleton_QUADC1__ { mesh m; bgeot::mesh_precomposite mp; };
148 
149  pintegration_method QUADC1_composite_int_method(im_param_list &params,
150  std::vector<dal::pstatic_stored_object> &dependencies) {
151 
152  just_for_singleton_QUADC1__ &jfs
153  = dal::singleton<just_for_singleton_QUADC1__>::instance();
154 
155  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
156  << params.size() << " should be 1.");
157  GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
158  pintegration_method pim = params[0].method();
159  GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
160 
161  jfs.m.clear();
162  size_type i0 = jfs.m.add_point(base_node(0.0, 0.0));
163  size_type i1 = jfs.m.add_point(base_node(1.0, 0.0));
164  size_type i2 = jfs.m.add_point(base_node(0.0, 1.0));
165  size_type i3 = jfs.m.add_point(base_node(1.0, 1.0));
166  size_type i4 = jfs.m.add_point(base_node(0.5, 0.5));
167  jfs.m.add_triangle(i1, i3, i4);
168  jfs.m.add_triangle(i2, i0, i4);
169  jfs.m.add_triangle(i3, i2, i4);
170  jfs.m.add_triangle(i0, i1, i4);
171  jfs.mp.initialise(jfs.m);
172 
173  mesh_im mi(jfs.m);
174  mi.set_integration_method(jfs.m.convex_index(), pim);
175 
176  pintegration_method
177  p = std::make_shared<integration_method>
178  (composite_approx_int_method(jfs.mp, mi,
179  bgeot::parallelepiped_of_reference(2)));
180  dependencies.push_back(p->approx_method()->ref_convex());
181  dependencies.push_back(p->approx_method()->pintegration_points());
182  return p;
183  }
184 
185  struct just_for_singleton_pyramidc__ { mesh m; bgeot::mesh_precomposite mp; };
186 
187  pintegration_method pyramid_composite_int_method(im_param_list &params,
188  std::vector<dal::pstatic_stored_object> &dependencies) {
189 
190  just_for_singleton_pyramidc__ &jfs
191  = dal::singleton<just_for_singleton_pyramidc__>::instance();
192 
193  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
194  << params.size() << " should be 1.");
195  GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
196  pintegration_method pim = params[0].method();
197  GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
198 
199  jfs.m.clear();
200  size_type i0 = jfs.m.add_point(base_node(-1.0, -1.0, 0.0));
201  size_type i1 = jfs.m.add_point(base_node( 1.0, -1.0, 0.0));
202  size_type i2 = jfs.m.add_point(base_node(-1.0, 1.0, 0.0));
203  size_type i3 = jfs.m.add_point(base_node( 1.0, 1.0, 0.0));
204  size_type i4 = jfs.m.add_point(base_node( 0.0, 0.0, 1.0));
205  jfs.m.add_tetrahedron(i0, i1, i2, i4);
206  jfs.m.add_tetrahedron(i1, i3, i2, i4);
207  jfs.mp.initialise(jfs.m);
208 
209  mesh_im mi(jfs.m);
210  mi.set_integration_method(jfs.m.convex_index(), pim);
211 
212  pintegration_method
213  p = std::make_shared<integration_method>
214  (composite_approx_int_method(jfs.mp, mi,
215  bgeot::pyramid_QK_of_reference(1)));
216  dependencies.push_back(p->approx_method()->ref_convex());
217  dependencies.push_back(p->approx_method()->pintegration_points());
218  return p;
219  }
220 
221 
222 
223 } /* end of namespace getfem. */
bgeot_poly_composite.h
Handle composite polynomials.
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:164
dal_naming_system.h
Naming system.
getfem_integration.h
Integration methods (exact and approximated) on convexes.
getfem_mesh_im.h
Define the getfem::mesh_im class (integration of getfem::mesh_fem).
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::pyramid_QK_of_reference
pconvex_ref pyramid_QK_of_reference(dim_type k)
pyramidal element of reference of degree k (k = 1 or 2 only)
Definition: bgeot_convex_ref.cc:512
bgeot::structured_mesh_for_convex
void structured_mesh_for_convex(pconvex_ref cvr, short_type k, pbasic_mesh &pm, pmesh_precomposite &pmp, bool force_simplexification)
This function returns a mesh in pm which contains a refinement of the convex cvr if force_simplexific...
Definition: bgeot_poly_composite.cc:498
bgeot::parallelepiped_of_reference
pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k)
parallelepiped of reference of dimension nc (and degree 1)
Definition: bgeot_convex_ref.cc:737
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:185
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46