GetFEM  5.5
bgeot_convex_ref.cc
1 /*===========================================================================
2 
3  Copyright (C) 2001-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 #include "getfem/dal_singleton.h"
22 #include "getfem/bgeot_convex_ref.h"
23 #include "getfem/bgeot_mesh_structure.h"
24 #include "getfem/bgeot_comma_init.h"
25 
26 namespace bgeot {
27 
28  // ******************************************************************
29  // Interface with qhull
30  // ******************************************************************
31 
32 # if !defined(GETFEM_HAVE_LIBQHULL_R_QHULL_RA_H)
33  void qhull_delaunay(const std::vector<base_node> &,
34  gmm::dense_matrix<size_type>&) {
35  GMM_ASSERT1(false, "Qhull header files not installed. "
36  "Install qhull library and reinstall GetFEM library.");
37  }
38 # else
39 
40  extern "C" { // old versions of this header lack a __cplusplus guard
41 # include <libqhull_r/qhull_ra.h>
42  }
43 
44  void qhull_delaunay(const std::vector<base_node> &pts,
45  gmm::dense_matrix<size_type>& simplexes) {
46  // cout << "running delaunay with " << pts.size() << " points\n";
47  size_type dim = pts[0].size(); /* points dimension. */
48  if (pts.size() <= dim) { gmm::resize(simplexes, dim+1, 0); return; }
49  if (pts.size() == dim+1) {
50  gmm::resize(simplexes, dim+1, 1);
51  for (size_type i=0; i <= dim; ++i) simplexes(i, 0) = i;
52  return;
53  }
54  std::vector<coordT> Pts(dim * pts.size());
55  for (size_type i=0; i < pts.size(); ++i)
56  gmm::copy(pts[i], gmm::sub_vector(Pts, gmm::sub_interval(i*dim, dim)));
57  boolT ismalloc=0; /* True if qhull should free points in
58  * qh_freeqhull() or reallocation */
59  /* Be Aware: option QJ could destabilizate all, it can break everything. */
60  /* option Qbb -> QbB (????) */
61  /* option flags for qhull, see qh_opt.htm */
62  char flags[]= "qhull QJ d Qbb Pp T0"; //QJ s i TO";//"qhull Tv";
63  FILE *outfile= 0; /* output from qh_produce_output()
64  * use NULL to skip qh_produce_output() */
65  FILE *errfile= stderr; /* error messages from qhull code */
66  int exitcode; /* 0 if no error from qhull */
67  facetT *facet; /* set by FORALLfacets */
68  int curlong, totlong; /* memory remaining after qh_memfreeshort */
69  vertexT *vertex, **vertexp;
70  qhT context = {};
71  qhT* qh = &context;
72  exitcode = qh_new_qhull (qh, int(dim), int(pts.size()), &Pts[0], ismalloc,
73  flags, outfile, errfile);
74  if (!exitcode) { /* if no error */
75  size_type nbf=0;
76  FORALLfacets { if (!facet->upperdelaunay) nbf++; }
77  gmm::resize(simplexes, dim+1, nbf);
78  /* 'qh facet_list' contains the convex hull */
79  nbf=0;
80  FORALLfacets {
81  if (!facet->upperdelaunay) {
82  size_type s=0;
83  FOREACHvertex_(facet->vertices) {
84  assert(s < (unsigned)(dim+1));
85  simplexes(s++,nbf) = qh_pointid(qh, vertex->point);
86  }
87  nbf++;
88  }
89  }
90  }
91  qh_freeqhull(qh, !qh_ALL);
92  qh_memfreeshort(qh, &curlong, &totlong);
93  if (curlong || totlong)
94  cerr << "qhull internal warning (main): did not free " << totlong <<
95  " bytes of long memory (" << curlong << " pieces)\n";
96 
97  }
98 
99 #endif
100 
101 
102 
103  size_type simplexified_tab(pconvex_structure cvs, size_type **tab);
104 
105  static void simplexify_convex(bgeot::convex_of_reference *cvr,
106  mesh_structure &m) {
107  pconvex_structure cvs = cvr->structure();
108  m.clear();
109  auto basic_cvs = basic_structure(cvs);
110  dim_type n = basic_cvs->dim();
111  std::vector<size_type> ipts(n+1);
112  if (basic_cvs->nb_points() == n + 1) {
113  for (size_type i = 0; i <= n; ++i) ipts[i] = i;
114  m.add_simplex(n, ipts.begin());
115  }
116  else {
117  size_type *tab;
118  size_type nb = simplexified_tab(basic_cvs, &tab);
119  if (nb) {
120  for (size_type nc = 0; nc < nb; ++nc) {
121  for (size_type i = 0; i <= n; ++i) ipts[i] = *tab++;
122  m.add_simplex(n, ipts.begin());
123  }
124  } else {
125 # if defined(GETFEM_HAVE_LIBQHULL_R_QHULL_RA_H)
126  gmm::dense_matrix<size_type> t;
127  qhull_delaunay(cvr->points(), t);
128  for (size_type nc = 0; nc < gmm::mat_ncols(t); ++nc) {
129  for (size_type i = 0; i <= n; ++i) ipts[i] = t(i,nc);
130  m.add_simplex(n, ipts.begin());
131  }
132 # endif
133  }
134  }
135  }
136 
137  /* ********************************************************************* */
138  /* Point tab storage. */
139  /* ********************************************************************* */
140 
141  struct stored_point_tab_key : virtual public dal::static_stored_object_key {
142  const stored_point_tab *pspt;
143  bool compare(const static_stored_object_key &oo) const override {
144  const stored_point_tab_key &o
145  = dynamic_cast<const stored_point_tab_key &>(oo);
146  const stored_point_tab &x = *pspt;
147  const stored_point_tab &y = *(o.pspt);
148  if (&x == &y) return false;
149  std::vector<base_node>::const_iterator it1 = x.begin(), it2 = y.begin();
150  base_node::const_iterator itn1, itn2, itne;
151  for ( ; it1 != x.end() && it2 != y.end() ; ++it1, ++it2) {
152  if ((*it1).size() < (*it2).size()) return true;
153  if ((*it1).size() > (*it2).size()) return false;
154  itn1 = (*it1).begin(); itne = (*it1).end(); itn2 = (*it2).begin();
155  for ( ; itn1 != itne ; ++itn1, ++itn2)
156  if (*itn1 < *itn2) return true;
157  else if (*itn1 > *itn2) return false;
158  }
159  if (it2 != y.end()) return true;
160  return false;
161  }
162  bool equal(const static_stored_object_key &oo) const override {
163  auto &o = dynamic_cast<const stored_point_tab_key &>(oo);
164  auto &x = *pspt;
165  auto &y = *(o.pspt);
166  if (&x == &y) return true;
167  if (x.size() != y.size()) return false;
168  auto it1 = x.begin();
169  auto it2 = y.begin();
170  base_node::const_iterator itn1, itn2, itne;
171  for ( ; it1 != x.end() && it2 != y.end() ; ++it1, ++it2) {
172  if ((*it1).size() != (*it2).size()) return false;
173  itn1 = (*it1).begin(); itne = (*it1).end(); itn2 = (*it2).begin();
174  for ( ; itn1 != itne ; ++itn1, ++itn2)
175  if (*itn1 != *itn2) return false;
176  }
177  return true;
178  }
179  stored_point_tab_key(const stored_point_tab *p) : pspt(p) {}
180  };
181 
182 
183  pstored_point_tab store_point_tab(const stored_point_tab &spt) {
184  dal::pstatic_stored_object_key
185  pk = std::make_shared<stored_point_tab_key>(&spt);
186  dal::pstatic_stored_object o = dal::search_stored_object(pk);
187  if (o) return std::dynamic_pointer_cast<const stored_point_tab>(o);
188  pstored_point_tab p = std::make_shared<stored_point_tab>(spt);
189  DAL_STORED_OBJECT_DEBUG_CREATED(p.get(), "Stored point tab");
190  dal::pstatic_stored_object_key
191  psp = std::make_shared<stored_point_tab_key>(p.get());
192  dal::add_stored_object(psp, p, dal::AUTODELETE_STATIC_OBJECT);
193  return p;
194  }
195 
196  convex_of_reference::convex_of_reference
197  (pconvex_structure cvs_, bool auto_basic_) :
198  convex<base_node>(move(cvs_)), basic_convex_ref_(0),
199  auto_basic(auto_basic_) {
200  DAL_STORED_OBJECT_DEBUG_CREATED(this, "convex of refrence");
201  psimplexified_convex = std::make_shared<mesh_structure>();
202  // dal::singleton<cleanup_simplexified_convexes>::instance()
203  // .push_back(psimplexified_convex);
204  }
205 
206  /* should be called on the basic_convex_ref */
207  const mesh_structure* convex_of_reference::simplexified_convex() const {
208  GMM_ASSERT1(auto_basic,
209  "always use simplexified_convex on the basic_convex_ref() "
210  "[this=" << nb_points() << ", basic="
211  << basic_convex_ref_->nb_points());
212  return psimplexified_convex.get();
213  }
214 
215  // Key type for static storing
216  class convex_of_reference_key : virtual public dal::static_stored_object_key{
217  int type; // 0 = simplex structure degree K
218  // 1 = equilateral simplex of ref.
219  // 2 = dummy
220  dim_type N; short_type K; short_type nf;
221  public :
222  bool compare(const static_stored_object_key &oo) const override{
223  const convex_of_reference_key &o
224  = dynamic_cast<const convex_of_reference_key &>(oo);
225  if (type < o.type) return true;
226  if (type > o.type) return false;
227  if (N < o.N) return true;
228  if (N > o.N) return false;
229  if (K < o.K) return true;
230  if (K > o.K) return false;
231  if (nf < o.nf) return true;
232  return false;
233  }
234  bool equal(const static_stored_object_key &oo) const override{
235  auto &o = dynamic_cast<const convex_of_reference_key &>(oo);
236  if (type != o.type) return false;
237  if (N != o.N) return false;
238  if (K != o.K) return false;
239  if (nf != o.nf) return false;
240  return true;
241  }
242  convex_of_reference_key(int t, dim_type NN, short_type KK = 0,
243  short_type nnf = 0)
244  : type(t), N(NN), K(KK), nf(nnf) {}
245  };
246 
247 
248  /* simplexes. */
249 
250  class K_simplex_of_ref_ : public convex_of_reference {
251  public :
252  scalar_type is_in(const base_node &pt) const {
253  // return a negative number if pt is in the convex
254  GMM_ASSERT1(pt.size() == cvs->dim(),
255  "K_simplex_of_ref_::is_in: Dimensions mismatch");
256  scalar_type e = -1.0, r = (pt.size() > 0) ? -pt[0] : 0.0;
257  base_node::const_iterator it = pt.begin(), ite = pt.end();
258  for (; it != ite; e += *it, ++it) r = std::max(r, -(*it));
259  return std::max(r, e);
260  }
261 
262  scalar_type is_in_face(short_type f, const base_node &pt) const {
263  // return zero if pt is in the face of the convex
264  // negative if the point is on the side of the face where the element is
265  GMM_ASSERT1(pt.size() == cvs->dim(),
266  "K_simplex_of_ref_::is_in_face: Dimensions mismatch");
267  if (f > 0) return -pt[f-1];
268  scalar_type e = -1.0;
269  base_node::const_iterator it = pt.begin(), ite = pt.end();
270  for (; it != ite; e += *it, ++it) {};
271  return e / sqrt(scalar_type(pt.size()));
272  }
273 
274  void project_into(base_node &pt) const {
275  if (auto_basic) {
276  GMM_ASSERT1(pt.size() == cvs->dim(),
277  "K_simplex_of_ref_::project_into: Dimensions mismatch");
278  scalar_type sum_coordinates = 0.0;
279  for (const auto &coord : pt) sum_coordinates += coord;
280  if (sum_coordinates > 1.0) gmm::scale(pt, 1.0 / sum_coordinates);
281  for (auto &coord : pt) {
282  if (coord < 0.0) coord = 0.0;
283  if (coord > 1.0) coord = 1.0;
284  }
285  } else
286  basic_convex_ref_->project_into(pt);
287  }
288 
289  K_simplex_of_ref_(dim_type NN, short_type KK) :
290  convex_of_reference(simplex_structure(NN, KK), (KK == 1) || (NN == 0))
291  {
292  if ((KK != 1) && (NN != 0)) basic_convex_ref_ = simplex_of_reference(NN, 1);
293  size_type R = cvs->nb_points();
294  convex<base_node>::points().resize(R);
295  normals_.resize(NN+1);
296  base_node null(NN); null.fill(0.0);
297  std::fill(normals_.begin(), normals_.end(), null);
298  std::fill(convex<base_node>::points().begin(),
299  convex<base_node>::points().end(), null);
300  for (size_type i = 1; i <= NN; ++i)
301  normals_[i][i-1] = -1.0;
302  if (NN > 0)
303  std::fill(normals_[0].begin(), normals_[0].end(),
304  scalar_type(1.0)/sqrt(scalar_type(NN)));
305  base_node c(NN); c.fill(0.0);
306 
307  if (KK == 0) {
308  c.fill(1.0/(NN+1));
309  convex<base_node>::points()[0] = c;
310  }
311  else {
312  size_type sum = 0, l;
313  for (size_type r = 0; r < R; ++r) {
314  convex<base_node>::points()[r] = c;
315  if (KK != 0 && NN > 0) {
316  l = 0; c[l] += 1.0 / scalar_type(KK); sum++;
317  while (sum > KK) {
318  sum -= int(floor(0.5+(c[l] * KK)));
319  c[l] = 0.0; l++; if (l == NN) break;
320  c[l] += 1.0 / scalar_type(KK); sum++;
321  }
322  }
323  }
324  }
325  ppoints = store_point_tab(convex<base_node>::points());
326  if (auto_basic) simplexify_convex(this, *psimplexified_convex);
327  }
328  };
329 
330  pconvex_ref simplex_of_reference(dim_type nc, short_type K) {
331  dal::pstatic_stored_object_key
332  pk = std::make_shared<convex_of_reference_key>(0, nc, K);
333  dal::pstatic_stored_object o = dal::search_stored_object(pk);
334  if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
335  pconvex_ref p = std::make_shared<K_simplex_of_ref_>(nc, K);
336  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
337  dal::PERMANENT_STATIC_OBJECT);
338  pconvex_ref p1 = basic_convex_ref(p);
339  if (p != p1) add_dependency(p, p1);
340  return p;
341  }
342 
343  /* ******************************************************************** */
344  /* Incomplete Q2 quadrilateral or hexahedral of reference. */
345  /* ******************************************************************** */
346  /* By Yao Koutsawa <[email protected]> 2012-12-10 */
347 
348  class Q2_incomplete_of_ref_ : public convex_of_reference {
349  public :
350  scalar_type is_in(const base_node& pt) const
351  { return basic_convex_ref_->is_in(pt); }
352  scalar_type is_in_face(short_type f, const base_node& pt) const
353  { return basic_convex_ref_->is_in_face(f, pt); }
354 
355  Q2_incomplete_of_ref_(dim_type nc) :
356  convex_of_reference(Q2_incomplete_structure(nc), false)
357  {
358  GMM_ASSERT1(nc == 2 || nc == 3, "Sorry exist only in dimension 2 or 3");
359  convex<base_node>::points().resize(cvs->nb_points());
360  normals_.resize(nc == 2 ? 4: 6);
361  basic_convex_ref_ = parallelepiped_of_reference(nc);
362 
363  if(nc==2) {
364  normals_[0] = { 1, 0};
365  normals_[1] = {-1, 0};
366  normals_[2] = { 0, 1};
367  normals_[3] = { 0,-1};
368 
369  convex<base_node>::points()[0] = base_node(0.0, 0.0);
370  convex<base_node>::points()[1] = base_node(0.5, 0.0);
371  convex<base_node>::points()[2] = base_node(1.0, 0.0);
372  convex<base_node>::points()[3] = base_node(0.0, 0.5);
373  convex<base_node>::points()[4] = base_node(1.0, 0.5);
374  convex<base_node>::points()[5] = base_node(0.0, 1.0);
375  convex<base_node>::points()[6] = base_node(0.5, 1.0);
376  convex<base_node>::points()[7] = base_node(1.0, 1.0);
377 
378  } else {
379  normals_[0] = { 1, 0, 0};
380  normals_[1] = {-1, 0, 0};
381  normals_[2] = { 0, 1, 0};
382  normals_[3] = { 0,-1, 0};
383  normals_[4] = { 0, 0, 1};
384  normals_[5] = { 0, 0,-1};
385 
386  convex<base_node>::points()[0] = base_node(0.0, 0.0, 0.0);
387  convex<base_node>::points()[1] = base_node(0.5, 0.0, 0.0);
388  convex<base_node>::points()[2] = base_node(1.0, 0.0, 0.0);
389  convex<base_node>::points()[3] = base_node(0.0, 0.5, 0.0);
390  convex<base_node>::points()[4] = base_node(1.0, 0.5, 0.0);
391  convex<base_node>::points()[5] = base_node(0.0, 1.0, 0.0);
392  convex<base_node>::points()[6] = base_node(0.5, 1.0, 0.0);
393  convex<base_node>::points()[7] = base_node(1.0, 1.0, 0.0);
394 
395  convex<base_node>::points()[8] = base_node(0.0, 0.0, 0.5);
396  convex<base_node>::points()[9] = base_node(1.0, 0.0, 0.5);
397  convex<base_node>::points()[10] = base_node(0.0, 1.0, 0.5);
398  convex<base_node>::points()[11] = base_node(1.0, 1.0, 0.5);
399 
400  convex<base_node>::points()[12] = base_node(0.0, 0.0, 1.0);
401  convex<base_node>::points()[13] = base_node(0.5, 0.0, 1.0);
402  convex<base_node>::points()[14] = base_node(1.0, 0.0, 1.0);
403  convex<base_node>::points()[15] = base_node(0.0, 0.5, 1.0);
404  convex<base_node>::points()[16] = base_node(1.0, 0.5, 1.0);
405  convex<base_node>::points()[17] = base_node(0.0, 1.0, 1.0);
406  convex<base_node>::points()[18] = base_node(0.5, 1.0, 1.0);
407  convex<base_node>::points()[19] = base_node(1.0, 1.0, 1.0);
408  }
409  ppoints = store_point_tab(convex<base_node>::points());
410  }
411  };
412 
413 
414  DAL_SIMPLE_KEY(Q2_incomplete_of_reference_key_, dim_type);
415 
416  pconvex_ref Q2_incomplete_of_reference(dim_type nc) {
417  dal::pstatic_stored_object_key
418  pk = std::make_shared<Q2_incomplete_of_reference_key_>(nc);
419  dal::pstatic_stored_object o = dal::search_stored_object(pk);
420  if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
421  pconvex_ref p = std::make_shared<Q2_incomplete_of_ref_>(nc);
422  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
423  dal::PERMANENT_STATIC_OBJECT);
424  pconvex_ref p1 = basic_convex_ref(p);
425  if (p != p1) add_dependency(p, p1);
426  return p;
427  }
428 
429  /* ******************************************************************** */
430  /* Pyramidal element of reference. */
431  /* ******************************************************************** */
432 
433  class pyramid_QK_of_ref_ : public convex_of_reference {
434  public :
435  scalar_type is_in_face(short_type f, const base_node& pt) const {
436  // return zero if pt is in the face of the convex
437  // negative if the point is on the side of the face where the element is
438  GMM_ASSERT1(pt.size() == 3, "Dimensions mismatch");
439  if (f == 0)
440  return -pt[2];
441  else
442  return gmm::vect_sp(normals_[f], pt) - sqrt(2.)/2.;
443  }
444 
445  scalar_type is_in(const base_node& pt) const {
446  // return a negative number if pt is in the convex
447  scalar_type r = is_in_face(0, pt);
448  for (short_type f = 1; f < 5; ++f) r = std::max(r, is_in_face(f, pt));
449  return r;
450  }
451 
452  void project_into(base_node &pt) const {
453  if (auto_basic) {
454  GMM_ASSERT1(pt.size() == 3, "Dimensions mismatch");
455  if (pt[2] < .0) pt[2] = 0.;
456  for (short_type f = 1; f < 5; ++f) {
457  scalar_type reldist = gmm::vect_sp(normals_[f], pt)*sqrt(2.);
458  if (reldist > 1.)
459  gmm::scale(pt, 1./reldist);
460  }
461  } else
462  basic_convex_ref_->project_into(pt);
463  }
464 
465  pyramid_QK_of_ref_(dim_type k) : convex_of_reference(pyramid_QK_structure(k), k == 1) {
466  GMM_ASSERT1(k == 1 || k == 2,
467  "Sorry exist only in degree 1 or 2, not " << k);
468 
469  convex<base_node>::points().resize(cvs->nb_points());
470  normals_.resize(cvs->nb_faces());
471  if (k != 1) basic_convex_ref_ = pyramid_QK_of_reference(1);
472 
473  normals_[0] = { 0., 0., -1.};
474  normals_[1] = { 0.,-1., 1.};
475  normals_[2] = { 1., 0., 1.};
476  normals_[3] = { 0., 1., 1.};
477  normals_[4] = {-1., 0., 1.};
478 
479  for (size_type i = 0; i < normals_.size(); ++i)
480  gmm::scale(normals_[i], 1. / gmm::vect_norm2(normals_[i]));
481 
482  if (k==1) {
483  convex<base_node>::points()[0] = base_node(-1.0, -1.0, 0.0);
484  convex<base_node>::points()[1] = base_node( 1.0, -1.0, 0.0);
485  convex<base_node>::points()[2] = base_node(-1.0, 1.0, 0.0);
486  convex<base_node>::points()[3] = base_node( 1.0, 1.0, 0.0);
487  convex<base_node>::points()[4] = base_node( 0.0, 0.0, 1.0);
488  } else if (k==2) {
489  convex<base_node>::points()[0] = base_node(-1.0, -1.0, 0.0);
490  convex<base_node>::points()[1] = base_node( 0.0, -1.0, 0.0);
491  convex<base_node>::points()[2] = base_node( 1.0, -1.0, 0.0);
492  convex<base_node>::points()[3] = base_node(-1.0, 0.0, 0.0);
493  convex<base_node>::points()[4] = base_node( 0.0, 0.0, 0.0);
494  convex<base_node>::points()[5] = base_node( 1.0, 0.0, 0.0);
495  convex<base_node>::points()[6] = base_node(-1.0, 1.0, 0.0);
496  convex<base_node>::points()[7] = base_node( 0.0, 1.0, 0.0);
497  convex<base_node>::points()[8] = base_node( 1.0, 1.0, 0.0);
498  convex<base_node>::points()[9] = base_node(-0.5, -0.5, 0.5);
499  convex<base_node>::points()[10] = base_node( 0.5, -0.5, 0.5);
500  convex<base_node>::points()[11] = base_node(-0.5, 0.5, 0.5);
501  convex<base_node>::points()[12] = base_node( 0.5, 0.5, 0.5);
502  convex<base_node>::points()[13] = base_node( 0.0, 0.0, 1.0);
503  }
504  ppoints = store_point_tab(convex<base_node>::points());
505  if (auto_basic) simplexify_convex(this, *psimplexified_convex);
506  }
507  };
508 
509 
510  DAL_SIMPLE_KEY(pyramid_QK_reference_key_, dim_type);
511 
512  pconvex_ref pyramid_QK_of_reference(dim_type k) {
513  dal::pstatic_stored_object_key
514  pk = std::make_shared<pyramid_QK_reference_key_>(k);
515  dal::pstatic_stored_object o = dal::search_stored_object(pk);
516  if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
517  pconvex_ref p = std::make_shared<pyramid_QK_of_ref_>(k);
518  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
519  dal::PERMANENT_STATIC_OBJECT);
520  pconvex_ref p1 = basic_convex_ref(p);
521  if (p != p1) add_dependency(p, p1);
522  return p;
523  }
524 
525 
526  /* ******************************************************************** */
527  /* Incomplete quadratic pyramidal element of reference. */
528  /* ******************************************************************** */
529 
530  class pyramid_Q2_incomplete_of_ref_ : public convex_of_reference {
531  public :
532  scalar_type is_in(const base_node& pt) const
533  { return basic_convex_ref_->is_in(pt); }
534  scalar_type is_in_face(short_type f, const base_node& pt) const
535  { return basic_convex_ref_->is_in_face(f, pt); }
536 
537  pyramid_Q2_incomplete_of_ref_() : convex_of_reference(pyramid_Q2_incomplete_structure(), false) {
538  convex<base_node>::points().resize(cvs->nb_points());
539  normals_.resize(cvs->nb_faces());
540  basic_convex_ref_ = pyramid_QK_of_reference(1);
541 
542  normals_ = basic_convex_ref_->normals();
543 
544  convex<base_node>::points()[0] = base_node(-1.0, -1.0, 0.0);
545  convex<base_node>::points()[1] = base_node( 0.0, -1.0, 0.0);
546  convex<base_node>::points()[2] = base_node( 1.0, -1.0, 0.0);
547  convex<base_node>::points()[3] = base_node(-1.0, 0.0, 0.0);
548  convex<base_node>::points()[4] = base_node( 1.0, 0.0, 0.0);
549  convex<base_node>::points()[5] = base_node(-1.0, 1.0, 0.0);
550  convex<base_node>::points()[6] = base_node( 0.0, 1.0, 0.0);
551  convex<base_node>::points()[7] = base_node( 1.0, 1.0, 0.0);
552  convex<base_node>::points()[8] = base_node(-0.5, -0.5, 0.5);
553  convex<base_node>::points()[9] = base_node( 0.5, -0.5, 0.5);
554  convex<base_node>::points()[10] = base_node(-0.5, 0.5, 0.5);
555  convex<base_node>::points()[11] = base_node( 0.5, 0.5, 0.5);
556  convex<base_node>::points()[12] = base_node( 0.0, 0.0, 1.0);
557 
558  ppoints = store_point_tab(convex<base_node>::points());
559  if (auto_basic) simplexify_convex(this, *psimplexified_convex);
560  }
561  };
562 
563 
564  DAL_SIMPLE_KEY(pyramid_Q2_incomplete_reference_key_, dim_type);
565 
566  pconvex_ref pyramid_Q2_incomplete_of_reference() {
567  dal::pstatic_stored_object_key
568  pk = std::make_shared<pyramid_Q2_incomplete_reference_key_>(0);
569  dal::pstatic_stored_object o = dal::search_stored_object(pk);
570  if (o)
571  return std::dynamic_pointer_cast<const convex_of_reference>(o);
572  else {
573  pconvex_ref p = std::make_shared<pyramid_Q2_incomplete_of_ref_>();
574  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
575  dal::PERMANENT_STATIC_OBJECT);
576  pconvex_ref p1 = basic_convex_ref(p);
577  if (p != p1) add_dependency(p, p1);
578  return p;
579  }
580  }
581 
582 
583  /* ******************************************************************** */
584  /* Incomplete quadratic triangular prism element of reference. */
585  /* ******************************************************************** */
586 
587  class prism_incomplete_P2_of_ref_ : public convex_of_reference {
588  public :
589  scalar_type is_in(const base_node& pt) const
590  { return basic_convex_ref_->is_in(pt); }
591  scalar_type is_in_face(short_type f, const base_node& pt) const
592  { return basic_convex_ref_->is_in_face(f, pt); }
593 
594  prism_incomplete_P2_of_ref_() : convex_of_reference(prism_incomplete_P2_structure(), false) {
595  convex<base_node>::points().resize(cvs->nb_points());
596  normals_.resize(cvs->nb_faces());
597  basic_convex_ref_ = prism_of_reference(3);
598 
599  normals_ = basic_convex_ref_->normals();
600 
601  convex<base_node>::points()[0] = base_node(0.0, 0.0, 0.0);
602  convex<base_node>::points()[1] = base_node(0.5, 0.0, 0.0);
603  convex<base_node>::points()[2] = base_node(1.0, 0.0, 0.0);
604  convex<base_node>::points()[3] = base_node(0.0, 0.5, 0.0);
605  convex<base_node>::points()[4] = base_node(0.5, 0.5, 0.0);
606  convex<base_node>::points()[5] = base_node(0.0, 1.0, 0.0);
607  convex<base_node>::points()[6] = base_node(0.0, 0.0, 0.5);
608  convex<base_node>::points()[7] = base_node(1.0, 0.0, 0.5);
609  convex<base_node>::points()[8] = base_node(0.0, 1.0, 0.5);
610  convex<base_node>::points()[9] = base_node(0.0, 0.0, 1.0);
611  convex<base_node>::points()[10] = base_node(0.5, 0.0, 1.0);
612  convex<base_node>::points()[11] = base_node(1.0, 0.0, 1.0);
613  convex<base_node>::points()[12] = base_node(0.0, 0.5, 1.0);
614  convex<base_node>::points()[13] = base_node(0.5, 0.5, 1.0);
615  convex<base_node>::points()[14] = base_node(0.0, 1.0, 1.0);
616 
617  ppoints = store_point_tab(convex<base_node>::points());
618  if (auto_basic) simplexify_convex(this, *psimplexified_convex);
619  }
620  };
621 
622 
623  DAL_SIMPLE_KEY(prism_incomplete_P2_reference_key_, dim_type);
624 
625  pconvex_ref prism_incomplete_P2_of_reference() {
626  dal::pstatic_stored_object_key
627  pk = std::make_shared<prism_incomplete_P2_reference_key_>(0);
628  dal::pstatic_stored_object o = dal::search_stored_object(pk);
629  if (o)
630  return std::dynamic_pointer_cast<const convex_of_reference>(o);
631  else {
632  pconvex_ref p = std::make_shared<prism_incomplete_P2_of_ref_>();
633  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
634  dal::PERMANENT_STATIC_OBJECT);
635  pconvex_ref p1 = basic_convex_ref(p);
636  if (p != p1) add_dependency(p, p1);
637  return p;
638  }
639  }
640 
641 
642  /* ******************************************************************** */
643  /* Products. */
644  /* ******************************************************************** */
645 
646  DAL_DOUBLE_KEY(product_ref_key_, pconvex_ref, pconvex_ref);
647 
648  struct product_ref_ : public convex_of_reference {
649  pconvex_ref cvr1, cvr2;
650 
651  scalar_type is_in(const base_node &pt) const {
652  dim_type n1 = cvr1->structure()->dim(), n2 = cvr2->structure()->dim();
653  base_node pt1(n1), pt2(n2);
654  GMM_ASSERT1(pt.size() == cvs->dim(),
655  "product_ref_::is_in: Dimension does not match");
656  std::copy(pt.begin(), pt.begin()+n1, pt1.begin());
657  std::copy(pt.begin()+n1, pt.end(), pt2.begin());
658  return std::max(cvr1->is_in(pt1), cvr2->is_in(pt2));
659  }
660 
661  scalar_type is_in_face(short_type f, const base_node &pt) const {
662  // ne controle pas si le point est dans le convexe mais si un point
663  // suppos\E9 appartenir au convexe est dans une face donn\E9e
664  dim_type n1 = cvr1->structure()->dim(), n2 = cvr2->structure()->dim();
665  base_node pt1(n1), pt2(n2);
666  GMM_ASSERT1(pt.size() == cvs->dim(), "Dimensions mismatch");
667  std::copy(pt.begin(), pt.begin()+n1, pt1.begin());
668  std::copy(pt.begin()+n1, pt.end(), pt2.begin());
669  if (f < cvr1->structure()->nb_faces()) return cvr1->is_in_face(f, pt1);
670  else return cvr2->is_in_face(short_type(f - cvr1->structure()->nb_faces()), pt2);
671  }
672 
673  void project_into(base_node &pt) const {
674  if (auto_basic) {
675  GMM_ASSERT1(pt.size() == cvs->dim(), "Dimensions mismatch");
676  dim_type n1 = cvr1->structure()->dim(), n2 = cvr2->structure()->dim();
677  base_node pt1(n1), pt2(n2);
678  std::copy(pt.begin(), pt.begin()+n1, pt1.begin());
679  std::copy(pt.begin()+n1, pt.end(), pt2.begin());
680  cvr1->project_into(pt1);
681  cvr2->project_into(pt2);
682  std::copy(pt1.begin(), pt1.end(), pt.begin());
683  std::copy(pt2.begin(), pt2.end(), pt.begin()+n1);
684  } else
685  basic_convex_ref_->project_into(pt);
686  }
687 
688  product_ref_(pconvex_ref a, pconvex_ref b) :
689  convex_of_reference(
690  convex_direct_product(*a, *b).structure(),
691  basic_convex_ref(a) == a && basic_convex_ref(b) == b)
692  {
693  if (a->structure()->dim() < b->structure()->dim())
694  GMM_WARNING1("Illegal convex: swap your operands: dim(cv1)=" <<
695  int(a->structure()->dim()) << " < dim(cv2)=" <<
696  int(b->structure()->dim()));
697  cvr1 = a; cvr2 = b;
698  *((convex<base_node> *)(this)) = convex_direct_product(*a, *b);
699  normals_.resize(cvs->nb_faces());
700  base_small_vector null(cvs->dim()); null.fill(0.0);
701  std::fill(normals_.begin(), normals_.end(), null);
702  for (size_type r = 0; r < cvr1->structure()->nb_faces(); r++)
703  std::copy(cvr1->normals()[r].begin(), cvr1->normals()[r].end(),
704  normals_[r].begin());
705  for (size_type r = 0; r < cvr2->structure()->nb_faces(); r++)
706  std::copy(cvr2->normals()[r].begin(), cvr2->normals()[r].end(),
707  normals_[r+cvr1->structure()->nb_faces()].begin()
708  + cvr1->structure()->dim());
709  ppoints = store_point_tab(convex<base_node>::points());
710 
711  if (basic_convex_ref(a) != a || basic_convex_ref(b) != b)
712  basic_convex_ref_ = convex_ref_product(basic_convex_ref(a),
713  basic_convex_ref(b));
714  if (auto_basic) simplexify_convex(this, *psimplexified_convex);
715  }
716  };
717 
718 
719  pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b) {
720  dal::pstatic_stored_object_key
721  pk = std::make_shared<product_ref_key_>(a, b);
722  dal::pstatic_stored_object o = dal::search_stored_object(pk);
723  if (o)
724  return std::dynamic_pointer_cast<const convex_of_reference>(o);
725  else {
726  pconvex_ref p = std::make_shared<product_ref_>(a, b);
727  dal::add_stored_object(pk, p, a, b,
728  convex_product_structure(a->structure(),
729  b->structure()),
730  p->pspt(), dal::PERMANENT_STATIC_OBJECT);
731  pconvex_ref p1 = basic_convex_ref(p);
732  if (p != p1) add_dependency(p, p1);
733  return p;
734  }
735  }
736 
737  pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k) {
738  if (nc <= 1) return simplex_of_reference(nc,k);
739  return convex_ref_product(parallelepiped_of_reference(dim_type(nc-1),k),
740  simplex_of_reference(k));
741  }
742 
743  pconvex_ref prism_of_reference(dim_type nc) {
744  if (nc <= 2) return parallelepiped_of_reference(nc);
745  else return convex_ref_product(simplex_of_reference(dim_type(nc-1)),
746  simplex_of_reference(1));
747  }
748 
749  /* equilateral ref convexes are used for estimation of convex quality */
750  class equilateral_simplex_of_ref_ : public convex_of_reference {
751  scalar_type r_inscr;
752  public:
753  scalar_type is_in(const base_node &pt) const {
754  GMM_ASSERT1(pt.size() == cvs->dim(), "Dimension does not match");
755  scalar_type d(0);
756  for (size_type f = 0; f < normals().size(); ++f) {
757  const base_node &x0 = (f ? convex<base_node>::points()[f-1]
758  : convex<base_node>::points().back());
759  scalar_type v = gmm::vect_sp(pt-x0, normals()[f]);
760  if (f == 0) d = v; else d = std::max(d,v);
761  }
762  return d;
763  }
764 
765  scalar_type is_in_face(short_type f, const base_node &pt) const {
766  GMM_ASSERT1(pt.size() == cvs->dim(), "Dimension does not match");
767  const base_node &x0 = (f ? convex<base_node>::points()[f-1]
768  : convex<base_node>::points().back());
769  return gmm::vect_sp(pt-x0, normals()[f]);
770  }
771 
772  void project_into(base_node &pt) const {
773  dim_type N = cvs->dim();
774  GMM_ASSERT1(pt.size() == N, "Dimension does not match");
775  base_node G(N); G.fill(0.);
776  for (const base_node &x : convex<base_node>::points())
777  G += x;
778  gmm::scale(G, scalar_type(1)/scalar_type(N+1));
779  for (size_type f = 0; f < normals().size(); ++f) {
780  scalar_type r = gmm::vect_sp(pt-G, normals()[f]);
781  if (r > r_inscr)
782  pt = G + r_inscr/r*(pt-G);
783  }
784 
785  }
786 
787  equilateral_simplex_of_ref_(size_type N) :
788  convex_of_reference(simplex_structure(dim_type(N), 1), true)
789  {
790  //https://math.stackexchange.com/questions/2739915/radius-of-inscribed-sphere-of-n-simplex
791  r_inscr = scalar_type(1)/sqrt(scalar_type(2*N)*scalar_type(N+1));
792 
793  pconvex_ref prev = equilateral_simplex_of_reference(dim_type(N-1));
794  convex<base_node>::points().resize(N+1);
795  normals_.resize(N+1);
796  base_node G(N); G.fill(0.);
797  for (short_type i=0; i < N+1; ++i) {
798  convex<base_node>::points()[i].resize(N);
799  if (i != N) {
800  std::copy(prev->convex<base_node>::points()[i].begin(),
801  prev->convex<base_node>::points()[i].end(),
802  convex<base_node>::points()[i].begin());
803  convex<base_node>::points()[i][N-1] = 0.;
804  } else {
805  convex<base_node>::points()[i] = scalar_type(1)/scalar_type(N) * G;
806  convex<base_node>::points()[i][N-1]
807  = sqrt(1. - gmm::vect_norm2_sqr(convex<base_node>::points()[i]));
808  }
809  G += convex<base_node>::points()[i];
810  }
811  gmm::scale(G, scalar_type(1)/scalar_type(N+1));
812  for (size_type f=0; f < N+1; ++f) {
813  normals_[f] = G - convex<base_node>::points()[f];
814  gmm::scale(normals_[f], 1/gmm::vect_norm2(normals_[f]));
815  }
816  ppoints = store_point_tab(convex<base_node>::points());
817  if (auto_basic) simplexify_convex(this, *psimplexified_convex);
818  }
819  };
820 
821  pconvex_ref equilateral_simplex_of_reference(dim_type nc) {
822  if (nc <= 1) return simplex_of_reference(nc);
823  dal::pstatic_stored_object_key
824  pk = std::make_shared<convex_of_reference_key>(1, nc);
825  dal::pstatic_stored_object o = dal::search_stored_object(pk);
826  if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
827  pconvex_ref p = std::make_shared<equilateral_simplex_of_ref_>(nc);
828  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
829  dal::PERMANENT_STATIC_OBJECT);
830  return p;
831  }
832 
833 
834  /* generic convex with n global nodes */
835 
836  class generic_dummy_ : public convex_of_reference {
837  public:
838  scalar_type is_in(const base_node &) const
839  { GMM_ASSERT1(false, "Information not available here"); }
840  scalar_type is_in_face(short_type, const base_node &) const
841  { GMM_ASSERT1(false, "Information not available here"); }
842  void project_into(base_node &) const
843  { GMM_ASSERT1(false, "Operation not available here"); }
844 
845  generic_dummy_(dim_type d, size_type n, short_type nf) :
846  convex_of_reference(generic_dummy_structure(d, n, nf), true)
847  {
848  convex<base_node>::points().resize(n);
849  normals_.resize(0);
850  base_node P(d);
851  std::fill(P.begin(), P.end(), scalar_type(1)/scalar_type(20));
852  std::fill(convex<base_node>::points().begin(), convex<base_node>::points().end(), P);
853  ppoints = store_point_tab(convex<base_node>::points());
854  }
855  };
856 
857  pconvex_ref generic_dummy_convex_ref(dim_type nc, size_type n,
858  short_type nf) {
859  dal::pstatic_stored_object_key
860  pk = std::make_shared<convex_of_reference_key>(2, nc, short_type(n), nf);
861  dal::pstatic_stored_object o = dal::search_stored_object(pk);
862  if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
863  pconvex_ref p = std::make_shared<generic_dummy_>(nc, n, nf);
864  dal::add_stored_object(pk, p, p->structure(), p->pspt(),
865  dal::PERMANENT_STATIC_OBJECT);
866  return p;
867  }
868 
869 
870 } /* end of namespace bgeot. */
bgeot_comma_init.h
convenient initialization of vectors via overload of "operator,".
bgeot_convex_ref.h
Reference convexes.
bgeot_mesh_structure.h
Mesh structure definition.
bgeot::convex_of_reference
Base class for reference convexes.
Definition: bgeot_convex_ref.h:91
bgeot::convex_of_reference::points
const stored_point_tab & points() const
return the vertices of the reference convex.
Definition: bgeot_convex_ref.h:122
bgeot::convex_of_reference::simplexified_convex
const mesh_structure * simplexified_convex() const
return a mesh structure composed of simplexes whose union is the reference convex.
Definition: bgeot_convex_ref.cc:207
bgeot::convex_of_reference::normals
const std::vector< base_small_vector > & normals() const
return the normal vector for each face.
Definition: bgeot_convex_ref.h:119
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:72
bgeot::mesh_structure::clear
void clear()
erase the mesh
Definition: bgeot_mesh_structure.cc:204
dal_singleton.h
A simple singleton implementation.
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:26
bgeot::prism_incomplete_P2_structure
pconvex_structure prism_incomplete_P2_structure()
Give a pointer on the 3D quadratic incomplete prism structure.
Definition: bgeot_convex_structure.cc:670
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:719
bgeot::pyramid_Q2_incomplete_structure
pconvex_structure pyramid_Q2_incomplete_structure()
Give a pointer on the 3D quadratic incomplete pyramid structure.
Definition: bgeot_convex_structure.cc:610
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_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:513
bgeot::equilateral_simplex_of_reference
pconvex_ref equilateral_simplex_of_reference(dim_type nc)
equilateral simplex (degree 1).
Definition: bgeot_convex_ref.cc:821
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::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:330
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:53
bgeot::generic_dummy_structure
pconvex_structure generic_dummy_structure(dim_type nc, size_type n, short_type nf)
Generic convex with n global nodes.
Definition: bgeot_convex_structure.cc:730
bgeot::pyramid_Q2_incomplete_of_reference
pconvex_ref pyramid_Q2_incomplete_of_reference()
incomplete quadratic pyramidal element of reference (13-node)
Definition: bgeot_convex_ref.cc:566
bgeot::prism_incomplete_P2_of_reference
pconvex_ref prism_incomplete_P2_of_reference()
incomplete quadratic prism element of reference (15-node)
Definition: bgeot_convex_ref.cc:625
bgeot::Q2_incomplete_structure
pconvex_structure Q2_incomplete_structure(dim_type nc)
Give a pointer on the structures of a incomplete Q2 quadrilateral/hexahedral of dimension d = 2 or 3.
Definition: bgeot_convex_structure.cc:433
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::generic_dummy_convex_ref
pconvex_ref generic_dummy_convex_ref(dim_type nc, size_type n, short_type nf)
generic convex with n global nodes
Definition: bgeot_convex_ref.cc:857
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:147
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::Q2_incomplete_of_reference
pconvex_ref Q2_incomplete_of_reference(dim_type nc)
incomplete Q2 quadrilateral/hexahedral of reference of dimension d = 2 or 3
Definition: bgeot_convex_ref.cc:416
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:373
bgeot::basic_structure
pconvex_structure basic_structure(pconvex_structure cv)
Original structure (if concerned)
Definition: bgeot_convex_structure.h:171
bgeot::basic_convex_ref
pconvex_ref basic_convex_ref(pconvex_ref cvr)
return the associated order 1 reference convex.
Definition: bgeot_convex_ref.h:137
bgeot::prism_of_reference
pconvex_ref prism_of_reference(dim_type nc)
prism of reference of dimension nc (and degree 1)
Definition: bgeot_convex_ref.cc:743
dal::add_stored_object
void add_stored_object(pstatic_stored_object_key k, pstatic_stored_object o, permanence perm)
Add an object with two optional dependencies.
Definition: dal_static_stored_objects.cc:283
dal::search_stored_object
pstatic_stored_object search_stored_object(pstatic_stored_object_key k)
Gives a pointer to an object from a key pointer.
Definition: dal_static_stored_objects.cc:176
bgeot::stored_point_tab
Point tab storage.
Definition: bgeot_convex_ref.h:49