GetFEM  5.4.2
getfem_regular_meshes.cc
1 /*===========================================================================
2 
3  Copyright (C) 1999-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 
24 
25 namespace getfem
26 {
27  void parallelepiped_regular_simplex_mesh_
28  (mesh &me, dim_type N, const base_node &org,
29  const base_small_vector *ivect, const size_type *iref) {
31  pararef = *(bgeot::parallelepiped_of_reference(N));
32 
33  if (N >= 5) GMM_WARNING1("CAUTION : Simplexification in dimension >= 5 "
34  "has not been tested and the resulting mesh "
35  "should be not conformal");
36 
37  const bgeot::mesh_structure &sl
38  = *(bgeot::parallelepiped_of_reference(N)->simplexified_convex());
39 
40  base_node a = org;
41  size_type i, nbpt = pararef.nb_points();
42 
43  for (i = 0; i < nbpt; ++i) {
44  gmm::clear(a);
45  for (dim_type n = 0; n < N; ++n)
46  gmm::add(gmm::scaled(ivect[n],pararef.points()[i][n]),a);
47  pararef.points()[i] = a;
48  }
49 
50  // bgeot::simplexify(cvt, sl, pararef.points(), N, me.eps());
51 
52  size_type nbs = sl.nb_convex();
53  std::vector<size_type> tab1(N+1), tab(N), tab3(nbpt);
54  size_type total = 0;
55  std::fill(tab.begin(), tab.end(), 0);
56  while (tab[N-1] != iref[N-1]) {
57  for (a = org, i = 0; i < N; i++)
58  gmm::add(gmm::scaled(ivect[i],scalar_type(tab[i])),a);
59  //a.addmul(scalar_type(tab[i]), ivect[i]);
60 
61  for (i = 0; i < nbpt; i++)
62  tab3[i] = me.add_point(a + pararef.points()[i]);
63 
64  for (i = 0; i < nbs; i++) {
65  const mesh::ind_cv_ct &tab2 = sl.ind_points_of_convex(i);
66  for (dim_type l = 0; l <= N; l++)
67  // tab1[l] = tab3[tab2[l]];
68  tab1[l] = tab3[(tab2[l]
69  + (((total & 1) && N != 3) ? (nbpt/2) : 0)) % nbpt];
70  me.add_simplex(N, tab1.begin());
71  }
72 
73  for (dim_type l = 0; l < N; l++) {
74  tab[l]++; total++;
75  if (l < N-1 && tab[l] >= iref[l]) { total -= tab[l]; tab[l] = 0; }
76  else break;
77  }
78  }
79  }
80 
81 
82  void parallelepiped_regular_prism_mesh_
83  (mesh &me, dim_type N, const base_node &org,
84  const base_small_vector *ivect, const size_type *iref) {
85  mesh aux;
86  parallelepiped_regular_simplex_mesh_(aux, dim_type(N-1), org, ivect, iref);
87  std::vector<base_node> ptab(2 * N);
88 
89  for (dal::bv_visitor cv(aux.convex_index()); !cv.finished(); ++cv) {
90  std::copy(aux.points_of_convex(cv).begin(),
91  aux.points_of_convex(cv).end(), ptab.begin());
92 
93  for (size_type k = 0; k < iref[N-1]; ++k) {
94 
95  for (dim_type j = 0; j < N; ++j) ptab[j+N] = ptab[j] + ivect[N-1];
96  me.add_prism_by_points(N, ptab.begin());
97 
98  std::copy(ptab.begin()+N, ptab.end(), ptab.begin());
99  }
100  }
101  }
102 
103  void parallelepiped_regular_mesh_
104  (mesh &me, dim_type N, const base_node &org,
105  const base_small_vector *ivect, const size_type *iref, bool linear_gt) {
107  pararef = *(bgeot::parallelepiped_of_reference(N));
108  base_node a = org;
109  size_type i, nbpt = pararef.nb_points();
110 
111  for (i = 0; i < nbpt; ++i) {
112  gmm::clear(a);
113  for (dim_type n = 0; n < N; ++n)
114  gmm::add(gmm::scaled(ivect[n],pararef.points()[i][n]),a);
115  //a.addmul(pararef.points()[i][n], ivect[n]);
116  pararef.points()[i] = a;
117  }
118 
119  std::vector<size_type> tab1(N+1), tab(N), tab3(nbpt);
120  size_type total = 0;
121  std::fill(tab.begin(), tab.end(), 0);
122  while (tab[N-1] != iref[N-1]) {
123  for (a = org, i = 0; i < N; i++)
124  gmm::add(gmm::scaled(ivect[i], scalar_type(tab[i])),a);
125  //a.addmul(scalar_type(tab[i]), ivect[i]);
126 
127  for (i = 0; i < nbpt; i++)
128  tab3[i] = me.add_point(a + pararef.points()[i]);
129  me.add_convex(linear_gt ?
130  bgeot::parallelepiped_linear_geotrans(N) :
131  bgeot::parallelepiped_geotrans(N, 1), tab3.begin());
132 
133  for (dim_type l = 0; l < N; l++) {
134  tab[l]++; total++;
135  if (l < N-1 && tab[l] >= iref[l]) { total -= tab[l]; tab[l] = 0; }
136  else break;
137  }
138  }
139  }
140 
141  void parallelepiped_regular_pyramid_mesh_
142  (mesh &me, const base_node &org,
143  const base_small_vector *ivect, const size_type *iref) {
144  dim_type N = 3;
146  pararef = *(bgeot::parallelepiped_of_reference(N));
147  base_node a = org, barycenter(N);
148  size_type i, nbpt = pararef.nb_points();
149 
150  for (i = 0; i < nbpt; ++i) {
151  gmm::clear(a);
152  for (dim_type n = 0; n < N; ++n)
153  gmm::add(gmm::scaled(ivect[n],pararef.points()[i][n]),a);
154  //a.addmul(pararef.points()[i][n], ivect[n]);
155  pararef.points()[i] = a;
156  barycenter += a;
157  }
158  barycenter /= double(nbpt);
159 
160  std::vector<size_type> tab1(N+1), tab(N), tab3(nbpt);
161  size_type total = 0;
162  std::fill(tab.begin(), tab.end(), 0);
163  while (tab[N-1] != iref[N-1]) {
164  for (a = org, i = 0; i < N; i++)
165  gmm::add(gmm::scaled(ivect[i], scalar_type(tab[i])),a);
166  //a.addmul(scalar_type(tab[i]), ivect[i]);
167 
168  for (i = 0; i < nbpt; i++)
169  tab3[i] = me.add_point(a + pararef.points()[i]);
170  size_type ib = me.add_point(a + barycenter);
171  me.add_pyramid(tab3[0],tab3[1],tab3[2],tab3[3],ib);
172  me.add_pyramid(tab3[4],tab3[6],tab3[5],tab3[7],ib);
173  me.add_pyramid(tab3[0],tab3[4],tab3[1],tab3[5],ib);
174  me.add_pyramid(tab3[1],tab3[5],tab3[3],tab3[7],ib);
175  me.add_pyramid(tab3[3],tab3[7],tab3[2],tab3[6],ib);
176  me.add_pyramid(tab3[2],tab3[6],tab3[0],tab3[4],ib);
177 
178  for (dim_type l = 0; l < N; l++) {
179  tab[l]++; total++;
180  if (l < N-1 && tab[l] >= iref[l]) { total -= tab[l]; tab[l] = 0; }
181  else break;
182  }
183  }
184  }
185 
186 
187  /* deformation inside a unit square -- ugly */
188  static base_node shake_func(const base_node& x) {
189  base_node z(x.size());
190  scalar_type c1 = 1., c2 = 1.;
191  for (size_type i=0; i < x.size(); ++i) {
192  c1*=(x[i]*(1.-x[i]));
193  c2*=(.5 - gmm::abs(x[i]-.5));
194  }
195  z[0] = x[0] + c1; // 1.
196  for (size_type i=1; i < x.size(); ++i) {
197  z[i] = x[i] + c2/3.; // /10
198  }
199  return z;
200  }
201 
202  static base_node radial_deformation(const base_node& x) {
203  GMM_ASSERT1(x.size() == 2, "For two-dimensional meshes only. \n");
204  base_node z(x.size());
205  z[0] = x[0] - 0.5 ;
206  z[1] = x[1] - 0.5 ;
207  scalar_type r = sqrt( z[0] * z[0] + z[1] * z[1] ) ;
208  scalar_type theta = atan2(z[1], z[0]);
209  if ( r < 0.5 - 1.e-6)
210 // theta += 1000. * gmm::sqrt(r) * (0.5 - r) * (0.5 - r) * (0.5 - r) * (0.5 - r) * gmm::sqrt(gmm::abs(0.1 - r)) * gmm::sqrt(gmm::abs(0.15 - r)) ;
211  theta += 10000. * gmm::sqrt(r) * (0.5 - r) * (0.5 - r) * (0.5 - r) * (0.5 - r) * (0.1 - r) * (0.15 - r) ;
212  z[0] = r * cos(theta) + 0.5;
213  z[1] = r * sin(theta) + 0.5;
214  return z;
215  }
216 
217  static void noise_unit_mesh(mesh& m, std::vector<size_type> nsubdiv,
219  size_type N = nsubdiv.size();
220  for (dal::bv_visitor ip(m.points().index()); !ip.finished(); ++ip) {
221  bool is_border = false;
222  base_node& P = m.points()[ip];
223  for (size_type i=0; i < N; ++i) {
224  if (gmm::abs(P[i]) < 1e-10 || gmm::abs(P[i]-1.) < 1e-10)
225  is_border = true;
226  }
227  if (!is_border) {
228  P = shake_func(P);
229  if (N == 2) P = radial_deformation(P) ;
230  for (size_type i=0; i < N; ++i)
231  P[i] += 0.*(double(1)/double(nsubdiv[i]* pgt->complexity()))
232  * gmm::random(double());
233  }
234  }
235  }
236 
237 
238  void regular_unit_mesh(mesh& m, std::vector<size_type> nsubdiv,
239  bgeot::pgeometric_trans pgt, bool noised) {
240  mesh msh;
241  dim_type N = dim_type(nsubdiv.size());
242  base_node org(N); gmm::clear(org);
243  std::vector<base_small_vector> vtab(N);
244  for (dim_type i = 0; i < N; i++) {
245  vtab[i] = base_small_vector(N); gmm::clear(vtab[i]);
246  (vtab[i])[i] = 1. / scalar_type(nsubdiv[i]) * 1.;
247  }
248  if (pgt->basic_structure() == bgeot::simplex_structure(N)) {
249  getfem::parallelepiped_regular_simplex_mesh
250  (msh, N, org, vtab.begin(), nsubdiv.begin());
251  } else if (pgt->basic_structure() == bgeot::parallelepiped_structure(N)) {
252  getfem::parallelepiped_regular_mesh
253  (msh, N, org, vtab.begin(), nsubdiv.begin());
254  } else if (pgt->basic_structure() == bgeot::prism_P1_structure(N)) {
255  getfem::parallelepiped_regular_prism_mesh
256  (msh, N, org, vtab.begin(), nsubdiv.begin());
257  } else if (pgt->basic_structure() == bgeot::pyramid_QK_structure(1)) {
258  getfem::parallelepiped_regular_pyramid_mesh
259  (msh, org, vtab.begin(), nsubdiv.begin());
260  } else {
261  GMM_ASSERT1(false, "cannot build a regular mesh for "
263  }
264 
265  m.clear();
266  /* build a mesh with a geotrans of degree K */
267  for (dal::bv_visitor cv(msh.convex_index()); !cv.finished(); ++cv) {
268  if (pgt == msh.trans_of_convex(cv)) {
269  m.add_convex_by_points(msh.trans_of_convex(cv),
270  msh.points_of_convex(cv).begin());
271  } else {
272  std::vector<base_node> pts(pgt->nb_points());
273  for (size_type i=0; i < pgt->nb_points(); ++i) {
274  pts[i] = msh.trans_of_convex(cv)->transform
275  (pgt->convex_ref()->points()[i], msh.points_of_convex(cv));
276  }
277  m.add_convex_by_points(pgt, pts.begin());
278  }
279  }
280 
281  /* apply a continuous deformation + some noise */
282  if (noised) noise_unit_mesh(m, nsubdiv, pgt);
283 
284  m.optimize_structure(false);
285  }
286 
287 
288 
289  void regular_mesh(mesh& m, const std::string &st) {
290  std::stringstream s(st);
291  bgeot::md_param PARAM;
292  PARAM.read_param_file(s);
293 
294  std::string GT = PARAM.string_value("GT");
295  GMM_ASSERT1(!GT.empty(), "regular mesh : you have at least to "
296  "specify the geometric transformation");
298 
299  size_type N = pgt->dim();
300  base_small_vector org(N); gmm::clear(org);
301 
302  const auto &o = PARAM.array_value("ORG");
303  if (o.size() > 0) {
304  GMM_ASSERT1(o.size() == N,
305  "ORG parameter should be an array of size " << N);
306  for (size_type i = 0; i < N; ++i) {
307  GMM_ASSERT1(o[i].type_of_param() == bgeot::md_param::REAL_VALUE,
308  "ORG should be a real array.");
309  org[i] = o[i].real();
310  }
311  }
312 
313  bool noised = (PARAM.int_value("NOISED") != 0);
314 
315  std::vector<size_type> nsubdiv(N);
316  gmm::fill(nsubdiv, 2);
317  const auto &ns = PARAM.array_value("NSUBDIV");
318  if (ns.size() > 0) {
319  GMM_ASSERT1(ns.size() == N,
320  "NSUBDIV parameter should be an array of size " << N);
321  for (size_type i = 0; i < N; ++i) {
322  GMM_ASSERT1(ns[i].type_of_param() == bgeot::md_param::REAL_VALUE,
323  "NSUBDIV should be an integer array");
324  nsubdiv[i] = size_type(ns[i].real()+0.5);
325  }
326  }
327 
328  base_small_vector sizes(N);
329  gmm::fill(sizes, 1.0);
330 
331  const auto &si = PARAM.array_value("SIZES");
332  if (si.size() > 0) {
333  GMM_ASSERT1(si.size() == N,
334  "SIZES parameter should be an array of size " << N);
335  for (size_type i = 0; i < N; ++i) {
336  GMM_ASSERT1(si[i].type_of_param() == bgeot::md_param::REAL_VALUE,
337  "SIZES should be a real array");
338  sizes[i] = si[i].real();
339  }
340  }
341 
342  regular_unit_mesh(m, nsubdiv, pgt, noised);
343 
344  base_matrix M(N,N);
345  for (size_type i=0; i < N; ++i) M(i,i) = sizes[i];
346  m.transformation(M);
347  m.translation(org);
348 
349  }
350 
351 
352  void regular_ball_mesh(mesh& m, const std::string &st) {
353  std::stringstream s(st);
354  bgeot::md_param PARAM;
355  PARAM.read_param_file(s);
356 
357  std::string GT = PARAM.string_value("GT");
358  GMM_ASSERT1(!GT.empty(), "regular ball mesh : you have at least to "
359  "specify the geometric transformation");
361 
362  size_type N = pgt->dim();
363  base_small_vector org(N);
364 
365  const auto &o = PARAM.array_value("ORG");
366  if (o.size() > 0) {
367  GMM_ASSERT1(o.size() == N,
368  "ORG parameter should be an array of size " << N);
369  for (size_type i = 0; i < N; ++i) {
370  GMM_ASSERT1(o[i].type_of_param() == bgeot::md_param::REAL_VALUE,
371  "ORG should be a real array");
372  org[i] = o[i].real();
373  }
374  }
375 
376  // "NOISED" applies only to the interior of all merged sub-meshes
377  bool noised = (PARAM.int_value("NOISED") != 0);
378 
379  size_type nsubdiv0(2), nsubdiv1(2);
380  const auto &ns = PARAM.array_value("NSUBDIV");
381  if (ns.size() > 0) {
382  GMM_ASSERT1(ns.size() == 2,
383  "NSUBDIV parameter should be an array of size 2");
384  for (size_type i = 0; i < 2; ++i)
385  GMM_ASSERT1(ns[i].type_of_param() == bgeot::md_param::REAL_VALUE,
386  "NSUBDIV should be an integer array");
387  nsubdiv0 = size_type(ns[0].real()+0.5);
388  nsubdiv1 = size_type(ns[1].real()+0.5);
389  }
390 
391  scalar_type radius(1), core_ratio(M_SQRT1_2);
392  const auto &si = PARAM.array_value("SIZES");
393  if (si.size() > 0) {
394  GMM_ASSERT1(si.size() == 1,
395  "SIZES parameter should be an array of size 1");
396  GMM_ASSERT1(si[0].type_of_param() == bgeot::md_param::REAL_VALUE,
397  "SIZES should be a real array");
398  radius = si[0].real();
399  }
400 
401  std::vector<size_type> nsubdiv(N);
402  gmm::fill(nsubdiv, nsubdiv0);
403  regular_unit_mesh(m, nsubdiv, pgt, noised);
404  std::vector<mesh> mm(N);
405  for (size_type i = 0; i < N; ++i) {
406  gmm::fill(nsubdiv, nsubdiv0);
407  nsubdiv[i] = nsubdiv1;
408  regular_unit_mesh(mm[i], nsubdiv, pgt, noised);
409  }
410 
411  base_matrix M(N,N);
412  for (size_type i=0; i < N; ++i) M(i,i) = core_ratio;
413  m.transformation(M);
414 
415  scalar_type peel_ratio = scalar_type(1)-core_ratio;
416  for (size_type i = 0; i < N; ++i) {
417  base_matrix MM(M);
418  MM(i,i) = peel_ratio;
419  mm[i].transformation(MM);
420 
421  std::vector<base_node> pts(mm[i].points().card(), base_node(N));
422  size_type j(0);
423  for (dal::bv_visitor pt(mm[i].points().index()); !pt.finished(); ++pt, ++j) {
424  pts[j] = mm[i].points()[pt];
425  for (unsigned k=0; k < N; ++k)
426  if (k != i) pts[j][k] += (pts[j][k]/core_ratio) * pts[j][i];
427  }
428  j = size_type(0);
429  for (dal::bv_visitor pt(mm[i].points().index()); !pt.finished(); ++pt, ++j)
430  mm[i].points()[pt] = pts[j];
431 
432  base_small_vector trsl(N);
433  trsl[i] = core_ratio;
434  mm[i].translation(trsl);
435  for (dal::bv_visitor cv(mm[i].convex_index()); !cv.finished(); ++cv)
436  m.add_convex_by_points(mm[i].trans_of_convex(cv),
437  mm[i].points_of_convex(cv).begin());
438  }
439 
440  std::vector<base_node> pts(m.points().card(), base_node(N));
441  size_type j(0);
442  for (dal::bv_visitor pt(m.points().index()); !pt.finished(); ++pt, ++j) {
443  pts[j] = m.points()[pt];
444  scalar_type maxcoord(0);
445  size_type kmax(0);
446  for (size_type k=0; k < N; ++k)
447  if (gmm::abs(pts[j][k]) > maxcoord) {
448  maxcoord = gmm::abs(pts[j][k]);
449  kmax = k;
450  }
451  if (maxcoord > 1e-10) {
452  scalar_type l(0), l0(0);
453  for (size_type k=0; k < N; ++k) {
454  if (k != kmax) {
455  scalar_type theta = M_PI_4 * pts[j][k] / maxcoord;
456  scalar_type c0 = std::min(scalar_type(1), maxcoord);
457  pts[j][k] = c0*tan(theta)* maxcoord + (scalar_type(1)-c0)*pts[j][k];
458  }
459  l += pts[j][k] * pts[j][k];
460  }
461  l = sqrt(l);
462  l0 = l/maxcoord;
463  scalar_type scale(radius);
464  scalar_type c(core_ratio);
465  c *= std::max(scalar_type(0.3),
466  (scalar_type(1) - sqrt(l0*l0 - scalar_type(1))));
467  if (l > c) {
468  scale -= (l - c) / (l0 - c) * (radius - radius/(l/maxcoord));
469  }
470  for (size_type k=0; k < N; ++k)
471  pts[j][k] *= scale;
472  }
473  }
474  j = size_type(0);
475  for (dal::bv_visitor pt(m.points().index()); !pt.finished(); ++pt, ++j)
476  m.points()[pt] = pts[j];
477  m.points().resort();
478 
479  size_type symmetries(PARAM.int_value("SYMMETRIES"));
480  symmetries = std::min(symmetries,N);
481 
482  for (size_type sym=0; sym < N-symmetries; ++sym) {
483  size_type sym0 = (sym+1) % N;
484  if (sym0 != 0) {
485  gmm::clear(M);
486  M(sym,sym0) = scalar_type(-1);
487  M(sym0,sym) = scalar_type(1);
488  for (size_type i=0; i < N; ++i)
489  if (i != sym && i != sym0) M(i,i) = scalar_type(1);
490  } else {
491  base_matrix M1(M), M2(M);
492  gmm::mult(M1,M2,M);
493  }
494  mesh m0;
495  m0.copy_from(m);
496  m0.transformation(M);
497  for (dal::bv_visitor cv(m0.convex_index()); !cv.finished(); ++cv)
498  m.add_convex_by_points(m0.trans_of_convex(cv),
499  m0.points_of_convex(cv).begin());
500  }
501 
502  m.translation(org);
503 
504  }
505  void regular_ball_shell_mesh(mesh &m, const std::string &st) {
506  std::stringstream s(st);
507  bgeot::md_param PARAM;
508  PARAM.read_param_file(s);
509 
510  std::string GT = PARAM.string_value("GT");
511  GMM_ASSERT1(!GT.empty(), "regular ball mesh : you have at least to "
512  "specify the geometric transformation");
514 
515  size_type N = pgt->dim();
516  base_small_vector org(N);
517  const auto &o = PARAM.array_value("ORG");
518  if (o.size() > 0) {
519  GMM_ASSERT1(o.size() == N,
520  "ORG parameter should be an array of size " << N);
521  for (size_type i = 0; i < N; ++i) {
522  GMM_ASSERT1(o[i].type_of_param() == bgeot::md_param::REAL_VALUE,
523  "ORG should be a real array");
524  org[i] = o[i].real();
525  }
526  }
527 
528  // "NOISED" applies only to the interior of all merged sub-meshes
529  bool noised = (PARAM.int_value("NOISED") != 0);
530 
531  size_type nsubdiv0(3), nsubdiv1(2);
532  const auto &ns = PARAM.array_value("NSUBDIV");
533  if (ns.size() > 0) {
534  GMM_ASSERT1(ns.size() == 2,
535  "NSUBDIV parameter should be an array of size 2");
536  for (size_type i = 0; i < 2; ++i)
537  GMM_ASSERT1(ns[i].type_of_param() == bgeot::md_param::REAL_VALUE,
538  "NSUBDIV should be an integer array");
539  nsubdiv0 = size_type(ns[0].real()+0.5);
540  nsubdiv1 = size_type(ns[1].real()+0.5);
541  }
542 
543  scalar_type radius(1), thickness(0.5);
544  const auto &si = PARAM.array_value("SIZES");
545  if (si.size() > 0) {
546  GMM_ASSERT1(si.size() == 2,
547  "SIZES parameter should be an array of size 2");
548  GMM_ASSERT1(si[0].type_of_param() == bgeot::md_param::REAL_VALUE &&
549  si[1].type_of_param() == bgeot::md_param::REAL_VALUE ,
550  "SIZES should be a real array");
551  radius = si[0].real();
552  thickness = si[1].real();
553  }
554 
555  std::vector<size_type> nsubdiv(N, nsubdiv0);
556  nsubdiv[N-1] = nsubdiv1;
557 
558  mesh m0;
559  regular_unit_mesh(m0, nsubdiv, pgt, noised);
560 
561  std::vector<base_node> pts(m0.points().card(), base_node(N));
562  size_type j(0);
563  for (dal::bv_visitor pt(m0.points().index()); !pt.finished(); ++pt, ++j) {
564  pts[j] = m0.points()[pt];
565  scalar_type l(1);
566  for (size_type k=0; k < N-1; ++k) {
567  pts[j][k] = tan(pts[j][k]*M_PI_4);
568  l += pts[j][k]*pts[j][k];
569  }
570  l = sqrt(l);
571  scalar_type r(radius-thickness+thickness*pts[j][N-1]);
572  for (size_type k=0; k < N-1; ++k)
573  pts[j][k] *= r/l;
574  pts[j][N-1] = r/l;
575  }
576 
577  j = size_type(0);
578  for (dal::bv_visitor pt(m0.points().index()); !pt.finished(); ++pt, ++j)
579  m0.points()[pt] = pts[j];
580  m0.points().resort();
581 
582  base_matrix M(N,N);
583  for (size_type i=0; i < N-1; ++i)
584  M(i,i+1) = scalar_type(1);
585  M(N-1,0) = scalar_type(1);
586 
587  for (size_type i = 0; i < N; ++i) {
588  m0.transformation(M);
589  for (dal::bv_visitor cv(m0.convex_index()); !cv.finished(); ++cv)
590  m.add_convex_by_points(m0.trans_of_convex(cv),
591  m0.points_of_convex(cv).begin());
592  }
593 
594  size_type symmetries(PARAM.int_value("SYMMETRIES"));
595  symmetries = std::min(symmetries,N);
596 
597  for (size_type sym=0; sym < N-symmetries; ++sym) {
598  size_type sym0 = (sym+1) % N;
599  if (sym0 != 0) {
600  gmm::clear(M);
601  M(sym,sym0) = scalar_type(-1);
602  M(sym0,sym) = scalar_type(1);
603  for (size_type i=0; i < N; ++i)
604  if (i != sym && i != sym0) M(i,i) = scalar_type(1);
605  } else {
606  base_matrix M1(M), M2(M);
607  gmm::mult(M1,M2,M);
608  }
609  m0.copy_from(m);
610  m0.transformation(M);
611  for (dal::bv_visitor cv(m0.convex_index()); !cv.finished(); ++cv)
612  m.add_convex_by_points(m0.trans_of_convex(cv),
613  m0.points_of_convex(cv).begin());
614  }
615 
616  m.translation(org);
617 
618  }
619 
620 
621 } /* end of namespace getfem. */
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1173
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::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem_regular_meshes.h
Build regular meshes.
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
getfem::mesh::optimize_structure
void optimize_structure(bool with_renumbering=true)
Pack the mesh : renumber convexes and nodes such that there is no holes in their numbering.
Definition: getfem_mesh.cc:223
getfem::regular_ball_mesh
void regular_ball_mesh(mesh &m, const std::string &st)
Build a regular mesh on a ball, parametrized by the string st.
Definition: getfem_regular_meshes.cc:352
bgeot::convex< base_node >
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::mesh_structure::ind_points_of_convex
const ind_set & ind_points_of_convex(size_type ic) const
Return a container to the list of points attached to convex ic.
Definition: bgeot_mesh_structure.h:107
getfem::regular_ball_shell_mesh
void regular_ball_shell_mesh(mesh &m, const std::string &st)
Build a regular mesh on a ball shell, parametrized by the string st.
Definition: getfem_regular_meshes.cc:505
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
getfem::regular_mesh
void regular_mesh(mesh &m, const std::string &st)
Build a regular mesh parametrized by the string st.
Definition: getfem_regular_meshes.cc:289
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:91
getfem::mesh::copy_from
void copy_from(const mesh &m)
Clone a mesh.
Definition: getfem_mesh.cc:487
bgeot::mesh_structure::nb_convex
size_type nb_convex() const
The total number of convexes in the mesh.
Definition: bgeot_mesh_structure.h:96
getfem::mesh::transformation
void transformation(const base_matrix &)
apply the given matrix transformation to each mesh node.
Definition: getfem_mesh.cc:255
getfem::mesh::translation
void translation(const base_small_vector &)
Apply the given translation to each mesh node.
Definition: getfem_mesh.cc:252
getfem::mesh::add_convex_by_points
size_type add_convex_by_points(bgeot::pgeometric_trans pgt, ITER ipts, const scalar_type tol=scalar_type(0))
Add a convex to the mesh, given a geometric transformation and a list of point coordinates.
Definition: getfem_mesh.h:558
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:747
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:73
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
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::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::regular_unit_mesh
void regular_unit_mesh(mesh &m, std::vector< size_type > nsubdiv, bgeot::pgeometric_trans pgt, bool noised=false)
Build a regular mesh of the unit square/cube/, etc.
Definition: getfem_regular_meshes.cc:238
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