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