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