GetFEM++  5.3
getfem_integration_composite.cc
1 /*===========================================================================
2 
3  Copyright (C) 2002-2017 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 
23 #include "getfem/bgeot_poly_composite.h"
24 #include "getfem/getfem_integration.h"
25 #include "getfem/getfem_mesh_im.h"
26 #include "getfem/dal_naming_system.h"
27 
28 namespace getfem {
29 
30  papprox_integration
31  composite_approx_int_method(const bgeot::mesh_precomposite &mp,
32  const mesh_im &mi,
33  bgeot::pconvex_ref cr) {
34  auto p = std::make_shared<approx_integration>(cr);
35  base_vector w;
36  for (dal::bv_visitor cv(mi.convex_index()); !cv.finished(); ++cv) {
37  pintegration_method pim = mi.int_method_of_element(cv);
38  bgeot::pgeometric_trans pgt = mi.linked_mesh().trans_of_convex(cv);
39  if (pim->type() != IM_APPROX || !(pgt->is_linear())) {
40  GMM_ASSERT1(false, "Approx integration and linear transformation "
41  "are required.");
42  }
43  papprox_integration pai = pim->approx_method();
44 
45  for (size_type j = 0; j < pai->nb_points_on_convex(); ++j) {
46  base_node pt = pgt->transform((*(pim->pintegration_points()))[j],
47  mi.linked_mesh().points_of_convex(cv));
48  p->add_point(pt, pai->coeff(j) * gmm::abs(mp.det[cv]));
49  }
50  for (short_type f = 0; f < pgt->structure()->nb_faces(); ++f) {
51 
52  base_node barycentre = gmm::mean_value(mi.linked_mesh().points_of_face_of_convex(cv, f).begin(),
53  mi.linked_mesh().points_of_face_of_convex(cv, f).end());
54  short_type f2 = short_type(-1);
55  for (short_type f3 = 0; f3 < cr->structure()->nb_faces(); ++f3) {
56  if (gmm::abs(cr->is_in_face(f3, barycentre)) < 1.0E-7)
57  { f2 = f3; break;}
58  }
59  if (f2 != short_type(-1)) {
60  w.resize(gmm::mat_nrows(mp.gtrans[cv]));
61  gmm::mult(mp.gtrans[cv], pgt->normals()[f], w);
62  scalar_type coeff_mul = gmm::abs(gmm::vect_norm2(w) * mp.det[cv]);
63  for (size_type j = 0; j < pai->nb_points_on_face(f); ++j) {
64  base_node pt = pgt->transform
65  (pai->point_on_face(f, j),
66  mi.linked_mesh().points_of_convex(cv));
67  p->add_point(pt, pai->coeff_on_face(f, j) * coeff_mul, f2);
68  }
69  }
70  }
71  }
72  p->valid_method();
73 
74  return p;
75  }
76 
77  typedef dal::naming_system<integration_method>::param_list im_param_list;
78 
79  pintegration_method structured_composite_int_method(im_param_list &params,
80  std::vector<dal::pstatic_stored_object> &dependencies) {
81  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
82  << params.size() << " should be 2.");
83  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 0,
84  "Bad type of parameters");
85  pintegration_method pim = params[0].method();
86  int k = int(::floor(params[1].num() + 0.01));
87  GMM_ASSERT1(pim->type() == IM_APPROX && k > 0 && k <= 500 &&
88  double(k) == params[1].num(), "Bad parameters");
89 
90  bgeot::pbasic_mesh pm;
91  bgeot::pmesh_precomposite pmp;
92 
93  structured_mesh_for_convex(pim->approx_method()->ref_convex(),
94  short_type(k), pm, pmp);
95  mesh m(*pm);
96  mesh_im mi(m);
97  mi.set_integration_method(pm->convex_index(), pim);
98 
99  pintegration_method
100  p = std::make_shared<integration_method>
101  (composite_approx_int_method(*pmp, mi,
102  pim->approx_method()->ref_convex()));
103  dependencies.push_back(p->approx_method()->ref_convex());
104  dependencies.push_back(p->approx_method()->pintegration_points());
105  return p;
106  }
107 
108 
109 
110  struct just_for_singleton_HCT__ { mesh m; bgeot::mesh_precomposite mp; };
111 
112  pintegration_method HCT_composite_int_method(im_param_list &params,
113  std::vector<dal::pstatic_stored_object> &dependencies) {
114 
115  just_for_singleton_HCT__ &jfs
116  = dal::singleton<just_for_singleton_HCT__>::instance();
117 
118  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
119  << params.size() << " should be 1.");
120  GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
121  pintegration_method pim = params[0].method();
122  GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
123 
124 
125  jfs.m.clear();
126  size_type i0 = jfs.m.add_point(base_node(1.0/3.0, 1.0/3.0));
127  size_type i1 = jfs.m.add_point(base_node(0.0, 0.0));
128  size_type i2 = jfs.m.add_point(base_node(1.0, 0.0));
129  size_type i3 = jfs.m.add_point(base_node(0.0, 1.0));
130  jfs.m.add_triangle(i0, i2, i3);
131  jfs.m.add_triangle(i0, i3, i1);
132  jfs.m.add_triangle(i0, i1, i2);
133  jfs.mp = bgeot::mesh_precomposite(jfs.m);
134 
135  mesh_im mi(jfs.m);
136  mi.set_integration_method(jfs.m.convex_index(), pim);
137 
138  pintegration_method
139  p = std::make_shared<integration_method>
140  (composite_approx_int_method(jfs.mp, mi,
141  pim->approx_method()->ref_convex()));
142  dependencies.push_back(p->approx_method()->ref_convex());
143  dependencies.push_back(p->approx_method()->pintegration_points());
144  return p;
145  }
146 
147 
148  struct just_for_singleton_QUADC1__ { mesh m; bgeot::mesh_precomposite mp; };
149 
150  pintegration_method QUADC1_composite_int_method(im_param_list &params,
151  std::vector<dal::pstatic_stored_object> &dependencies) {
152 
153  just_for_singleton_QUADC1__ &jfs
154  = dal::singleton<just_for_singleton_QUADC1__>::instance();
155 
156  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
157  << params.size() << " should be 1.");
158  GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
159  pintegration_method pim = params[0].method();
160  GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
161 
162  jfs.m.clear();
163  size_type i0 = jfs.m.add_point(base_node(0.0, 0.0));
164  size_type i1 = jfs.m.add_point(base_node(1.0, 0.0));
165  size_type i2 = jfs.m.add_point(base_node(0.0, 1.0));
166  size_type i3 = jfs.m.add_point(base_node(1.0, 1.0));
167  size_type i4 = jfs.m.add_point(base_node(0.5, 0.5));
168  jfs.m.add_triangle(i1, i3, i4);
169  jfs.m.add_triangle(i2, i0, i4);
170  jfs.m.add_triangle(i3, i2, i4);
171  jfs.m.add_triangle(i0, i1, i4);
172  jfs.mp = bgeot::mesh_precomposite(jfs.m);
173 
174  mesh_im mi(jfs.m);
175  mi.set_integration_method(jfs.m.convex_index(), pim);
176 
177  pintegration_method
178  p = std::make_shared<integration_method>
179  (composite_approx_int_method(jfs.mp, mi,
180  bgeot::parallelepiped_of_reference(2)));
181  dependencies.push_back(p->approx_method()->ref_convex());
182  dependencies.push_back(p->approx_method()->pintegration_points());
183  return p;
184  }
185 
186  struct just_for_singleton_pyramidc__ { mesh m; bgeot::mesh_precomposite mp; };
187 
188  pintegration_method pyramid_composite_int_method(im_param_list &params,
189  std::vector<dal::pstatic_stored_object> &dependencies) {
190 
191  just_for_singleton_pyramidc__ &jfs
192  = dal::singleton<just_for_singleton_pyramidc__>::instance();
193 
194  GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
195  << params.size() << " should be 1.");
196  GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
197  pintegration_method pim = params[0].method();
198  GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
199 
200  jfs.m.clear();
201  size_type i0 = jfs.m.add_point(base_node(-1.0, -1.0, 0.0));
202  size_type i1 = jfs.m.add_point(base_node( 1.0, -1.0, 0.0));
203  size_type i2 = jfs.m.add_point(base_node(-1.0, 1.0, 0.0));
204  size_type i3 = jfs.m.add_point(base_node( 1.0, 1.0, 0.0));
205  size_type i4 = jfs.m.add_point(base_node( 0.0, 0.0, 1.0));
206  jfs.m.add_tetrahedron(i0, i1, i2, i4);
207  jfs.m.add_tetrahedron(i1, i3, i2, i4);
208  jfs.mp = bgeot::mesh_precomposite(jfs.m);
209 
210  mesh_im mi(jfs.m);
211  mi.set_integration_method(jfs.m.convex_index(), pim);
212 
213  pintegration_method
214  p = std::make_shared<integration_method>
215  (composite_approx_int_method(jfs.mp, mi,
216  bgeot::pyramid_QK_of_reference(1)));
217  dependencies.push_back(p->approx_method()->ref_convex());
218  dependencies.push_back(p->approx_method()->pintegration_points());
219  return p;
220  }
221 
222 
223 
224 } /* end of namespace getfem. */
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:470
bgeot_poly_composite.h
Handle composite polynomials.
getfem_mesh_im.h
Define the getfem::mesh_im class (integration of getfem::mesh_fem).
getfem_integration.h
Integration methods (exact and approximated) on convexes.
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:146
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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:454
dal_naming_system.h
Naming system.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_assembling.h:48
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:79
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:680
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:183