gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
automaticMeshSizeField.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributed by Arthur Bawin
7 
8 #include "automaticMeshSizeField.h"
9 #include "GModel.h"
10 #include "GRegion.h"
11 #include "GFace.h"
12 #include "GEntity.h"
13 #include "MVertex.h"
14 #include "MTetrahedron.h"
15 #include "SBoundingBox3d.h"
16 #include "GmshMessage.h"
17 #include "curvature.h"
18 #include "Numeric.h"
19 #include "robustPredicates.h"
20 
21 #include <queue>
22 #include <inttypes.h>
23 
24 #ifdef HAVE_HXT
25 extern "C" {
26 #include "hxt_tools.h"
27 #include "hxt_edge.h"
28 #include "hxt_bbox.h"
29 #include "hxt_tetMesh.h"
30 #include "hxt_tetUtils.h"
31 #include "hxt_tetFlag.h"
32 #include "hxt_tetDelaunay.h"
33 #include "hxt_tetRefine.h"
34 }
35 #endif
36 
37 #if defined(HAVE_P4EST)
38 #include <p8est_search.h>
39 #endif
40 
41 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
42 static inline void norme2(double v[3], double *norme2)
43 {
44  *norme2 = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
45 }
46 
47 static inline bool isPoint(double x1, double y1, double z1, double x2,
48  double y2, double z2, double tol)
49 {
50  return (fabs(x2 - x1) < tol && fabs(y2 - y1) < tol && fabs(z2 - z1) < tol);
51 }
52 
53 void writeNodalCurvature(double *nodalCurvature, int size, const char *filename)
54 {
55  FILE *f = fopen(filename, "w");
56  if(f == nullptr) {
57  printf("Erreur : fileOutput == NULL\n");
58  exit(-1);
59  }
60 
61  for(int i = 0; i < size; ++i) {
62  fprintf(f, "%f %f %f - %d\n", nodalCurvature[6 * i + 0],
63  nodalCurvature[6 * i + 1], nodalCurvature[6 * i + 2], i);
64  fprintf(f, "%f %f %f\n", nodalCurvature[6 * i + 3],
65  nodalCurvature[6 * i + 4], nodalCurvature[6 * i + 5]);
66  }
67  fclose(f);
68 }
69 
70 static HXTStatus getAllFacesOfAllRegions(std::vector<GRegion *> &regions,
71  HXTMesh *m,
72  std::vector<GFace *> &allFaces)
73 {
74  std::set<GFace *, GEntityPtrLessThan> allFacesSet;
75  if(m) {
76  m->brep.numVolumes = regions.size();
77  HXT_CHECK(hxtAlignedMalloc(&m->brep.numSurfacesPerVolume,
78  m->brep.numVolumes * sizeof(uint32_t)));
79  }
80  uint32_t to_alloc = 0;
81  for(std::size_t i = 0; i < regions.size(); i++) {
82  std::vector<GFace *> const &f = regions[i]->faces();
83  std::vector<GFace *> const &f_e = regions[i]->embeddedFaces();
84  if(m) {
85  m->brep.numSurfacesPerVolume[i] = f.size() + f_e.size();
86  to_alloc += m->brep.numSurfacesPerVolume[i];
87  }
88  allFacesSet.insert(f.begin(), f.end());
89  allFacesSet.insert(f_e.begin(), f_e.end());
90  }
91  allFaces.insert(allFaces.begin(), allFacesSet.begin(), allFacesSet.end());
92 
93  if(!m) return HXT_STATUS_OK;
94 
95  HXT_CHECK(
96  hxtAlignedMalloc(&m->brep.surfacesPerVolume, to_alloc * sizeof(uint32_t)));
97 
98  uint32_t counter = 0;
99  for(std::size_t i = 0; i < regions.size(); i++) {
100  std::vector<GFace *> const &f = regions[i]->faces();
101  std::vector<GFace *> const &f_e = regions[i]->embeddedFaces();
102  for(size_t j = 0; j < f.size(); j++)
103  m->brep.surfacesPerVolume[counter++] = f[j]->tag();
104  for(size_t j = 0; j < f_e.size(); j++)
105  m->brep.surfacesPerVolume[counter++] = f_e[j]->tag();
106  }
107 
108  // printf("volume 0 has %d faces\n",m->brep.numSurfacesPerVolume[0]);
109  // for (int
110  // i=0;i<m->brep.numSurfacesPerVolume[0];i++)printf("%d",m->brep.surfacesPerVolume[i]);
111  // printf("\n");
112 
113  return HXT_STATUS_OK;
114 }
115 
116 static HXTStatus getAllEdgesOfAllFaces(std::vector<GFace *> &faces, HXTMesh *m,
117  std::vector<GEdge *> &allEdges)
118 {
119  if(m) {
120  m->brep.numSurfaces = faces.size();
121  HXT_CHECK(hxtAlignedMalloc(&m->brep.numCurvesPerSurface,
122  m->brep.numSurfaces * sizeof(uint32_t)));
123  }
124  uint32_t to_alloc = 0;
125 
126  std::set<GEdge *, GEntityPtrLessThan> allEdgesSet;
127  for(std::size_t i = 0; i < faces.size(); i++) {
128  std::vector<GEdge *> const &f = faces[i]->edges();
129  std::vector<GEdge *> const &f_e = faces[i]->embeddedEdges();
130  if(m) {
131  m->brep.numCurvesPerSurface[i] = f.size() + f_e.size();
132  to_alloc += m->brep.numCurvesPerSurface[i];
133  }
134  allEdgesSet.insert(f.begin(), f.end());
135  allEdgesSet.insert(f_e.begin(), f_e.end());
136  }
137  allEdges.insert(allEdges.begin(), allEdgesSet.begin(), allEdgesSet.end());
138 
139  if(!m) return HXT_STATUS_OK;
140 
141  HXT_CHECK(
142  hxtAlignedMalloc(&m->brep.curvesPerSurface, to_alloc * sizeof(uint32_t)));
143 
144  uint32_t counter = 0;
145  for(std::size_t i = 0; i < faces.size(); i++) {
146  std::vector<GEdge *> const &f = faces[i]->edges();
147  std::vector<GEdge *> const &f_e = faces[i]->embeddedEdges();
148  for(size_t j = 0; j < f.size(); j++)
149  m->brep.curvesPerSurface[counter++] = f[j]->tag();
150  for(size_t j = 0; j < f_e.size(); j++)
151  m->brep.curvesPerSurface[counter++] = f_e[j]->tag();
152  }
153  return HXT_STATUS_OK;
154 }
155 
156 HXTStatus Gmsh2Hxt(std::vector<GRegion *> &regions, HXTMesh *m,
157  std::map<MVertex *, uint32_t> &v2c,
158  std::vector<MVertex *> &c2v);
159 
160 HXTStatus Gmsh2Hxt(std::vector<GFace *> &faces, HXTMesh *m,
161  std::map<MVertex *, uint32_t> &v2c,
162  std::vector<MVertex *> &c2v)
163 {
164  std::vector<GEdge *> edges;
165  HXT_CHECK(getAllEdgesOfAllFaces(faces, m, edges));
166  std::set<MVertex *> all;
167 
168  uint64_t ntri = 0;
169  uint64_t nedg = 0;
170 
171  for(size_t j = 0; j < edges.size(); j++) {
172  GEdge *ge = edges[j];
173  nedg += ge->lines.size();
174  for(size_t i = 0; i < ge->lines.size(); i++) {
175  all.insert(ge->lines[i]->getVertex(0));
176  all.insert(ge->lines[i]->getVertex(1));
177  }
178  }
179 
180  for(size_t j = 0; j < faces.size(); j++) {
181  GFace *gf = faces[j];
182  ntri += gf->triangles.size();
183  for(size_t i = 0; i < gf->triangles.size(); i++) {
184  all.insert(gf->triangles[i]->getVertex(0));
185  all.insert(gf->triangles[i]->getVertex(1));
186  all.insert(gf->triangles[i]->getVertex(2));
187  }
188  }
189 
190  m->vertices.num = m->vertices.size = all.size();
191  HXT_CHECK(
192  hxtAlignedMalloc(&m->vertices.coord, 4 * m->vertices.num * sizeof(double)));
193 
194  size_t count = 0;
195  c2v.resize(all.size());
196  // for(auto it = all.begin(); it != all.end(); it++)
197  // {
198  // m->vertices.coord[4 * count + 0] = (*it)->x();
199  // m->vertices.coord[4 * count + 1] = (*it)->y();
200  // m->vertices.coord[4 * count + 2] = (*it)->z();
201  // m->vertices.coord[4 * count + 3] = 0.0;
202  // v2c[*it] = count;
203  // c2v[count++] = *it;
204  // }
205  for(MVertex *v : all) {
206  m->vertices.coord[4 * count + 0] = v->x();
207  m->vertices.coord[4 * count + 1] = v->y();
208  m->vertices.coord[4 * count + 2] = v->z();
209  m->vertices.coord[4 * count + 3] = 0;
210  // if(CTX::instance()->mesh.lcFromPoints) { // size on embedded points in
211  // volume
212  // auto it = vlc.find(v);
213  // if(it != vlc.end())
214  // m->vertices.coord[4 * count + 3] = it->second;
215  // }
216  v2c[v] = count;
217  c2v[count++] = v;
218  }
219  all.clear();
220 
221  m->lines.num = m->lines.size = nedg;
222  uint64_t index = 0;
223 
224  HXT_CHECK(
225  hxtAlignedMalloc(&m->lines.node, (m->lines.num) * 2 * sizeof(uint32_t)));
226  HXT_CHECK(
227  hxtAlignedMalloc(&m->lines.color, (m->lines.num) * sizeof(uint32_t)));
228 
229  for(size_t j = 0; j < edges.size(); j++) {
230  GEdge *ge = edges[j];
231  for(size_t i = 0; i < ge->lines.size(); i++) {
232  m->lines.node[2 * index + 0] = v2c[ge->lines[i]->getVertex(0)];
233  m->lines.node[2 * index + 1] = v2c[ge->lines[i]->getVertex(1)];
234  m->lines.color[index] = ge->tag();
235  index++;
236  }
237  }
238 
239  m->triangles.num = m->triangles.size = ntri;
240  HXT_CHECK(hxtAlignedMalloc(&m->triangles.node,
241  (m->triangles.num) * 3 * sizeof(uint32_t)));
242  HXT_CHECK(hxtAlignedMalloc(&m->triangles.color,
243  (m->triangles.num) * sizeof(uint32_t)));
244 
245  index = 0;
246  for(size_t j = 0; j < faces.size(); j++) {
247  GFace *gf = faces[j];
248  for(size_t i = 0; i < gf->triangles.size(); i++) {
249  m->triangles.node[3 * index + 0] = v2c[gf->triangles[i]->getVertex(0)];
250  m->triangles.node[3 * index + 1] = v2c[gf->triangles[i]->getVertex(1)];
251  m->triangles.node[3 * index + 2] = v2c[gf->triangles[i]->getVertex(2)];
252  m->triangles.color[index] = gf->tag();
253  index++;
254  }
255  }
256  return HXT_STATUS_OK;
257 }
258 
259 // HXTStatus Gmsh2Hxt(std::vector<GRegion *> &regions, HXTMesh *m,
260 // std::map<MVertex *, uint32_t> &v2c, std::vector<MVertex *> &c2v);
261 
262 // HXTStatus Gmsh2Hxt(std::vector<GFace *> &faces, HXTMesh *m, std::map<MVertex
263 // *, uint32_t> &v2c, std::vector<MVertex *> &c2v);
264 
265 HXTStatus Gmsh2HxtLocal(std::vector<GFace *> &faces, HXTMesh *m,
266  std::map<MVertex *, uint32_t> &v2c,
267  std::vector<MVertex *> &c2v)
268 {
269  std::vector<GEdge *> edges;
270  HXT_CHECK(getAllEdgesOfAllFaces(faces, m, edges));
271  std::set<MVertex *> all;
272 
273  uint64_t ntri = 0;
274  uint64_t nedg = 0;
275 
276  for(size_t j = 0; j < edges.size(); j++) {
277  GEdge *ge = edges[j];
278  nedg += ge->lines.size();
279  for(size_t i = 0; i < ge->lines.size(); i++) {
280  // all.insert(ge->lines[i]->getVertex(0));
281  // all.insert(ge->lines[i]->getVertex(1));
282  }
283  }
284 
285  for(size_t j = 0; j < faces.size(); j++) {
286  GFace *gf = faces[j];
287  ntri += gf->triangles.size();
288  for(size_t i = 0; i < gf->triangles.size(); i++) {
289  all.insert(gf->triangles[i]->getVertex(0));
290  all.insert(gf->triangles[i]->getVertex(1));
291  all.insert(gf->triangles[i]->getVertex(2));
292  }
293  }
294 
295  // printf("%d vertices %d triangles\n",all.size(),ntri);
296 
297  m->vertices.num = m->vertices.size = all.size();
298  HXT_CHECK(
299  hxtAlignedMalloc(&m->vertices.coord, 4 * m->vertices.num * sizeof(double)));
300 
301  size_t count = 0;
302  c2v.resize(all.size());
303  for(auto it = all.begin(); it != all.end(); it++) {
304  m->vertices.coord[4 * count + 0] = (*it)->x();
305  m->vertices.coord[4 * count + 1] = (*it)->y();
306  m->vertices.coord[4 * count + 2] = (*it)->z();
307  m->vertices.coord[4 * count + 3] = 0.0;
308  v2c[*it] = count;
309  c2v[count++] = *it;
310  }
311  all.clear();
312 
313  m->lines.num = m->lines.size = nedg;
314  uint64_t index = 0;
315 
316  HXT_CHECK(
317  hxtAlignedMalloc(&m->lines.node, (m->lines.num) * 2 * sizeof(uint32_t)));
318  HXT_CHECK(
319  hxtAlignedMalloc(&m->lines.color, (m->lines.num) * sizeof(uint32_t)));
320 
321  for(size_t j = 0; j < edges.size(); j++) {
322  GEdge *ge = edges[j];
323  for(size_t i = 0; i < ge->lines.size(); i++) {
324  m->lines.node[2 * index + 0] = v2c[ge->lines[i]->getVertex(0)];
325  m->lines.node[2 * index + 1] = v2c[ge->lines[i]->getVertex(1)];
326  m->lines.color[index] = ge->tag();
327  index++;
328  }
329  }
330 
331  m->triangles.num = m->triangles.size = ntri;
332  HXT_CHECK(hxtAlignedMalloc(&m->triangles.node,
333  (m->triangles.num) * 3 * sizeof(uint32_t)));
334  HXT_CHECK(hxtAlignedMalloc(&m->triangles.color,
335  (m->triangles.num) * sizeof(uint32_t)));
336 
337  index = 0;
338  for(size_t j = 0; j < faces.size(); j++) {
339  GFace *gf = faces[j];
340  for(size_t i = 0; i < gf->triangles.size(); i++) {
341  m->triangles.node[3 * index + 0] = v2c[gf->triangles[i]->getVertex(0)];
342  m->triangles.node[3 * index + 1] = v2c[gf->triangles[i]->getVertex(1)];
343  m->triangles.node[3 * index + 2] = v2c[gf->triangles[i]->getVertex(2)];
344  m->triangles.color[index] = gf->tag();
345  index++;
346  }
347  }
348  return HXT_STATUS_OK;
349 }
350 
351 HXTStatus Gmsh2HxtGlobal(std::vector<GFace *> &faces, HXTMesh *m,
352  std::map<MVertex *, uint32_t> &v2c,
353  std::vector<MVertex *> &c2v)
354 {
355  std::vector<GEdge *> edges;
356  HXT_CHECK(getAllEdgesOfAllFaces(faces, m, edges));
357  std::set<MVertex *> all;
358 
359  uint64_t cumsum = 0;
360  uint64_t cumsumNoEdges = 0;
361  for(size_t j = 0; j < edges.size(); j++) {
362  GEdge *ge = edges[j];
363  cumsum += ge->lines.size();
364  for(size_t i = 0; i < ge->lines.size(); i++) {
365  all.insert(ge->lines[i]->getVertex(0));
366  all.insert(ge->lines[i]->getVertex(1));
367  }
368  }
369 
370  for(size_t j = 0; j < faces.size(); j++) {
371  GFace *gf = faces[j];
372  cumsum += gf->triangles.size();
373  cumsumNoEdges += gf->triangles.size();
374  for(size_t i = 0; i < gf->triangles.size(); i++) {
375  all.insert(gf->triangles[i]->getVertex(0));
376  all.insert(gf->triangles[i]->getVertex(1));
377  all.insert(gf->triangles[i]->getVertex(2));
378  }
379  }
380 
381  c2v.resize(cumsum);
382  all.clear();
383 
384  cumsum = 0;
385  size_t count_c2v2 = 0;
386  for(size_t j = 0; j < faces.size(); ++j) {
387  // all propre à la face
388  GFace *gf = faces[j];
389  cumsum += gf->triangles.size();
390  for(size_t i = 0; i < gf->triangles.size(); i++) {
391  all.insert(gf->triangles[i]->getVertex(0));
392  all.insert(gf->triangles[i]->getVertex(1));
393  all.insert(gf->triangles[i]->getVertex(2));
394  }
395  size_t count = 0;
396  for(auto it = all.begin(); it != all.end(); it++) {
397  v2c[*it] = count++;
398  c2v[count_c2v2++] = *it;
399  }
400  all.clear();
401  }
402 
403  return HXT_STATUS_OK;
404 }
405 
406 /* ======================================================================================
407  Functions from hxt_octree
408  ======================================================================================
409  */
410 
411 static bool rtreeCallback(uint64_t id, void *ctx)
412 {
413  std::vector<uint64_t> *vec = reinterpret_cast<std::vector<uint64_t> *>(ctx);
414  vec->push_back(id);
415  return true;
416 }
417 
418 /* Create a p4est connectivity structure for the given bounding box. */
419 static p4est_connectivity_t *
420 p8est_connectivity_new_cube(ForestOptions *forestOptions)
421 {
422  const p4est_topidx_t num_vertices = 8;
423  const p4est_topidx_t num_trees = 1;
424  const p4est_topidx_t num_ett = 0;
425  const p4est_topidx_t num_ctt = 0;
426 
427  double centreX = (forestOptions->bbox[0] + forestOptions->bbox[3]) / 2.0;
428  double centreY = (forestOptions->bbox[1] + forestOptions->bbox[4]) / 2.0;
429  double centreZ = (forestOptions->bbox[2] + forestOptions->bbox[5]) / 2.0;
430  double cX = (forestOptions->bbox[3] - forestOptions->bbox[0]) / 2.0;
431  double cY = (forestOptions->bbox[4] - forestOptions->bbox[1]) / 2.0;
432  double cZ = (forestOptions->bbox[5] - forestOptions->bbox[2]) / 2.0;
433 
434  double scalingFactor =
435  1.3; // The octree is this times bigger than the surface mesh's bounding box
436  double c = scalingFactor * fmax(fmax(cX, cY), cZ);
437 
438  // TODO : Compute any bounding box, not necessarily aligned with the axes
439  const double vertices[8 * 3] = {
440  centreX - c, centreY - c, centreZ - c, centreX + c, centreY - c,
441  centreZ - c, centreX - c, centreY + c, centreZ - c, centreX + c,
442  centreY + c, centreZ - c, centreX - c, centreY - c, centreZ + c,
443  centreX + c, centreY - c, centreZ + c, centreX - c, centreY + c,
444  centreZ + c, centreX + c, centreY + c, centreZ + c,
445  };
446  const p4est_topidx_t tree_to_vertex[1 * 8] = {0, 1, 2, 3, 4, 5, 6, 7};
447  const p4est_topidx_t tree_to_tree[1 * 6] = {0, 0, 0, 0, 0, 0};
448  const int8_t tree_to_face[1 * 6] = {0, 1, 2, 3, 4, 5};
449 
450  return p4est_connectivity_new_copy(num_vertices, num_trees, 0, 0, vertices,
451  tree_to_vertex, tree_to_tree, tree_to_face,
452  nullptr, &num_ett, nullptr, nullptr,
453  nullptr, &num_ctt, nullptr, nullptr);
454 }
455 
456 static p4est_connectivity_t *
457 p8est_connectivity_new_square(ForestOptions *forestOptions)
458 {
459  const p4est_topidx_t num_vertices = 8;
460  const p4est_topidx_t num_trees = 1;
461  const p4est_topidx_t num_ett = 0;
462  const p4est_topidx_t num_ctt = 0;
463 
464  double centreX = (forestOptions->bbox[0] + forestOptions->bbox[3]) / 2.0;
465  double centreY = (forestOptions->bbox[1] + forestOptions->bbox[4]) / 2.0;
466  double cX = (forestOptions->bbox[3] - forestOptions->bbox[0]) / 2.0;
467  double cY = (forestOptions->bbox[4] - forestOptions->bbox[1]) / 2.0;
468 
469  double scalingFactor =
470  1.5; // The octree is this times bigger than the surface mesh's bounding box
471  double c = scalingFactor * fmax(cX, cY);
472 
473  // TODO : Compute any bounding box, not necessarily aligned with the axes
474  const double vertices[8 * 3] = {
475  centreX - c, centreY - c, 0.0, centreX + c, centreY - c, 0.0,
476  centreX - c, centreY + c, 0.0, centreX + c, centreY + c, 0.0,
477  centreX - c, centreY - c, 0.0, centreX + c, centreY - c, 0.0,
478  centreX - c, centreY + c, 0.0, centreX + c, centreY + c, 0.0,
479  };
480  const p4est_topidx_t tree_to_vertex[1 * 8] = {0, 1, 2, 3, 4, 5, 6, 7};
481  const p4est_topidx_t tree_to_tree[1 * 6] = {0, 0, 0, 0, 0, 0};
482  const int8_t tree_to_face[1 * 6] = {0, 1, 2, 3, 4, 5};
483 
484  return p4est_connectivity_new_copy(num_vertices, num_trees, 0, 0, vertices,
485  tree_to_vertex, tree_to_tree, tree_to_face,
486  nullptr, &num_ett, nullptr, nullptr,
487  nullptr, &num_ctt, nullptr, nullptr);
488 }
489 static inline double bulkSize(double x, double y, double z, double hBulk)
490 {
491  return hBulk;
492 }
493 
494 /* Fills xyz[] with the coordinates of the center of the given tree cell. */
495 static inline void getCellCenter(p4est_t *p4est, p4est_topidx_t which_tree,
496  p4est_quadrant_t *q, double xyz[3])
497 {
498  p4est_qcoord_t half_length = P4EST_QUADRANT_LEN(q->level) / 2;
499  p4est_qcoord_to_vertex(p4est->connectivity, which_tree, q->x + half_length,
500  q->y + half_length, q->z + half_length, xyz);
501 }
502 
503 /* Fills min[] & max[] with the coordinates of the cell viewed as a bounding
504  * box. */
505 static inline void getCellBBox(p4est_t *p4est, p4est_topidx_t which_tree,
506  p4est_quadrant_t *q, double min[3],
507  double max[3])
508 {
509  p4est_qcoord_t length = P4EST_QUADRANT_LEN(q->level);
510  p4est_qcoord_to_vertex(p4est->connectivity, which_tree, q->x, q->y, q->z,
511  min);
512  p4est_qcoord_to_vertex(p4est->connectivity, which_tree, q->x + length,
513  q->y + length, q->z + length, max);
514 }
515 
516 /* Fills h with the dimension of the given cell (length of a cell edge). */
517 static void getCellSize(p4est_t *p4est, p4est_topidx_t which_tree,
518  p4est_quadrant_t *q, double *h)
519 {
520  double min[3], max[3];
521  p4est_qcoord_t length = P4EST_QUADRANT_LEN(q->level);
522  p4est_qcoord_to_vertex(p4est->connectivity, which_tree, q->x, q->y, q->z,
523  min);
524  p4est_qcoord_to_vertex(p4est->connectivity, which_tree, q->x + length,
525  q->y + length, q->z + length, max);
526  // All cell edges are supposed to be the same length h (-:
527  *h = fmax(max[0] - min[0], fmax(max[1] - min[1], max[2] - min[2]));
528 }
529 
530 /* Callback used by p4est to initialize the user_data on each tree cell. */
531 static inline void initializeCell(p4est_t *p4est, p4est_topidx_t which_tree,
532  p4est_quadrant_t *q)
533 {
534  ForestOptions *forestOptions = (ForestOptions *)p4est->user_pointer;
535  size_data_t *data = (size_data_t *)q->p.user_data;
536 
537  double center[3];
538  getCellCenter(p4est, which_tree, q, center);
539 
540  // Set cell size
541  data->size = forestOptions->sizeFunction(center[0], center[1], center[2],
542  forestOptions->hbulk);
543  // data->M = SMetric3(1. / (forestOptions->hbulk * forestOptions->hbulk));
544  data->M = SMetric3(1. / (0.4 * 0.4));
545 
546  // Set size gradient to zero
547  for(int i = 0; i < P4EST_DIM; ++i) data->ds[i] = 0.0;
548  // Set cell dimension (edge length)
549  getCellSize(p4est, which_tree, q, &(data->h));
550 
551  data->hasIntersection = false;
552 }
553 
554 /* Creates (allocates) the forestOptions structure. */
555 HXTStatus forestOptionsCreate(ForestOptions **forestOptions)
556 {
557  HXT_CHECK(hxtAlignedMalloc(forestOptions, sizeof(ForestOptions)));
558  if(*forestOptions == nullptr) return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
559  return HXT_STATUS_OK;
560 }
561 
562 /* Destroys forestOptions. */
563 HXTStatus forestOptionsDelete(ForestOptions **forestOptions)
564 {
565  HXT_CHECK(hxtFree(forestOptions));
566  return HXT_STATUS_OK;
567 }
568 
569 HXTStatus loadGlobalData(ForestOptions *forestOptions, const char *filename)
570 {
571  FILE *f = fopen(filename, "r");
572  char buf[BUFSIZ] = {""};
573  if(fgets(buf, BUFSIZ, f) == nullptr) {
574  return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
575  }
576  sscanf(buf, "%lf %lf", &forestOptions->hmin, &forestOptions->hmax);
577  fclose(f);
578  Msg::Info("Loaded global data from %s", filename);
579  Msg::Info("Min size = %f", forestOptions->hmin);
580  Msg::Info("Max size = %f", forestOptions->hmax);
581  return HXT_STATUS_OK;
582 }
583 
584 /* Creates a (sequential) forest structure by loading a p4est file. */
585 HXTStatus forestLoad(Forest **forest, const char *forestFile,
586  const char *dataFile, ForestOptions *forestOptions)
587 {
588  if(forestFile == nullptr) return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
589 
590  HXT_CHECK(hxtMalloc(forest, sizeof(Forest)));
591  if(*forest == nullptr) return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
592 
593  HXT_CHECK(loadGlobalData(forestOptions, dataFile));
594 
595  p4est_connectivity_t *connect;
596  sc_MPI_Comm mpicomm = sc_MPI_COMM_WORLD;
597  int load_data = true;
598  int autopartition = true;
599  int broadcasthead = true;
600 
601  (*forest)->p4est = p4est_load_ext(forestFile, mpicomm, sizeof(size_data_t),
602  load_data, autopartition, broadcasthead,
603  (void *)forestOptions, &connect);
604 
605  ForestOptions *fO = (ForestOptions *)(*forest)->p4est->user_pointer;
606  if(fO == nullptr) return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
607 
608  return HXT_STATUS_OK;
609 }
610 
611 /* Creates a (sequential) forest structure from the forestOptions information.
612  The forest is not refined; it consists of the root octant. */
613 HXTStatus forestCreate(int argc, char **argv, Forest **forest,
614  const char *filename, ForestOptions *forestOptions)
615 {
616  HXT_CHECK(hxtMalloc(forest, sizeof(Forest)));
617  if(*forest == nullptr) return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
618 
619  int mpiret;
620  sc_MPI_Comm mpicomm;
621  p4est_connectivity_t *connect;
622 
623  /* Initialize MPI; see sc_mpi.h.
624  * If configure --enable-mpi is given these are true MPI calls.
625  * Else these are dummy functions that simulate a single-processor run. */
626  mpiret = sc_MPI_Init(&argc, &argv);
627  SC_CHECK_MPI(mpiret);
628  mpicomm = sc_MPI_COMM_WORLD;
629 
630  /* These functions are optional. If called they store the MPI rank as a
631  * static variable so subsequent global p4est log messages are only issued
632  * from processor zero. Here we turn off most of the logging; see sc.h. */
633  // sc_init(mpicomm, 1, 1, NULL, SC_LP_ESSENTIAL);
634  // p4est_init(NULL, SC_LP_PRODUCTION);
635 
636  /* Create a connectivity from the bounding box */
637  if(forestOptions->dim == 3) {
638  connect = p8est_connectivity_new_cube(forestOptions);
639  }
640  else {
641  connect = p8est_connectivity_new_square(forestOptions);
642  }
643 
644  if(connect == nullptr) return HXT_ERROR(HXT_STATUS_OUT_OF_MEMORY);
645 
646  // #ifdef P4EST_WITH_METIS
647  // // Use metis (if p4est is compiled with the flag '--with-metis') to
648  // // * reorder the connectivity for better parititioning of the forest
649  // // * across processors.
650  // p4est_connectivity_reorder(mpicomm, 0, conn, P4EST_CONNECT_FACE);
651  // #endif /* P4EST_WITH_METIS */
652 
653  // Assign bulkSize callback if no sizeFunction is specified.
654  if(forestOptions->sizeFunction == nullptr)
655  forestOptions->sizeFunction = &bulkSize;
656 
657  (*forest)->p4est = p4est_new(mpicomm, connect, sizeof(size_data_t),
658  initializeCell, (void *)forestOptions);
659  (*forest)->forestOptions = forestOptions;
660 
661  return HXT_STATUS_OK;
662 }
663 
664 /* Deletes the forest structure. */
665 HXTStatus forestDelete(Forest **forest)
666 {
667  /* Destroy the p4est structure. */
668  p4est_connectivity_destroy((*forest)->p4est->connectivity);
669  p4est_destroy((*forest)->p4est);
670  /* Verify that allocations internal to p4est and sc do not leak memory.
671  * This should be called if sc_init () has been called earlier. */
672  // sc_finalize();
673  /* This is standard MPI programs. Without --enable-mpi, this is a dummy. */
674  int mpiret = sc_MPI_Finalize();
675  SC_CHECK_MPI(mpiret);
676 
677  HXT_CHECK(hxtFree(forest));
678 
679  return HXT_STATUS_OK;
680 }
681 
682 /* ========================================================================================================
683  FOREST REFINEMENT
684  ========================================================================================================
685  */
686 
687 /* Callback used by octreeRefineToBulkSize; returns 1 if the cells need
688  * refinement, 0 otherwise. */
689 static inline int refineToBulkSizeCallback(p4est_t *p4est,
690  p4est_topidx_t which_tree,
691  p4est_quadrant_t *q)
692 {
693  ForestOptions *forestOptions = (ForestOptions *)p4est->user_pointer;
694  size_data_t *data = (size_data_t *)q->p.user_data;
695  return data->h > forestOptions->hbulk;
696 }
697 
698 /* Used by curvatureRefine; returns 1 if the cell should be
699  refined according to the surface mesh curvature, 0 otherwise. */
700 static int curvatureRefineCallback(p4est_t *p4est, p4est_topidx_t which_tree,
701  p4est_quadrant_t *q)
702 {
703  ForestOptions *forestOptions = (ForestOptions *)p4est->user_pointer;
704  double h;
705  getCellSize(p4est, which_tree, q, &h);
706  double center[3];
707  getCellCenter(p4est, which_tree, q, center);
708 
709  if(forestOptions->dim == 3) {
710  double min[3], max[3];
711  getCellBBox(p4est, which_tree, q, min, max);
712 
713  std::vector<uint64_t> candidates;
714  forestOptions->triRTree->Search(min, max, rtreeCallback, &candidates);
715 
716  if(!candidates.empty()) {
717  double kmax = -1.e22; // To get min curvature size
718  double kmin = 1.e22;
719  double hf = DBL_MAX; // To get min feature size
720  for(auto tri = candidates.begin(); tri != candidates.end(); ++tri) {
721  for(int i = 0; i < 3; ++i) {
722  int node =
723  forestOptions->mesh->triangles.node[(size_t)3 * (*tri) + i];
724 
725  double *v1 = forestOptions->nodalCurvature + 6 * node;
726  double *v2 = forestOptions->nodalCurvature + 6 * node + 3;
727 
728  double k1, k2;
729  norme2(v1, &k1);
730  norme2(v2, &k2);
731 
732  kmax = fmax(kmax, fmax(k1, k2));
733  kmin = fmin(kmin, fmin(k1, k2));
734 
735  hf = fmin(hf, (*forestOptions->featureSizeAtVertices)[node]);
736  }
737  }
738 
739  // There is curvature and/or feature size : the cell size should
740  // accurately represent the surface and the feature size, so as not
741  // propagate small feature size far from the feature.
742  // Cell is refined according to the chosen density of nodes.
743  double hc = 2 * M_PI / (forestOptions->nodePerTwoPi * kmax);
744  double nElemPerCell = 2;
745 
746  if(h > fmin(hc / nElemPerCell, hf) && h >= forestOptions->hmin) {
747  return 1;
748  }
749  else {
750  return 0;
751  }
752  }
753  else { // candidates.empty()
754  // If the cell has no intersection with the surface mesh, it is not
755  // refined.
756  return 0;
757  }
758  }
759  else {
760  // 2D curvature refinement is still a work in progress
761  // double kmax = -1.e22;
762  // double x, y;
763  // bool inbox;
764 
765  // // Calculer kmax sur toutes les courbes
766  // int n = 10;
767  // double par[n];
768  // for(uint16_t i = 0; i < n; ++i) par[i] = i*2.*M_PI/n;
769 
770  // for(uint16_t i = 0; i < (*forestOptions->curvFunctions).size(); ++i){
771  // for(uint16_t j = 0; j < n; ++j){
772  // x = (*forestOptions->xFunctions)[i](par[j]);
773  // y = (*forestOptions->yFunctions)[i](par[j]);
774  // inbox = (x <= center[0] + h/2.) && (x >= center[0] - h/2.);
775  // inbox &= (y <= center[1] + h/2.) && (y >= center[1] - h/2.);
776  // if(inbox) kmax = fmax(kmax,
777  // fabs((*forestOptions->curvFunctions)[i](par[j])));
778  // }
779  // }
780 
781  // if(kmax > 0){
782  // if(kmax < 1e-3){
783  // return 0;
784  // } else{
785  // double hc = 2*M_PI/(forestOptions->nodePerTwoPi * kmax);
786  // int nElemPerCell = 1;
787  // if(h > hc/nElemPerCell && h >= forestOptions->hmin){
788  // return 1;
789  // } else{
790  // return 0;
791  // }
792  // }
793  // } else{
794  return 0;
795  // }
796  }
797 }
798 
799 static int coarsenCallback(p4est_t *p4est, p4est_topidx_t which_tree,
800  p4est_quadrant_t *children[])
801 {
802  ForestOptions *forestOptions = (ForestOptions *)p4est->user_pointer;
803 
804  for(int n = 0; n < P4EST_CHILDREN; ++n) {
805  size_data_t *data = (size_data_t *)children[n]->p.user_data;
806 
807  double min[3], max[3];
808  getCellBBox(p4est, which_tree, children[n], min, max);
809 
810  std::vector<uint64_t> candidates;
811  forestOptions->triRTree->Search(min, max, rtreeCallback, &candidates);
812 
813  // Cells are not merged if any one of them touches the surface mesh.
814  if(!candidates.empty()) return 0;
815  // Cells are not merged if the resulting cell size is > than hbulk.
816  if(2.0 * data->h > forestOptions->hbulk) return 0;
817  }
818 
819  return 1;
820 }
821 
822 static void assignSizeAfterRefinement(p4est_iter_volume_info_t *info,
823  void *user_data)
824 {
825  p4est_t *p4est = info->p4est;
826  p4est_quadrant_t *q = info->quad;
827  p4est_topidx_t which_tree = info->treeid;
828  size_data_t *data = (size_data_t *)q->p.user_data;
829  ForestOptions *forestOptions = (ForestOptions *)user_data;
830 
831  double h;
832  getCellSize(p4est, which_tree, q, &h);
833  double center[3];
834  getCellCenter(p4est, which_tree, q, center);
835 
836  if(forestOptions->dim == 3) {
837  double min[3], max[3];
838  getCellBBox(p4est, which_tree, q, min, max);
839 
840  std::vector<uint64_t> candidates;
841  forestOptions->triRTree->Search(min, max, rtreeCallback, &candidates);
842 
843  if(!candidates.empty()) {
844  double kmax = -1.0e22;
845  double hf = DBL_MAX;
846  for(auto tri = candidates.begin(); tri != candidates.end(); ++tri) {
847  for(int i = 0; i < 3; ++i) {
848  int node =
849  forestOptions->mesh->triangles.node[(size_t)3 * (*tri) + i];
850  double *v1 = forestOptions->nodalCurvature + 6 * node;
851  double *v2 = forestOptions->nodalCurvature + 6 * node + 3;
852  double k1, k2;
853  norme2(v1, &k1);
854  norme2(v2, &k2);
855  kmax = fmax(kmax, fmax(k1, k2));
856  hf = fmin(hf, (*forestOptions->featureSizeAtVertices)[node]);
857  }
858  }
859  data->size =
860  fmax(forestOptions->hmin,
861  fmin(forestOptions->hmax,
862  fmin(hf, 2 * M_PI / (forestOptions->nodePerTwoPi * kmax))));
863  // data->size = kmax;
864  }
865  else {
866  data->size =
867  fmax(forestOptions->hmin, fmin(forestOptions->hmax, data->size));
868  }
869  }
870  else {
871  // 2D
872  // double kmax = -1.e22;
873  // double x, y;
874  // bool inbox;
875 
876  // // Calculer kmax sur toutes les courbes
877  // int n = 10;
878  // double par[n];
879  // for(uint16_t i = 0; i < n; ++i) par[i] = i*2.*M_PI/n;
880 
881  // for(uint16_t i = 0; i < (*forestOptions->curvFunctions).size(); ++i){
882  // for(uint16_t j = 0; j < n; ++j){
883  // x = (*forestOptions->xFunctions)[i](par[j]);
884  // y = (*forestOptions->yFunctions)[i](par[j]);
885  // inbox = (x <= center[0] + h/2.) && (x >= center[0] - h/2.);
886  // inbox &= (y <= center[1] + h/2.) && (y >= center[1] - h/2.);
887  // if(inbox) kmax = fmax(kmax,
888  // fabs((*forestOptions->curvFunctions)[i](par[j])));
889  // }
890  // }
891 
892  // if(kmax > 0){
893  // double hc = 2*M_PI/(forestOptions->nodePerTwoPi * kmax);
894  // data->size = fmin(data->size, hc);
895  // }
896  }
897 }
898 
899 #if 0 // TODO remove once this is used
900 static void assignSizeAfterRefinementAniso(p4est_iter_volume_info_t * info, void *user_data){
901  p4est_t *p4est = info->p4est;
902  p4est_quadrant_t *q = info->quad;
903  p4est_topidx_t which_tree = info->treeid;
904  size_data_t *data = (size_data_t *) q->p.user_data;
905  ForestOptions *forestOptions = (ForestOptions *) user_data;
906 
907  double h;
908  getCellSize(p4est, which_tree, q, &h);
909  double center[3];
910  getCellCenter(p4est, which_tree, q, center);
911 
912  // SVector3 t1(1.0,0.0,0.0);
913  // SVector3 t2(0.0,1.0,0.0);
914  // SVector3 t3(0.0,0.0,1.0);
915  // double h3 = 0.1;
916  // double h2 = 0.1 + 3.*fabs(center[1]);
917  // double h1 = 0.1 + 60.*fabs(center[1]);
918  // data->M = SMetric3(1.0/(h1*h1), 1.0/(h2*h2), 1.0/(h3*h3), t1, t2, t3);
919 
920  if(forestOptions->dim == 3){
921  double min[3], max[3];
922  getCellBBox(p4est, which_tree, q, min, max);
923 
924  std::vector<uint64_t> candidates;
925  forestOptions->triRTree->Search(min, max, rtreeCallback, &candidates);
926 
927  SVector3 t1, t2, n;
928 
929  if(!candidates.empty()){
930  double kmax1 = -1.0e22;
931  double kmax2 = -1.0e22;
932  double hf = DBL_MAX;
933  for(auto tri = candidates.begin(); tri != candidates.end(); ++tri){
934  for(int i = 0; i < 3; ++i){
935  int node = forestOptions->mesh->triangles.node[(size_t) 3*(*tri)+i];
936  double *v1 = forestOptions->nodalCurvature + 6*node;
937  double *v2 = forestOptions->nodalCurvature + 6*node + 3;
938 
939  // double x = forestOptions->mesh->vertices.coord[(size_t) 4*node+0];
940  // double y = forestOptions->mesh->vertices.coord[(size_t) 4*node+1];
941  // double z = forestOptions->mesh->vertices.coord[(size_t) 4*node+2];
942 
943  t1 = SVector3(v1[0], v1[1], v1[2]);
944  t1 = t1.unit();
945  t2 = SVector3(v2[0], v2[1], v2[2]);
946  t2 = t2.unit();
947  n = crossprod(t1,t2);
948  n = n.unit();
949 
950  data->t1 = t1;
951  data->t2 = t2;
952  data->n = n;
953 
954  // fprintf(forestOptions->file2, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n",
955  // x, y, z, x+t1[0], y+t1[1], z+t1[2], 1.0, 1.0);
956  // fprintf(forestOptions->file2, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n",
957  // x, y, z, x+t2[0], y+t2[1], z+t2[2], 1.0, 1.0);
958  // fprintf(forestOptions->file2, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n",
959  // x, y, z, x+n[0], y+n[1], z+n[2], 1.0, 1.0);
960 
961  double k1, k2;
962  norme2(v1, &k1);
963  norme2(v2, &k2);
964  kmax1 = fmax(kmax1,k1);
965  kmax2 = fmax(kmax2,k2);
966  hf = fmin(hf, (*forestOptions->featureSizeAtVertices)[node]);
967  }
968  }
969  double hc1 = fmax(forestOptions->hmin, fmin(forestOptions->hmax, 2*M_PI/(forestOptions->nodePerTwoPi * kmax1)) );
970  double hc2 = fmax(forestOptions->hmin, fmin(forestOptions->hmax, 2*M_PI/(forestOptions->nodePerTwoPi/4.0 * kmax2)) );
971  hf = fmax(forestOptions->hmin, fmin(forestOptions->hmax, hf));
972  // SMetric3 m(1.0/(hf*hf), 1.0/(hc1*hc1), 1.0/(hc2*hc2), n, t1, t2);
973  SMetric3 m(1.0/(hf*hf), 1.0/(hc1*hc1), 1.0/(hc2*hc2), n, t1, t2);
974  data->M = m;
975  data->size = fmax(forestOptions->hmin, fmin(forestOptions->hmax, fmin(hf, 2*M_PI/(forestOptions->nodePerTwoPi * fmax(kmax1,kmax2) )) ));
976  }
977  else{
978  data->size = fmax(forestOptions->hmin, fmin(forestOptions->hmax, data->size));
979  }
980  } else{
981  // // 2D : To do
982  // double kmax = -1.e22;
983  // double x, y;
984  // bool inbox;
985 
986  // // Calculer kmax sur toutes les courbes
987  // int n = 10;
988  // double par[n];
989  // for(uint16_t i = 0; i < n; ++i) par[i] = i*2.*M_PI/n;
990 
991  // for(uint16_t i = 0; i < (*forestOptions->curvFunctions).size(); ++i){
992  // for(uint16_t j = 0; j < n; ++j){
993  // x = (*forestOptions->xFunctions)[i](par[j]);
994  // y = (*forestOptions->yFunctions)[i](par[j]);
995  // inbox = (x <= center[0] + h/2.) && (x >= center[0] - h/2.);
996  // inbox &= (y <= center[1] + h/2.) && (y >= center[1] - h/2.);
997  // if(inbox) kmax = fmax(kmax, fabs((*forestOptions->curvFunctions)[i](par[j])));
998  // }
999  // }
1000 
1001  // if(kmax > 0){
1002  // double hc = 2*M_PI/(forestOptions->nodePerTwoPi * kmax);
1003  // data->size = fmin(data->size, hc);
1004  // }
1005  }
1006 }
1007 #endif // TODO remove once this is used
1008 
1009 HXTStatus forestRefine(Forest *forest)
1010 {
1011  /* Refine recursively the tree until all cells' size is at most hbulk. */
1012  p4est_refine_ext(forest->p4est, 1, P4EST_QMAXLEVEL, refineToBulkSizeCallback,
1013  initializeCell, nullptr);
1014  // Refine recursively with respect to the curvature
1015  p4est_refine_ext(forest->p4est, 1, P4EST_QMAXLEVEL, curvatureRefineCallback,
1016  initializeCell, nullptr);
1017  // Coarsen
1018  p4est_coarsen_ext(forest->p4est, 1, 0, coarsenCallback, initializeCell,
1019  nullptr);
1020  // Balance the octree to get 2:1 ratio between adjacent cells
1021  p4est_balance_ext(forest->p4est, P4EST_CONNECT_FACE, initializeCell, nullptr);
1022  // Assign size on the new cells
1023 
1024  forest->forestOptions->file1 = fopen("nonNorme.pos", "w");
1025  if(forest->forestOptions->file1 == nullptr)
1026  return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
1027  forest->forestOptions->file2 = fopen("norme.pos", "w");
1028  if(forest->forestOptions->file2 == nullptr)
1029  return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
1030 
1031  fprintf(forest->forestOptions->file1, "View \"nonNorme\" {\n");
1032  fprintf(forest->forestOptions->file2, "View \"norme\" {\n");
1033 
1034  p4est_iterate(forest->p4est, nullptr, forest->forestOptions,
1035  assignSizeAfterRefinement, nullptr, nullptr, nullptr);
1036  // p4est_iterate(forest->p4est, NULL, forest->forestOptions,
1037  // assignSizeAfterRefinementAniso, NULL, NULL, NULL);
1038 
1039  fprintf(forest->forestOptions->file1, "};");
1040  fclose(forest->forestOptions->file1);
1041  fprintf(forest->forestOptions->file2, "};");
1042  fclose(forest->forestOptions->file2);
1043 
1044  return HXT_STATUS_OK;
1045 }
1046 
1047 /* ========================================================================================================
1048  SIZE GRADIENT COMPUTATION & SMOOTHING
1049  ========================================================================================================
1050  */
1051 static inline void resetCell(p4est_iter_volume_info_t *info, void *user_data)
1052 {
1053  size_data_t *data = (size_data_t *)info->quad->p.user_data;
1054  // Reset gradient.
1055  data->ds[0] = 0.0;
1056  data->ds[1] = 0.0;
1057  data->ds[2] = 0.0;
1058 }
1059 
1060 static void computeGradientCenter(p4est_iter_face_info_t *info, void *user_data)
1061 {
1062  p4est_iter_face_side_t *side[2];
1063  sc_array_t *sides = &(info->sides);
1064  size_data_t *data;
1065  size_data_t *data_opp;
1066  double s_avg, s_sum;
1067  int which_face_opp;
1068  // Index of current face on the opposite cell (0 if current is 1 and vice
1069  // versa).
1070  int iOpp;
1071 
1072  side[0] = p4est_iter_fside_array_index_int(sides, 0);
1073  side[1] = p4est_iter_fside_array_index_int(sides, 1);
1074 
1075  if(sides->elem_count == 2) {
1076  for(int i = 0; i < 2; i++) {
1077  iOpp = 1 - i;
1078  which_face_opp =
1079  side[iOpp]->face; /* 0,1 == -+x, 2,3 == -+y, 4,5 == -+z */
1080 
1081  // s_avg is the average size on the P4EST_HALF opposite cells
1082  s_avg = s_sum = 0;
1083 
1084  // Current cells are hanging
1085  if(side[i]->is_hanging) {
1086  for(int j = 0; j < P4EST_HALF; j++) {
1087  data = (size_data_t *)side[i]->is.hanging.quad[j]->p.user_data;
1088  s_avg += data->size;
1089  }
1090  s_sum = s_avg;
1091  s_avg /= P4EST_HALF;
1092 
1093  data_opp = (size_data_t *)side[iOpp]->is.full.quad->p.user_data;
1094 
1095  switch(which_face_opp) {
1096  case 0:
1097  data_opp->ds[0] -=
1098  0.5 * (s_avg - data_opp->size) / (data_opp->h / 2. + data->h / 2.);
1099  break;
1100  case 1:
1101  data_opp->ds[0] +=
1102  0.5 * (s_avg - data_opp->size) / (data_opp->h / 2. + data->h / 2.);
1103  break;
1104  case 2:
1105  data_opp->ds[1] -=
1106  0.5 * (s_avg - data_opp->size) / (data_opp->h / 2. + data->h / 2.);
1107  break;
1108  case 3:
1109  data_opp->ds[1] +=
1110  0.5 * (s_avg - data_opp->size) / (data_opp->h / 2. + data->h / 2.);
1111  break;
1112  case 4:
1113  data_opp->ds[2] -=
1114  0.5 * (s_avg - data_opp->size) / (data_opp->h / 2. + data->h / 2.);
1115  break;
1116  case 5:
1117  data_opp->ds[2] +=
1118  0.5 * (s_avg - data_opp->size) / (data_opp->h / 2. + data->h / 2.);
1119  break;
1120  }
1121  }
1122  // Current cell is full
1123  else {
1124  data = (size_data_t *)side[i]->is.full.quad->p.user_data;
1125  if(side[iOpp]->is_hanging) {
1126  // Current full - Opposite hanging
1127  double s_avg_opp = 0;
1128  for(int j = 0; j < P4EST_HALF; ++j) {
1129  data_opp =
1130  (size_data_t *)side[iOpp]->is.hanging.quad[j]->p.user_data;
1131  s_avg_opp += data_opp->size;
1132  }
1133  s_avg_opp /= P4EST_HALF;
1134 
1135  for(int j = 0; j < P4EST_HALF; ++j) {
1136  data_opp =
1137  (size_data_t *)side[iOpp]->is.hanging.quad[j]->p.user_data;
1138  switch(which_face_opp) {
1139  case 0:
1140  data_opp->ds[0] -= 0.5 * (data->size - data_opp->size) /
1141  (data_opp->h / 2. + data->h / 2.);
1142  break;
1143  case 1:
1144  data_opp->ds[0] += 0.5 * (data->size - data_opp->size) /
1145  (data_opp->h / 2. + data->h / 2.);
1146  break;
1147  case 2:
1148  data_opp->ds[1] -= 0.5 * (data->size - data_opp->size) /
1149  (data_opp->h / 2. + data->h / 2.);
1150  break;
1151  case 3:
1152  data_opp->ds[1] += 0.5 * (data->size - data_opp->size) /
1153  (data_opp->h / 2. + data->h / 2.);
1154  break;
1155  case 4:
1156  data_opp->ds[2] -= 0.5 * (data->size - data_opp->size) /
1157  (data_opp->h / 2. + data->h / 2.);
1158  break;
1159  case 5:
1160  data_opp->ds[2] += 0.5 * (data->size - data_opp->size) /
1161  (data_opp->h / 2. + data->h / 2.);
1162  break;
1163  }
1164  }
1165  }
1166  else {
1167  // Current full - Opposite full
1168  data_opp = (size_data_t *)side[iOpp]->is.full.quad->p.user_data;
1169  switch(which_face_opp) {
1170  case 0:
1171  data_opp->ds[0] -= 0.5 * (data->size - data_opp->size) /
1172  (data_opp->h / 2. + data->h / 2.);
1173  break;
1174  case 1:
1175  data_opp->ds[0] += 0.5 * (data->size - data_opp->size) /
1176  (data_opp->h / 2. + data->h / 2.);
1177  break;
1178  case 2:
1179  data_opp->ds[1] -= 0.5 * (data->size - data_opp->size) /
1180  (data_opp->h / 2. + data->h / 2.);
1181  break;
1182  case 3:
1183  data_opp->ds[1] += 0.5 * (data->size - data_opp->size) /
1184  (data_opp->h / 2. + data->h / 2.);
1185  break;
1186  case 4:
1187  data_opp->ds[2] -= 0.5 * (data->size - data_opp->size) /
1188  (data_opp->h / 2. + data->h / 2.);
1189  break;
1190  case 5:
1191  data_opp->ds[2] += 0.5 * (data->size - data_opp->size) /
1192  (data_opp->h / 2. + data->h / 2.);
1193  break;
1194  }
1195  }
1196  }
1197  }
1198  }
1199  // Nothing to do on the boundaries (0 flux)
1200 }
1201 
1202 HXTStatus forestComputeGradient(Forest *forest)
1203 {
1204  // Iterate on each cell to reset size gradient and half lengths.
1205  p4est_iterate(forest->p4est, nullptr, nullptr, resetCell, nullptr, nullptr,
1206  nullptr);
1207  // Compute gradient at cell center using finite differences
1208  p4est_iterate(forest->p4est, nullptr, nullptr, nullptr, computeGradientCenter,
1209  nullptr, nullptr);
1210 
1211  return HXT_STATUS_OK;
1212 }
1213 
1214 static inline void getMaxGradient(p4est_iter_volume_info_t *info,
1215  void *user_data)
1216 {
1217  p4est_quadrant_t *q = info->quad;
1218  size_data_t *data = (size_data_t *)q->p.user_data;
1219  double *gradMax = static_cast<double *>(user_data);
1220  gradMax[0] = SC_MAX(fabs(data->ds[0]), gradMax[0]);
1221  gradMax[1] = SC_MAX(fabs(data->ds[1]), gradMax[1]);
1222  gradMax[2] = SC_MAX(fabs(data->ds[2]), gradMax[2]);
1223 }
1224 
1225 static inline void getMinSize(p4est_iter_volume_info_t *info, void *user_data)
1226 {
1227  p4est_quadrant_t *q = info->quad;
1228  size_data_t *data = (size_data_t *)q->p.user_data;
1229  double *minsize = (double *)user_data;
1230  *minsize = fmin(*minsize, data->size);
1231 }
1232 
1233 HXTStatus forestGetMaxGradient(Forest *forest, double *gradMax)
1234 {
1235  p4est_iterate(forest->p4est, nullptr, (void *)gradMax, getMaxGradient,
1236  nullptr, nullptr, nullptr);
1237  return HXT_STATUS_OK;
1238 }
1239 
1240 HXTStatus forestGetMinSize(Forest *forest, double *minsize)
1241 {
1242  double minSize = 1e22;
1243  p4est_iterate(forest->p4est, nullptr, (void *)&minSize, getMinSize, nullptr,
1244  nullptr, nullptr);
1245  *minsize = minSize;
1246  return HXT_STATUS_OK;
1247 }
1248 
1249 static void limitSize(p4est_iter_face_info_t *info, void *user_data)
1250 {
1251  p4est_iter_face_side_t *side[2];
1252  sc_array_t *sides = &(info->sides);
1253  size_data_t *data;
1254  size_data_t *data_opp1;
1255  int which_dir;
1256  int which_face_opp;
1257  int iOpp;
1258 
1259  ForestOptions *forestOptions = (ForestOptions *)user_data;
1260  double alpha = forestOptions->gradation - 1.0;
1261  double tol = 1e-3;
1262 
1263  side[0] = p4est_iter_fside_array_index_int(sides, 0);
1264  side[1] = p4est_iter_fside_array_index_int(sides, 1);
1265 
1266  if(sides->elem_count == 2) {
1267  for(int i = 0; i < 2; ++i) {
1268  iOpp = 1 - i;
1269  which_dir = side[i]->face / 2; // Direction x (0), y (1) ou z(2)
1270  which_face_opp = side[iOpp]->face;
1271 
1272  if(side[i]->is_hanging) {
1273  // Current hanging - Opposes full
1274  data_opp1 = (size_data_t *)side[iOpp]->is.full.quad->p.user_data;
1275 
1276  for(int j = 0; j < P4EST_HALF; ++j) {
1277  data = (size_data_t *)side[i]->is.hanging.quad[j]->p.user_data;
1278 
1279  if(fabs(data->ds[which_dir]) > alpha + tol) {
1280  if(data->size > data_opp1->size) {
1281  switch(which_face_opp) {
1282  case 0:
1283  data->size = fmin(
1284  data->size, data_opp1->size +
1285  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1286  break;
1287  case 1:
1288  data->size = fmin(
1289  data->size, data_opp1->size +
1290  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1291  break;
1292  case 2:
1293  data->size = fmin(
1294  data->size, data_opp1->size +
1295  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1296  break;
1297  case 3:
1298  data->size = fmin(
1299  data->size, data_opp1->size +
1300  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1301  break;
1302  case 4:
1303  data->size = fmin(
1304  data->size, data_opp1->size +
1305  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1306  break;
1307  case 5:
1308  data->size = fmin(
1309  data->size, data_opp1->size +
1310  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1311  break;
1312  }
1313  }
1314  else {
1315  switch(which_face_opp) {
1316  case 0:
1317  data_opp1->size = fmin(
1318  data_opp1->size,
1319  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1320  break;
1321  case 1:
1322  data_opp1->size = fmin(
1323  data_opp1->size,
1324  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1325  break;
1326  case 2:
1327  data_opp1->size = fmin(
1328  data_opp1->size,
1329  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1330  break;
1331  case 3:
1332  data_opp1->size = fmin(
1333  data_opp1->size,
1334  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1335  break;
1336  case 4:
1337  data_opp1->size = fmin(
1338  data_opp1->size,
1339  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1340  break;
1341  case 5:
1342  data_opp1->size = fmin(
1343  data_opp1->size,
1344  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1345  break;
1346  }
1347  }
1348  } // if ds > alpha-1
1349  } // for j hanging
1350  } // if hanging
1351  else {
1352  data = (size_data_t *)side[i]->is.full.quad->p.user_data;
1353 
1354  if(fabs(data->ds[which_dir]) > alpha + tol) {
1355  if(side[iOpp]->is_hanging) {
1356  // Current full - Oppose hanging
1357  for(int j = 0; j < P4EST_HALF; ++j) {
1358  data_opp1 =
1359  (size_data_t *)side[iOpp]->is.hanging.quad[j]->p.user_data;
1360  if(data->size > data_opp1->size) {
1361  switch(which_face_opp) {
1362  case 0:
1363  data->size = fmin(
1364  data->size, data_opp1->size +
1365  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1366  break;
1367  case 1:
1368  data->size = fmin(
1369  data->size, data_opp1->size +
1370  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1371  break;
1372  case 2:
1373  data->size = fmin(
1374  data->size, data_opp1->size +
1375  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1376  break;
1377  case 3:
1378  data->size = fmin(
1379  data->size, data_opp1->size +
1380  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1381  break;
1382  case 4:
1383  data->size = fmin(
1384  data->size, data_opp1->size +
1385  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1386  break;
1387  case 5:
1388  data->size = fmin(
1389  data->size, data_opp1->size +
1390  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1391  break;
1392  }
1393  }
1394  else {
1395  switch(which_face_opp) {
1396  case 0:
1397  data_opp1->size = fmin(
1398  data_opp1->size,
1399  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1400  break;
1401  case 1:
1402  data_opp1->size = fmin(
1403  data_opp1->size,
1404  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1405  break;
1406  case 2:
1407  data_opp1->size = fmin(
1408  data_opp1->size,
1409  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1410  break;
1411  case 3:
1412  data_opp1->size = fmin(
1413  data_opp1->size,
1414  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1415  break;
1416  case 4:
1417  data_opp1->size = fmin(
1418  data_opp1->size,
1419  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1420  break;
1421  case 5:
1422  data_opp1->size = fmin(
1423  data_opp1->size,
1424  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1425  break;
1426  }
1427  }
1428  }
1429  }
1430  else {
1431  // Current full - Oppose full
1432  data_opp1 = (size_data_t *)side[iOpp]->is.full.quad->p.user_data;
1433  if(data->size > data_opp1->size) {
1434  switch(which_face_opp) {
1435  case 0:
1436  data->size = fmin(
1437  data->size, data_opp1->size +
1438  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1439  break;
1440  case 1:
1441  data->size = fmin(
1442  data->size, data_opp1->size +
1443  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1444  break;
1445  case 2:
1446  data->size = fmin(
1447  data->size, data_opp1->size +
1448  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1449  break;
1450  case 3:
1451  data->size = fmin(
1452  data->size, data_opp1->size +
1453  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1454  break;
1455  case 4:
1456  data->size = fmin(
1457  data->size, data_opp1->size +
1458  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1459  break;
1460  case 5:
1461  data->size = fmin(
1462  data->size, data_opp1->size +
1463  (alpha) * (data_opp1->h / 2. + data->h / 2.));
1464  break;
1465  }
1466  }
1467  else {
1468  switch(which_face_opp) {
1469  case 0:
1470  data_opp1->size = fmin(
1471  data_opp1->size,
1472  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1473  break;
1474  case 1:
1475  data_opp1->size = fmin(
1476  data_opp1->size,
1477  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1478  break;
1479  case 2:
1480  data_opp1->size = fmin(
1481  data_opp1->size,
1482  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1483  break;
1484  case 3:
1485  data_opp1->size = fmin(
1486  data_opp1->size,
1487  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1488  break;
1489  case 4:
1490  data_opp1->size = fmin(
1491  data_opp1->size,
1492  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1493  break;
1494  case 5:
1495  data_opp1->size = fmin(
1496  data_opp1->size,
1497  data->size + (alpha) * (data_opp1->h / 2. + data->h / 2.));
1498  break;
1499  }
1500  }
1501  }
1502  } // if gradient trop grand
1503  } // else
1504  }
1505  }
1506 }
1507 
1508 HXTStatus forestLimitSize(Forest *forest)
1509 {
1510  p4est_iterate(forest->p4est, nullptr, (void *)forest->forestOptions, nullptr,
1511  limitSize, nullptr, nullptr);
1512  return HXT_STATUS_OK;
1513 }
1514 
1515 HXTStatus forestSizeSmoothing(Forest *forest)
1516 {
1517  double gradMax[3], tol = 1e-3, gradLinf = 1e22;
1518  // double minSize = DBL_MAX;
1519  int iter = 0, nmax = 100;
1520 
1521  // HXT_CHECK(forestGetMinSize(forest, &minSize));
1522  // printf("Minsize before smoothing = %f\n", minSize);
1523 
1524  while(iter++ < nmax &&
1525  gradLinf > tol + forest->forestOptions->gradation - 1.0) {
1526  HXT_CHECK(forestComputeGradient(forest)); // Compute gradient
1527  HXT_CHECK(forestLimitSize(
1528  forest)); // Smooth large sizes to limit gradient to foresOptions->gradMax
1529  for(int i = 0; i < 3; ++i) gradMax[i] = -1e22; // Stopping criterion
1530  HXT_CHECK(forestGetMaxGradient(forest, gradMax));
1531  gradLinf = fmax(fabs(gradMax[0]), fmax(fabs(gradMax[1]), fabs(gradMax[2])));
1532  // HXT_CHECK(forestGetMinSize(forest, &minSize));
1533  // printf("Minsize = %f\n", minSize);
1534  }
1535  Msg::Info("Max gradient after smoothing : (%f - %f - %f)", fabs(gradMax[0]),
1536  fabs(gradMax[1]), fabs(gradMax[2]));
1537  return HXT_STATUS_OK;
1538 }
1539 
1540 /* ========================================================================================================
1541  SEARCH AND REPLACE
1542  ========================================================================================================
1543  */
1544 
1545 static int searchAndAssignConstant(p4est_t *p4est, p4est_topidx_t which_tree,
1546  p4est_quadrant_t *q,
1547  p4est_locidx_t local_num, void *point)
1548 {
1549  bool in_box, is_leaf = local_num >= 0;
1550  size_data_t *data = (size_data_t *)q->p.user_data;
1551  size_point_t *p = (size_point_t *)point;
1552  ForestOptions *forestOptions = (ForestOptions *)p4est->user_pointer;
1553  // We have to recompute the cell dimension h for the root (non leaves) octants
1554  // because it seems to be undefined. Otherwise it's contained in
1555  // q->p.user_data.
1556  double h, center[3];
1557  if(!is_leaf)
1558  getCellSize(p4est, which_tree, q, &h);
1559  else
1560  h = data->h;
1561  getCellCenter(p4est, which_tree, q, center);
1562 
1563  // double epsilon = 1e-13;
1564  // in_box = (p->x < center[0] + h/2. + epsilon) && (p->x > center[0] - h/2. -
1565  // epsilon); in_box &= (p->y < center[1] + h/2. + epsilon) && (p->y >
1566  // center[1] - h/2. - epsilon); in_box &= (p->z < center[2] + h/2. + epsilon)
1567  // && (p->z > center[2] - h/2. - epsilon);
1568  in_box = (p->x <= center[0] + h / 2.) && (p->x >= center[0] - h / 2.);
1569  in_box &= (p->y <= center[1] + h / 2.) && (p->y >= center[1] - h / 2.);
1570  in_box &= (p->z <= center[2] + h / 2.) && (p->z >= center[2] - h / 2.);
1571 
1572  // A point can be on the exact boundary of two cells, hence we take the min.
1573  if(in_box && is_leaf) {
1574  p->size = fmin(p->size, data->size);
1575  p->size = fmax(forestOptions->hmin, fmin(forestOptions->hmax, p->size));
1576  p->isFound = true;
1577  }
1578 
1579  return in_box;
1580 }
1581 
1582 static int searchAndAssignConstantAniso(p4est_t *p4est,
1583  p4est_topidx_t which_tree,
1584  p4est_quadrant_t *q,
1585  p4est_locidx_t local_num, void *point)
1586 {
1587  bool in_box, is_leaf = local_num >= 0;
1588  size_data_t *data = (size_data_t *)q->p.user_data;
1589  size_point_t *p = (size_point_t *)point;
1590  // ForestOptions *forestOptions = (ForestOptions *) p4est->user_pointer;
1591  // We have to recompute the cell dimension h for the root (non leaves) octants
1592  // because it seems to be undefined. Otherwise it's contained in
1593  // q->p.user_data.
1594  double h, center[3];
1595  if(!is_leaf)
1596  getCellSize(p4est, which_tree, q, &h);
1597  else
1598  h = data->h;
1599  getCellCenter(p4est, which_tree, q, center);
1600 
1601  double epsilon = 1e-13;
1602  in_box = (p->x < center[0] + h / 2. + epsilon) &&
1603  (p->x > center[0] - h / 2. - epsilon);
1604  in_box &= (p->y < center[1] + h / 2. + epsilon) &&
1605  (p->y > center[1] - h / 2. - epsilon);
1606  in_box &= (p->z < center[2] + h / 2. + epsilon) &&
1607  (p->z > center[2] - h / 2. - epsilon);
1608  // in_box = (p->x <= center[0] + h/2.) && (p->x >= center[0] - h/2.);
1609  // in_box &= (p->y <= center[1] + h/2.) && (p->y >= center[1] - h/2.);
1610  // in_box &= (p->z <= center[2] + h/2.) && (p->z >= center[2] - h/2.);
1611 
1612  // A point can be on the exact boundary of two cells, hence we take the min.
1613  if(in_box && is_leaf) {
1614  // p->size = fmin(p->size, data->size);
1615  // p->size = fmax(forestOptions->hmin, fmin(forestOptions->hmax, p->size) );
1616  p->m = data->M;
1617  // p->m.print("p->m");
1618  p->isFound = true;
1619  }
1620 
1621  return in_box;
1622 }
1623 
1624 static int searchAndAssignLinear(p4est_t *p4est, p4est_topidx_t which_tree,
1625  p4est_quadrant_t *q, p4est_locidx_t local_num,
1626  void *point)
1627 {
1628  bool in_box = false, is_leaf = local_num >= 0;
1629  size_data_t *data = (size_data_t *)q->p.user_data;
1630  size_point_t *p = (size_point_t *)point;
1631  ForestOptions *forestOptions = (ForestOptions *)p4est->user_pointer;
1632  // We have to recompute the cell dimension h for the root (non leaves) octants
1633  // because it seems to be undefined. Otherwise it's contained in
1634  // q->p.user_data.
1635  double h, center[3];
1636  // if(!is_leaf) getCellSize(p4est, which_tree, q, &h);
1637  // else h = data->h;
1638  getCellSize(p4est, which_tree, q, &h);
1639  getCellCenter(p4est, which_tree, q, center);
1640 
1641  double epsilon = 1e-10;
1642  // in_box = (p->x < center[0] + h/2. + epsilon) && (p->x > center[0] - h/2. -
1643  // epsilon); in_box &= (p->y < center[1] + h/2. + epsilon) && (p->y >
1644  // center[1] - h/2. - epsilon); in_box &= (p->z < center[2] + h/2. + epsilon)
1645  // && (p->z > center[2] - h/2. - epsilon);
1646 
1647  // This misses some points...
1648  // in_box = (p->x <= center[0] + h/2.) && (p->x >= center[0] - h/2.);
1649  // in_box &= (p->y <= center[1] + h/2.) && (p->y >= center[1] - h/2.);
1650  // in_box &= (p->z <= center[2] + h/2.) && (p->z >= center[2] - h/2.);
1651 
1652  SPoint3 C(center), P(p->x, p->y, p->z);
1653  SVector3 dir(C, P);
1654  SVector3 dx(1., 0., 0.);
1655  SVector3 dy(0., 1., 0.);
1656  SVector3 dz(0., 0., 1.);
1657 
1658  in_box = fabs(dot(dir, dx)) <= (h / 2. + epsilon) &&
1659  fabs(dot(dir, dy)) <= (h / 2. + epsilon) &&
1660  fabs(dot(dir, dz)) <= (h / 2. + epsilon);
1661 
1662  // A point can be on the exact boundary of two cells, hence we take the min.
1663  if(in_box && is_leaf) {
1664  p->size = fmin(p->size, data->size + data->ds[0] * (p->x - center[0]) +
1665  data->ds[1] * (p->y - center[1]) +
1666  data->ds[2] * (p->z - center[2]));
1667  p->size = fmax(forestOptions->hmin, fmin(forestOptions->hmax, p->size));
1668  p->isFound = true;
1669  }
1670 
1671  p->parcourus++;
1672 
1673  // if(is_leaf && !p->isFound && fabs(p->x) < 1e-4){
1674  // if(in_box && !p->isFound && fabs(p->x) < 1e-4){
1675  // printf("Point (%4.16e,%4.16e,%4.16e)\n", p->x, p->y, p->z);
1676  // printf("h = %4.16e, center = (%4.16e, %4.16e, %4.16e)\n", h, center[0],
1677  // center[1], center[2]); printf("%4.16e <= x <= %4.16e (%d)\n",
1678  // center[0]-h/2., center[0]+h/2., (p->x <= center[0] + h/2.) && (p->x >=
1679  // center[0] - h/2.) ); printf("%4.16e <= y <= %4.16e (%d)\n",
1680  // center[1]-h/2., center[1]+h/2., (p->y <= center[1] + h/2.) && (p->y >=
1681  // center[1] - h/2.) ); printf("%4.16e <= z <= %4.16e (%d)\n",
1682  // center[2]-h/2., center[2]+h/2., (p->z <= center[2] + h/2.) && (p->z >=
1683  // center[2] - h/2.) );
1684  // }
1685 
1686  return in_box;
1687 }
1688 
1689 // static int replace(p4est_t * p4est, p4est_topidx_t which_tree,
1690 // p4est_quadrant_t * q, p4est_locidx_t local_num, void *point){
1691 // bool in_box, is_leaf = local_num >= 0;
1692 // size_data_t *data = (size_data_t *) q->p.user_data;
1693 // size_point_t *p = (size_point_t *) point;
1694 // double h, center[3];
1695 // if(!is_leaf) getCellSize(p4est, which_tree, q, &h);
1696 // else h = data->h;
1697 // getCellCenter(p4est, which_tree, q, center);
1698 
1699 // in_box = (p->x <= center[0] + h/2.) && (p->x >= center[0] - h/2.);
1700 // in_box &= (p->y <= center[1] + h/2.) && (p->y >= center[1] - h/2.);
1701 // in_box &= (p->z <= center[2] + h/2.) && (p->z >= center[2] - h/2.);
1702 
1703 // if(in_box && is_leaf){
1704 // data->size = fmin(data->size, p->size);
1705 // }
1706 
1707 // return in_box;
1708 // }
1709 
1710 /* Search for a single point in the tree structure and returns its size.
1711  See searchAndAssign for the detailed comments. */
1712 HXTStatus forestSearchOne(Forest *forest, double x, double y, double z,
1713  double *size, int linear)
1714 {
1715  sc_array_t *points = sc_array_new_size(sizeof(size_point_t), 1);
1716 
1717  size_point_t *p = (size_point_t *)sc_array_index(points, 0);
1718  p->x = x;
1719  p->y = y;
1720  p->z = z;
1721  p->size = 1.0e22;
1722  p->isFound = false;
1723  p->parcourus = 0;
1724 
1725  if(linear) {
1726  p4est_search(forest->p4est, nullptr, searchAndAssignLinear, points);
1727  }
1728  else {
1729  p4est_search(forest->p4est, nullptr, searchAndAssignConstant, points);
1730  }
1731 
1732  if(!p->isFound) {
1733  Msg::Info("(%+.4f,%+.4f,%+.4f) was not found in the meshsize field 8-|", x,
1734  y, z);
1735  Msg::Info("Octants parcourus : %d\n", p->parcourus);
1736  }
1737  *size = p->size;
1738 
1739  sc_array_destroy(points);
1740 
1741  return HXT_STATUS_OK;
1742 }
1743 
1744 HXTStatus forestSearchOneAniso(Forest *forest, double x, double y, double z,
1745  SMetric3 &m, int linear)
1746 {
1747  sc_array_t *points = sc_array_new_size(sizeof(size_point_t), 1);
1748 
1749  size_point_t *p = (size_point_t *)sc_array_index(points, 0);
1750  p->x = x;
1751  p->y = y;
1752  p->z = z;
1753  p->m = SMetric3(1.0);
1754  p->isFound = false;
1755 
1756  p4est_search(forest->p4est, nullptr, searchAndAssignConstantAniso, points);
1757 
1758  if(!p->isFound)
1759  Msg::Info("Point (%f,%f,%f) n'a pas été trouvé dans l'octree 8-|", x, y, z);
1760  else
1761  m = p->m;
1762 
1763  // m.print("apres search");
1764 
1765  sc_array_destroy(points);
1766 
1767  return HXT_STATUS_OK;
1768 }
1769 
1770 // Not finished
1771 HXTStatus hxtForestSearch(Forest *forest, std::vector<double> *x,
1772  std::vector<double> *y, std::vector<double> *z,
1773  std::vector<double> *size)
1774 {
1775  // Array of size_point_t to search in the tree
1776  sc_array_t *points = sc_array_new_size(sizeof(size_point_t), (*x).size());
1777  size_point_t *p;
1778 
1779  for(size_t i = 0; i < (*x).size(); ++i) {
1780  p = (size_point_t *)sc_array_index(points, i);
1781  p->x = (*x)[i];
1782  p->y = (*y)[i];
1783  p->z = (*z)[i];
1784  p->size = -1.0;
1785  }
1786 
1787  // Search on all cells
1788  p4est_search(forest->p4est, nullptr, searchAndAssignLinear, points);
1789 
1790  // Get the sizes
1791  for(size_t i = 0; i < (*x).size(); ++i) {
1792  p = (size_point_t *)sc_array_index(points, i);
1793  (*size)[i] = p->size;
1794  }
1795 
1796  // Clean up
1797  sc_array_destroy(points);
1798 
1799  return HXT_STATUS_OK;
1800 }
1801 
1802 /* ========================================================================================================
1803  CLOSE SURFACES DETECTION
1804  ========================================================================================================
1805  */
1806 
1807 // To quickly sort 4 integers
1808 static void sort4(int *d)
1809 {
1810 #define SWAP(x, y) \
1811  if(d[y] < d[x]) { \
1812  int tmp = d[x]; \
1813  d[x] = d[y]; \
1814  d[y] = tmp; \
1815  }
1816  SWAP(0, 1);
1817  SWAP(2, 3);
1818  SWAP(0, 2);
1819  SWAP(1, 3);
1820  SWAP(1, 2);
1821 #undef SWAP
1822 }
1823 
1824 // Checks whether two MTetrahedra have a common face without creating MFaces
1825 // (slow) Returns the index of the common face in t1 if any, -1 otherwise
1826 static int commonFaceTetFast(MTetrahedron *t1, MTetrahedron *t2)
1827 {
1828  int t10 = t1->getVertex(0)->getNum();
1829  int t11 = t1->getVertex(1)->getNum();
1830  int t12 = t1->getVertex(2)->getNum();
1831  int t13 = t1->getVertex(3)->getNum();
1832  int t20 = t2->getVertex(0)->getNum();
1833  int t21 = t2->getVertex(1)->getNum();
1834  int t22 = t2->getVertex(2)->getNum();
1835  int t23 = t2->getVertex(3)->getNum();
1836 
1837  bool b0 = (t10 == t20) || (t10 == t21) || (t10 == t22) || (t10 == t23);
1838  bool b1 = (t11 == t20) || (t11 == t21) || (t11 == t22) || (t11 == t23);
1839  bool b2 = (t12 == t20) || (t12 == t21) || (t12 == t22) || (t12 == t23);
1840  bool b3 = (t13 == t20) || (t13 == t21) || (t13 == t22) || (t13 == t23);
1841 
1842  if(b0 + b1 + b2 + b3 < 3) { return -1; }
1843  else {
1844  int v1[4] = {t10, t11, t12, t13};
1845  int v1cpy[4] = {t10, t11, t12, t13};
1846  int v2[4] = {t20, t21, t22, t23};
1847  sort4(v1);
1848  sort4(v2);
1849  t10 = v1[0];
1850  t11 = v1[1];
1851  t12 = v1[2];
1852  t13 = v1[3];
1853  t20 = v2[0];
1854  t21 = v2[1];
1855  t22 = v2[2];
1856  t23 = v2[3];
1857 
1858  bool b00 = (t11 == t21) && (t12 == t22) && (t13 == t23);
1859  bool b01 = (t11 == t20) && (t12 == t22) && (t13 == t23);
1860  bool b02 = (t11 == t20) && (t12 == t21) && (t13 == t23);
1861  bool b03 = (t11 == t20) && (t12 == t21) && (t13 == t22);
1862 
1863  bool b10 = (t10 == t21) && (t12 == t22) && (t13 == t23);
1864  bool b11 = (t10 == t20) && (t12 == t22) && (t13 == t23);
1865  bool b12 = (t10 == t20) && (t12 == t21) && (t13 == t23);
1866  bool b13 = (t10 == t20) && (t12 == t21) && (t13 == t22);
1867 
1868  bool b20 = (t10 == t21) && (t11 == t22) && (t12 == t23);
1869  bool b21 = (t10 == t20) && (t11 == t22) && (t12 == t23);
1870  bool b22 = (t10 == t20) && (t11 == t21) && (t12 == t23);
1871  bool b23 = (t10 == t20) && (t11 == t21) && (t12 == t22);
1872 
1873  bool b30 = (t10 == t21) && (t11 == t22) && (t13 == t23);
1874  bool b31 = (t10 == t20) && (t11 == t22) && (t13 == t23);
1875  bool b32 = (t10 == t20) && (t11 == t21) && (t13 == t23);
1876  bool b33 = (t10 == t20) && (t11 == t21) && (t13 == t22);
1877 
1878  int missing = -1; // The vertex that is missing from the common face
1879  if(b00 || b01 || b02 || b03)
1880  missing = 0; // return 3;
1881  else if(b10 || b11 || b12 || b13)
1882  missing = 1; // return 2;
1883  else if(b20 || b21 || b22 || b23)
1884  missing = 3; // return 0;
1885  else if(b30 || b31 || b32 || b33)
1886  missing = 2; // return 1;
1887 
1888  if(missing >= 0) {
1889  if(v1cpy[0] == v1[missing])
1890  return 3;
1891  else if(v1cpy[1] == v1[missing])
1892  return 2;
1893  else if(v1cpy[2] == v1[missing])
1894  return 1;
1895  else if(v1cpy[3] == v1[missing])
1896  return 0;
1897  }
1898 
1899  return -1;
1900  }
1901 }
1902 
1903 // Only used to draw facets of the medial axis
1904 static bool sortClockwise(SPoint3 a, SPoint3 b, SPoint3 center, SVector3 normal)
1905 {
1906  // If dot(n, cross(A-C, B-C)) is positive, B is counterclockwise from A; if
1907  // it's negative, B is clockwise from A.
1908  SVector3 tmp = crossprod(SVector3(center, a), SVector3(center, b));
1909  return dot(normal, tmp) <= 0;
1910 }
1911 
1912 // Computes the medial axis of the geometry and assigns the feature size at
1913 // vertices
1914 HXTStatus featureSize(Forest *forest)
1915 {
1916  HXTMesh *mesh = forest->forestOptions->mesh;
1917  int nLayersPerGap = forest->forestOptions->nodePerGap;
1918  double hmin = forest->forestOptions->hmin;
1919  double hmax = forest->forestOptions->hmax;
1920 
1921  std::vector<MVertex *> allVertices;
1922  std::vector<double> sizeAtVertices(mesh->vertices.num, DBL_MAX);
1923  for(size_t i = 0; i < mesh->vertices.num; ++i) {
1924  allVertices.push_back(new MVertex(mesh->vertices.coord[4 * i + 0],
1925  mesh->vertices.coord[4 * i + 1],
1926  mesh->vertices.coord[4 * i + 2]));
1927  }
1928 
1929  int firstVertex = allVertices[0]->getNum();
1930 
1931  // All tets
1932  uint64_t count = 0;
1933  std::vector<MTetrahedron *> allTets;
1934  std::vector<std::vector<uint64_t> > tetIncidents;
1935  std::vector<std::vector<MEdge> > edgIncidents;
1936  for(uint32_t i = 0; i < mesh->vertices.num; ++i) {
1937  std::vector<uint64_t> vTet;
1938  tetIncidents.push_back(vTet);
1939  std::vector<MEdge> vEdg;
1940  edgIncidents.push_back(vEdg);
1941  }
1942 
1943  for(size_t i = 0; i < mesh->tetrahedra.num; ++i) {
1944  if(mesh->tetrahedra.node[4 * i + 3] != HXT_GHOST_VERTEX) {
1945  allTets.push_back(
1946  new MTetrahedron(allVertices[mesh->tetrahedra.node[4 * i + 0]],
1947  allVertices[mesh->tetrahedra.node[4 * i + 1]],
1948  allVertices[mesh->tetrahedra.node[4 * i + 2]],
1949  allVertices[mesh->tetrahedra.node[4 * i + 3]]));
1950 
1951  for(size_t j = 0; j < 4; ++j) {
1952  tetIncidents[mesh->tetrahedra.node[4 * i + j]].push_back(count);
1953  }
1954  for(size_t j = 0; j < 6; ++j) {
1955  MEdge e = allTets[count]->getEdge(j);
1956  edgIncidents[allTets[count]->getEdge(j).getVertex(0)->getNum() -
1957  firstVertex]
1958  .push_back(e);
1959  edgIncidents[allTets[count]->getEdge(j).getVertex(1)->getNum() -
1960  firstVertex]
1961  .push_back(e);
1962  }
1963  ++count;
1964  }
1965  }
1966 
1967  std::set<MEdge, MEdgeLessThan> axis;
1968  int elemDrawn = 0;
1969 
1970  FILE *file = fopen("medialAxis_toDraw.pos", "w");
1971  if(file == nullptr) return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
1972  FILE *file2 = fopen("keptEdges.pos", "w");
1973  if(file2 == nullptr) return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
1974 
1975  bool draw = true;
1976  if(draw) {
1977  fprintf(file, "View \"medialAxis\" {\n");
1978  fprintf(file2, "View \"keptEdges\" {\n");
1979  }
1980 
1981  int indFace;
1982 
1983  for(size_t i = 0; i < mesh->vertices.num; ++i) {
1984  // Pôle de p : circumcenter le plus loin parmi les tets incidents à p
1985  SPoint3 pole(0., 0., 0.), tmp(0., 0., 0.),
1986  p(mesh->vertices.coord[4 * i + 0], mesh->vertices.coord[4 * i + 1],
1987  mesh->vertices.coord[4 * i + 2]);
1988  double d = 0.;
1989 
1990  for(size_t j = 0; j < tetIncidents[i].size(); ++j) {
1991  tmp = allTets[tetIncidents[i][j]]->circumcenter();
1992  if(p.distance(tmp) > d) pole = tmp;
1993  d = fmax(d, p.distance(tmp));
1994  }
1995 
1996  // Pole vector
1997  SPoint3 vp = pole - p;
1998  double D = -(vp[0] * p[0] + vp[1] * p[1] + vp[2] * p[2]);
1999  SPoint3 p1(
2000  0., 0.,
2001  -D / vp[2]); // 2 points sur le plan qui passe par p et de normale vp
2002  SPoint3 p2(0., -D / vp[1], 0.);
2003 
2004  std::vector<MFace> up; // umbrella
2005  // Boucle sur les voronoi edges (paires de centres)
2006  double orientj, orientk;
2007  for(size_t j = 0; j < tetIncidents[i].size(); ++j) {
2008  uint64_t tetj = tetIncidents[i][j];
2009  SPoint3 cj = allTets[tetj]->circumcenter();
2010  for(size_t k = j; k < tetIncidents[i].size(); ++k) {
2011  uint64_t tetk = tetIncidents[i][k];
2012  if(tetj != tetk) {
2013  indFace = commonFaceTetFast(allTets[tetj], allTets[tetk]);
2014  if(indFace >= 0) {
2015  SPoint3 ck = allTets[tetk]->circumcenter();
2016  orientj = robustPredicates::orient3d((double *)p, (double *)p1,
2017  (double *)p2, (double *)cj);
2018  orientk = robustPredicates::orient3d((double *)p, (double *)p1,
2019  (double *)p2, (double *)ck);
2020  if(orientj * orientk < 0) {
2021  up.push_back(allTets[tetj]->getFace(indFace));
2022  // break; ?
2023  }
2024  }
2025  }
2026  }
2027  }
2028 
2029  double theta = M_PI / 8., rho = 8., maxAngle, minRatio, localAngle, alpha0,
2030  alpha1;
2031  std::vector<MEdge> checkedEdges;
2032  bool checked;
2033  for(size_t j = 0; j < edgIncidents[i].size();
2034  ++j) { // boucler sur les aretes incidentes à p
2035 
2036  MEdge e = edgIncidents[i][j];
2037  checked = false;
2038  for(size_t k = 0; k < checkedEdges.size(); ++k) {
2039  if(e == checkedEdges[k]) {
2040  checked = true;
2041  break;
2042  }
2043  }
2044 
2045  if(checked) { continue; }
2046  else {
2047  checkedEdges.push_back(e);
2048  }
2049 
2050  maxAngle = 0.0;
2051  minRatio = DBL_MAX;
2052 
2053  uint32_t v0 = e.getVertex(0)->getNum() - firstVertex;
2054  uint32_t v1 = e.getVertex(1)->getNum() - firstVertex;
2055  if(v0 == i || v1 == i) {
2056  for(size_t l = 0; l < up.size(); ++l) {
2057  // Angle condition
2058  localAngle = angle(e.tangent(), up[l].normal());
2059  localAngle = fmin(localAngle, fabs(M_PI - localAngle));
2060  maxAngle = fmax(maxAngle, localAngle);
2061 
2062  // Ratio condition
2063  MTriangle tri(up[l].getVertex(0), up[l].getVertex(1),
2064  up[l].getVertex(2));
2065  minRatio = fmin(minRatio, e.length() / tri.getOuterRadius());
2066  }
2067 
2068  if(maxAngle < M_PI / 2. - theta || minRatio > rho) {
2069  double *n0 = forest->forestOptions->nodeNormals + 3 * v0;
2070  double *n1 = forest->forestOptions->nodeNormals + 3 * v1;
2071  alpha0 = angle(SVector3(n0), e.tangent());
2072  alpha1 = angle(SVector3(n1), e.tangent());
2073 
2074  if(fmin(alpha0, fabs(M_PI - alpha0)) < M_PI / 8. &&
2075  fmin(alpha1, fabs(M_PI - alpha1)) < M_PI / 8.) {
2076  // Add edge to the set (axis, though actually unused), modifiy size
2077  // at its extrmities and draw the dual facet
2078  auto ret = axis.insert(e);
2079 
2080  if(ret.second) {
2081  double h = e.length() / nLayersPerGap;
2082  h = fmax(h, hmin);
2083  h = fmin(h, hmax);
2084  sizeAtVertices[v0] = fmin(h, sizeAtVertices[v0]);
2085  sizeAtVertices[v1] = fmin(h, sizeAtVertices[v1]);
2086 
2087  fprintf(file2, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n",
2088  e.getVertex(0)->x(), e.getVertex(0)->y(),
2089  e.getVertex(0)->z(), e.getVertex(1)->x(),
2090  e.getVertex(1)->y(), e.getVertex(1)->z(), e.length(),
2091  e.length());
2092 
2093  if(draw) {
2094  // Turning around the edge to draw the facet
2095  std::vector<SPoint3> centers;
2096  for(size_t jj = 0; jj < tetIncidents[i].size(); ++jj) {
2097  uint64_t tetj = tetIncidents[i][jj];
2098  for(size_t b = 0; b < 6; ++b) {
2099  if(allTets[tetj]->getEdge(b) == e)
2100  centers.push_back(allTets[tetj]->circumcenter());
2101  }
2102  }
2103  if(centers.size() > 2) {
2104  // Center of the face
2105  SPoint3 c(0., 0., 0.);
2106  for(size_t a = 0; a < centers.size(); ++a) {
2107  c += centers[a];
2108  }
2109  c /= centers.size();
2110  SVector3 normal =
2111  crossprod(SVector3(c, centers[0]), SVector3(c, centers[1]));
2112  normal.normalize();
2113  // Sort clockwise around center
2114  sort(centers.begin(), centers.end(),
2115  [c, normal](SPoint3 a, SPoint3 b) {
2116  return sortClockwise(a, b, c, normal);
2117  });
2118  // Draw
2119  for(size_t a = 1; a < centers.size() - 1; ++a) {
2120  fprintf(file,
2121  "ST(%f,%f,%f, %f,%f,%f, %f,%f,%f){%d,%d,%d};\n",
2122  centers[0][0], centers[0][1], centers[0][2],
2123  centers[a][0], centers[a][1], centers[a][2],
2124  centers[a + 1][0], centers[a + 1][1],
2125  centers[a + 1][2], 1, 1, 1);
2126  }
2127  fprintf(file, "ST(%f,%f,%f, %f,%f,%f, %f,%f,%f){%d,%d,%d};\n",
2128  centers[0][0], centers[0][1], centers[0][2],
2129  centers[centers.size() - 1][0],
2130  centers[centers.size() - 1][1],
2131  centers[centers.size() - 1][2], centers[1][0],
2132  centers[1][1], centers[1][2], 1, 1, 1);
2133  ++elemDrawn;
2134  }
2135  }
2136  } // if edge was inserted
2137  } // if edges does not have a too large angle with normals to its
2138  // extremities
2139  } // if edge passes conditions
2140  } // if i is an edge vertex
2141  } // for incident edges
2142  } // for vertices.num
2143 
2144  if(draw) {
2145  fprintf(file, "};");
2146  fclose(file);
2147  }
2148 
2149  // // Assign feature size in the octree cells
2150  // sc_array_t *points = sc_array_new_size(sizeof(size_point_t),
2151  // mesh->vertices.num); size_point_t *p_tmp; for(size_t i = 0; i <
2152  // mesh->vertices.num; ++i){
2153  // p_tmp = (size_point_t *) sc_array_index(points, i);
2154  // p_tmp->x = mesh->vertices.coord[(size_t) 4*i+0];
2155  // p_tmp->y = mesh->vertices.coord[(size_t) 4*i+1];
2156  // p_tmp->z = mesh->vertices.coord[(size_t) 4*i+2];
2157  // p_tmp->size = sizeAtVertices[i];
2158  // }
2159  // p4est_search(forest->p4est, NULL, replace, points);
2160 
2161  for(size_t i = 0; i < mesh->vertices.num; ++i) {
2162  (*forest->forestOptions->featureSizeAtVertices)[i] = sizeAtVertices[i];
2163  }
2164 
2165  return HXT_STATUS_OK;
2166 }
2167 
2168 /* ========================================================================================================
2169  ESTIMATE NUMBER OF TETRAHEDRA IN THE VOLUME MESH
2170  ========================================================================================================
2171  */
2172 
2173 static void elementEstimate(p4est_iter_volume_info_t *info, void *user_data)
2174 {
2175  p4est_quadrant_t *q = info->quad;
2176  size_data_t *data = (size_data_t *)q->p.user_data;
2177  p4est_t *p4est = info->p4est;
2178  p4est_topidx_t which_tree = info->treeid;
2179 
2180  double center[3];
2181  getCellCenter(p4est, which_tree, q, center);
2182 
2183  double octantVolume = data->h * data->h * data->h;
2184  double tetVolume = data->size * data->size * data->size * sqrt(2) / 12.0;
2185 
2186  *((double *)user_data) += octantVolume / tetVolume;
2187 }
2188 
2189 HXTStatus hxtOctreeElementEstimation(Forest *forest, double *elemEstimate)
2190 {
2191  p4est_iterate(forest->p4est, nullptr, (void *)elemEstimate, elementEstimate,
2192  nullptr, nullptr, nullptr);
2193  return HXT_STATUS_OK;
2194 }
2195 
2196 /* ========================================================================================================
2197  INTERPOLATION OF CURVATURE DIRECTIONS
2198  ========================================================================================================
2199  */
2200 
2201 int intersections;
2202 
2203 static void markIntersectingCells(p4est_iter_volume_info_t *info,
2204  void *user_data)
2205 {
2206  ForestOptions *forestOptions = (ForestOptions *)info->p4est->user_pointer;
2207  size_data_t *data = (size_data_t *)info->quad->p.user_data;
2208  if(forestOptions->dim == 3) {
2209  double min[3], max[3];
2210  getCellBBox(info->p4est, info->treeid, info->quad, min, max);
2211  std::vector<uint64_t> candidates;
2212  forestOptions->triRTree->Search(min, max, rtreeCallback, &candidates);
2213  if(!candidates.empty()) {
2214  data->hasIntersection = true;
2215  ++intersections;
2216  }
2217  }
2218 }
2219 
2220 static void pushInterpolationData(p4est_iter_volume_info_t *info,
2221  void *user_data)
2222 {
2223  size_data_t *data = (size_data_t *)info->quad->p.user_data;
2224  if(!data->hasIntersection) {
2225  p4est_quadrant_t *q = info->quad;
2226  p4est_t *p4est = info->p4est;
2227  p4est_topidx_t which_tree = info->treeid;
2228 
2229  std::vector<interpolation_data_t> *cellCenters =
2230  (std::vector<interpolation_data_t> *)user_data;
2231 
2232  double center[3];
2233  getCellCenter(p4est, which_tree, q, center);
2234  interpolation_data_t intdata = {
2235  {center[0], center[1], center[2]}, SVector3(0.), SVector3(0.), 0.};
2236  (*cellCenters).push_back(intdata);
2237  }
2238 }
2239 
2240 int cellCounter;
2241 
2242 static void addDistanceContribution(p4est_iter_volume_info_t *info,
2243  void *user_data)
2244 {
2245  size_data_t *data = (size_data_t *)info->quad->p.user_data;
2246  if(data->hasIntersection) {
2247  p4est_quadrant_t *q = info->quad;
2248  p4est_t *p4est = info->p4est;
2249  p4est_topidx_t which_tree = info->treeid;
2250  std::vector<interpolation_data_t> *cellCenters =
2251  (std::vector<interpolation_data_t> *)user_data;
2252  double center[3];
2253  getCellCenter(p4est, which_tree, q, center);
2254  double dist, xc, yc, zc;
2255  printf("working at cell %d/%d\n", cellCounter++, intersections);
2256 
2257  for(size_t i = 0; i < (*cellCenters).size(); ++i) {
2258  xc = (*cellCenters)[i].center[0];
2259  yc = (*cellCenters)[i].center[1];
2260  zc = (*cellCenters)[i].center[2];
2261  dist = sqrt((center[0] - xc) * (center[0] - xc) +
2262  (center[1] - yc) * (center[1] - yc) +
2263  (center[2] - zc) * (center[2] - zc));
2264  (*cellCenters)[i].denom += 1. / dist;
2265  (*cellCenters)[i].t1 += data->t1 * (1. / dist);
2266  (*cellCenters)[i].t2 += data->t2 * (1. / dist);
2267  }
2268  }
2269 }
2270 
2271 static void assignDirections(p4est_iter_volume_info_t *info, void *user_data)
2272 {
2273  size_data_t *data = (size_data_t *)info->quad->p.user_data;
2274  if(!data->hasIntersection) {
2275  p4est_quadrant_t *q = info->quad;
2276  p4est_t *p4est = info->p4est;
2277  p4est_topidx_t which_tree = info->treeid;
2278  std::vector<interpolation_data_t> *cellCenters =
2279  (std::vector<interpolation_data_t> *)user_data;
2280  double center[3];
2281  getCellCenter(p4est, which_tree, q, center);
2282 
2283  double xc, yc, zc;
2284  for(size_t i = 0; i < (*cellCenters).size(); ++i) {
2285  xc = (*cellCenters)[i].center[0];
2286  yc = (*cellCenters)[i].center[1];
2287  zc = (*cellCenters)[i].center[2];
2288  if(fabs(center[0] - xc) < 1e-13 && fabs(center[1] - yc) < 1e-13 &&
2289  fabs(center[2] - zc) < 1e-13) {
2290  data->t1 = (*cellCenters)[i].t1 * (1. / (*cellCenters)[i].denom);
2291  data->t2 = (*cellCenters)[i].t2 * (1. / (*cellCenters)[i].denom);
2292  data->n = crossprod(data->t1, data->t2);
2293  data->n = data->n.unit();
2294  break;
2295  }
2296  }
2297  }
2298 }
2299 
2300 static void drawDirections(p4est_iter_volume_info_t *info, void *user_data)
2301 {
2302  size_data_t *data = (size_data_t *)info->quad->p.user_data;
2303  ForestOptions *forestOptions = (ForestOptions *)info->p4est->user_pointer;
2304  p4est_quadrant_t *q = info->quad;
2305  p4est_t *p4est = info->p4est;
2306  p4est_topidx_t which_tree = info->treeid;
2307 
2308  double center[3];
2309  getCellCenter(p4est, which_tree, q, center);
2310 
2311  double x = center[0];
2312  double y = center[1];
2313  double z = center[2];
2314  fprintf(forestOptions->file3, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n", x, y, z,
2315  x + data->t1[0], y + data->t1[1], z + data->t1[2], 1.0, 1.0);
2316  fprintf(forestOptions->file3, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n", x, y, z,
2317  x + data->t2[0], y + data->t2[1], z + data->t2[2], 1.0, 1.0);
2318  fprintf(forestOptions->file3, "SL(%f,%f,%f, %f,%f,%f){%f,%f};\n", x, y, z,
2319  x + data->n[0], y + data->n[1], z + data->n[2], 1.0, 1.0);
2320 }
2321 
2322 HXTStatus forestInterpolateDirections(Forest *forest)
2323 {
2324  intersections = 0;
2325  p4est_iterate(forest->p4est, nullptr, nullptr, markIntersectingCells, nullptr,
2326  nullptr, nullptr);
2327 
2328  std::vector<interpolation_data_t> cellCenters;
2329  cellCounter = 0;
2330  p4est_iterate(forest->p4est, nullptr, &cellCenters, pushInterpolationData,
2331  nullptr, nullptr, nullptr);
2332  p4est_iterate(forest->p4est, nullptr, &cellCenters, addDistanceContribution,
2333  nullptr, nullptr, nullptr);
2334  p4est_iterate(forest->p4est, nullptr, &cellCenters, assignDirections, nullptr,
2335  nullptr, nullptr);
2336 
2337  forest->forestOptions->file3 = fopen("directions.pos", "w");
2338  if(forest->forestOptions->file3 == nullptr)
2339  return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
2340 
2341  fprintf(forest->forestOptions->file3, "View \"directions\" {\n");
2342 
2343  p4est_iterate(forest->p4est, nullptr, nullptr, drawDirections, nullptr,
2344  nullptr, nullptr);
2345 
2346  fprintf(forest->forestOptions->file3, "};");
2347  fclose(forest->forestOptions->file3);
2348 
2349  return HXT_STATUS_OK;
2350 }
2351 
2352 /* ========================================================================================================
2353  EXPORT
2354  ========================================================================================================
2355  */
2356 
2357 HXTStatus saveGlobalData(Forest *forest, const char *filename)
2358 {
2359  FILE *f = fopen(filename, "w");
2360  if(f == nullptr) return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
2361  fprintf(f, "%f %f\n", forest->forestOptions->hmin,
2362  forest->forestOptions->hmax);
2363  fclose(f);
2364  Msg::Info("Saved global size field data to %s", filename);
2365  return HXT_STATUS_OK;
2366 }
2367 
2368 static void exportToHexCallback(p4est_iter_volume_info_t *info, void *user_data)
2369 {
2370  p4est_quadrant_t *q = info->quad;
2371  size_data_t *data = (size_data_t *)q->p.user_data;
2372  p4est_t *p4est = info->p4est;
2373  p4est_topidx_t which_tree = info->treeid;
2374 
2375  FILE *f = (FILE *)user_data;
2376  double center[3], x[8], y[8], z[8];
2377  getCellCenter(p4est, which_tree, q, center);
2378  double h = data->h / 2.0, s = data->size, epsilon = 1e-12;
2379  x[0] = x[3] = x[4] = x[7] = center[0] - h - epsilon;
2380  x[1] = x[2] = x[5] = x[6] = center[0] + h + epsilon;
2381  y[0] = y[1] = y[4] = y[5] = center[1] - h - epsilon;
2382  y[2] = y[3] = y[6] = y[7] = center[1] + h + epsilon;
2383  z[0] = z[1] = z[2] = z[3] = center[2] - h - epsilon;
2384  z[4] = z[5] = z[6] = z[7] = center[2] + h + epsilon;
2385 
2386  fprintf(f,
2387  "SH(%f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f,%f,%f,%f, %f,%f,%f, "
2388  "%f,%f,%f, %f,%f,%f){%f,%f,%f,%f,%f,%f,%f,%f};\n",
2389  x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], x[3], y[3],
2390  z[3], x[4], y[4], z[4], x[5], y[5], z[5], x[6], y[6], z[6], x[7],
2391  y[7], z[7], s, s, s, s, s, s, s, s);
2392  // fprintf(f, "SH(%f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f,%f,%f,%f, %f,%f,%f,
2393  // %f,%f,%f, %f,%f,%f){%d, %d, %d, %d, %d, %d, %d, %d};\n",
2394  // x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], x[3], y[3], z[3],
2395  // x[4], y[4], z[4], x[5], y[5], z[5], x[6], y[6], z[6], x[7], y[7], z[7],
2396  // data->hasIntersection, data->hasIntersection, data->hasIntersection,
2397  // data->hasIntersection, data->hasIntersection, data->hasIntersection,
2398  // data->hasIntersection, data->hasIntersection);
2399 }
2400 
2401 static void exportToQuadCallback(p4est_iter_volume_info_t *info,
2402  void *user_data)
2403 {
2404  p4est_quadrant_t *q = info->quad;
2405  size_data_t *data = (size_data_t *)q->p.user_data;
2406 
2407  p4est_t *p4est = info->p4est;
2408  p4est_topidx_t which_tree = info->treeid;
2409 
2410  FILE *f = (FILE *)user_data;
2411 
2412  double center[3], x[8], y[8], z[8];
2413  getCellCenter(p4est, which_tree, q, center);
2414 
2415  double h = data->h / 2.0, s = data->size, epsilon = 1e-12;
2416  x[0] = x[3] = center[0] - h - epsilon;
2417  x[1] = x[2] = center[0] + h + epsilon;
2418  y[0] = y[1] = center[1] - h - epsilon;
2419  y[2] = y[3] = center[1] + h + epsilon;
2420  z[0] = z[1] = 0.0;
2421  z[2] = z[3] = 0.0;
2422 
2423  fprintf(f, "SQ(%f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f){%f,%f,%f,%f};\n", x[0],
2424  y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], x[3], y[3], z[3], s,
2425  s, s, s);
2426 }
2427 
2428 HXTStatus forestExport(Forest *forest, const char *forestFile)
2429 {
2430  FILE *f = fopen(forestFile, "w");
2431  if(f == nullptr) return HXT_ERROR(HXT_STATUS_FILE_CANNOT_BE_OPENED);
2432 
2433  fprintf(f, "View \"sizeField\" {\n");
2434  if(forest->forestOptions->dim == 3) {
2435  p4est_iterate(forest->p4est, nullptr, (void *)f, exportToHexCallback,
2436  nullptr, nullptr, nullptr);
2437  }
2438  else if(forest->forestOptions->dim == 2) {
2439  p4est_iterate(forest->p4est, nullptr, (void *)f, exportToQuadCallback,
2440  nullptr, nullptr, nullptr);
2441  }
2442  fprintf(f, "};");
2443  fclose(f);
2444  return HXT_STATUS_OK;
2445 }
2446 
2447 HXTStatus forestSave(Forest *forest, const char *forestFile,
2448  const char *dataFile)
2449 {
2450  HXT_CHECK(saveGlobalData(forest, dataFile));
2451  p4est_save_ext(forestFile, forest->p4est, true, false);
2452  return HXT_STATUS_OK;
2453 }
2454 
2455 #endif
2456 
2457 /* ======================================================================================
2458  End functions from hxt_octree
2459  ======================================================================================
2460  */
2461 
2462 double automaticMeshSizeField::operator()(double X, double Y, double Z,
2463  GEntity *ge)
2464 {
2465  double val = 1.e17;
2466 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
2467  HXTStatus s = forestSearchOne(forest, X, Y, Z, &val, true);
2468  if(s == HXT_STATUS_OK) { return val; }
2469  else
2470  Msg::Error("Cannot find point %g %g %g in the octree", X, Y, Z);
2471 #else
2472  Msg::Error("Gmsh has to be compiled with HXT and P4EST for using "
2473  "automaticMeshSizeField");
2474 #endif
2475  return val;
2476 }
2477 
2478 void automaticMeshSizeField::operator()(double X, double Y, double Z,
2479  SMetric3 &m, GEntity *ge)
2480 {
2481 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
2482  HXTStatus s = forestSearchOneAniso(forest, X, Y, Z, m, false);
2483  if(fabs(m.determinant()) < 1e-13) m = SMetric3();
2484  m.print("m");
2485  if(s != HXT_STATUS_OK)
2486  Msg::Error("Cannot find point %g %g %g in the octree", X, Y, Z);
2487 #else
2488  Msg::Error("Gmsh has to be compiled with HXT and P4EST for using "
2489  "automaticMeshSizeField");
2490 #endif
2491 }
2492 
2493 automaticMeshSizeField::~automaticMeshSizeField()
2494 {
2495 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
2496  if(forest) forestDelete(&forest);
2497  if(forestOptions) forestOptionsDelete(&forestOptions);
2498 #endif
2499 }
2500 
2501 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
2502 HXTStatus automaticMeshSizeField::updateHXT()
2503 {
2504  if(!updateNeeded) return HXT_STATUS_OK;
2505 
2506  if(forestOptions) HXT_CHECK(forestOptionsDelete(&forestOptions));
2507  if(forest) HXT_CHECK(forestDelete(&forest));
2508 
2509  updateNeeded = false;
2510 
2511  if(!_forestFile.empty()) {
2512  // Load .p4est file if given a valid file name
2513  Msg::Info("Loading size field from %s", _forestFile.c_str());
2514  HXT_CHECK(forestOptionsCreate(&forestOptions));
2515  size_t lastindex = _forestFile.find_last_of(".");
2516  std::string root = _forestFile.substr(0, lastindex);
2517  std::string forestFile = root + ".p4est";
2518  std::string dataFile = root + ".data";
2519  HXT_CHECK(
2520  forestLoad(&forest, forestFile.c_str(), dataFile.c_str(), forestOptions));
2521  }
2522  else {
2523  // Compute the size field otherwise
2524  int dim = GModel::current()->getDim();
2525 
2526  HXT_CHECK(forestOptionsCreate(&forestOptions));
2527 
2528  Msg::Info("Gradation = %f\n", _gradation);
2529  Msg::Info("Node density = %d\n", _nPointsPerCircle);
2530  if(dim == 3) {
2531  if(_nPointsPerGap > 0) {
2532  Msg::Info("Layers per gap = %d\n", _nPointsPerGap);
2533  }
2534  else {
2535  Msg::Info("Layers per gap = %d : not detecting features.\n",
2536  _nPointsPerGap);
2537  }
2538  }
2539 
2540  // The bounding box of the mesh/model
2541  double bbox_vertices[6];
2542 
2543  RTree<uint64_t, double, 3> triRTree;
2544  HXTMesh *mesh;
2545  double *nodalCurvature;
2546  // double *nodeNormals;
2547  // std::vector<std::function<double(double)>> curvFunctions;
2548  // std::vector<std::function<double(double)>> xFunctions;
2549  // std::vector<std::function<double(double)>> yFunctions;
2550 
2551  int debug = true;
2552 
2553  GModel *gm = GModel::current();
2554 
2555  size_t numVertices = gm->getNumMeshVertices();
2556  HXT_CHECK(hxtMalloc(&nodalCurvature, 6 * numVertices * sizeof(double)));
2557  std::vector<double> nodeNormals;
2558  for(size_t i = 0; i < 6 * numVertices; ++i) { nodalCurvature[i] = NAN; }
2559 
2560  if(dim == 3) {
2561  // Get all faces of the model
2562  std::vector<GFace *> faces;
2563  std::vector<GRegion *> regions;
2564  for(auto it = GModel::current()->firstRegion();
2565  it != GModel::current()->lastRegion(); ++it) {
2566  regions.push_back(*it);
2567  }
2568  if(!regions.empty()) {
2569  HXT_CHECK(getAllFacesOfAllRegions(regions, nullptr, faces));
2570  }
2571  else {
2572  Msg::Info("No volume in the model : looping over faces instead.");
2573  for(auto it = GModel::current()->firstFace();
2574  it != GModel::current()->lastFace(); ++it) {
2575  faces.push_back(*it);
2576  }
2577  }
2578 
2579  if(regions.empty()) {
2580  Msg::Error("Erreur : Pas de volume dans le modèle.");
2581  }
2582 
2583  // Create global HXT mesh structure
2584  HXT_CHECK(hxtMeshCreate(&mesh));
2585  std::map<MVertex *, uint32_t> v2c;
2586  std::vector<MVertex *> c2v;
2587  Gmsh2Hxt(faces, mesh, v2c, c2v);
2588  // Gmsh2Hxt(regions, mesh, v2c, c2v);
2589 
2590  // if(!regions.empty()){
2591  // Msg::Info("Volume found.");
2592  // HXT_CHECK( getAllFacesOfAllRegions(regions, NULL, faces) );
2593  // } else{
2594  // Msg::Info("No volume in the model : looping over faces instead.");
2595  // regions[0] = new GRegion(GModel::current(),-1);
2596  // for(auto it = GModel::current()->firstFace(); it !=
2597  // GModel::current()->lastFace(); ++it){
2598  // faces.push_back(*it);
2599  // regions[0]->addEmbeddedFace(*it);
2600  // }
2601  // }
2602 
2603  // // Create global HXT mesh structure
2604  // HXT_CHECK(hxtMeshCreate(&mesh));
2605  // std::map<MVertex *, uint32_t> v2c;
2606  // std::vector<MVertex *> c2v;
2607  // // Gmsh2Hxt(faces, mesh, v2c, c2v);
2608  // Gmsh2Hxt(regions, mesh, v2c, c2v);
2609 
2610  if(mesh->vertices.num == 0) {
2611  Msg::Error("Surface mesh is empty");
2612  HXT_CHECK(hxtMeshDelete(&mesh));
2613  Msg::Exit(1);
2614  }
2615 
2616  // HXT_CHECK(hxtMalloc(&nodalCurvature,6*mesh->vertices.num*sizeof(double)));
2617  // for(uint64_t i = 0; i < 6*mesh->vertices.num; ++i){ nodalCurvature[i] =
2618  // NAN; printf("%d\n",i);}
2619 
2620  // Create HXT mesh structure for each GFace
2621  std::vector<HXTMesh *> faceMeshes;
2622  for(size_t i = 0; i < faces.size(); ++i) {
2623  HXTMesh *meshFace;
2624  HXT_CHECK(hxtMeshCreate(&meshFace));
2625  std::vector<GFace *> oneFace;
2626  oneFace.push_back(faces[i]);
2627  std::map<MVertex *, uint32_t> v2cLoc;
2628  std::vector<MVertex *> c2vLoc;
2629  Gmsh2HxtLocal(oneFace, meshFace, v2cLoc, c2vLoc);
2630  faceMeshes.push_back(meshFace);
2631  }
2632 
2633  std::map<MVertex *, uint32_t> v2c2;
2634  std::vector<MVertex *> c2v2;
2635  Gmsh2HxtGlobal(faces, nullptr, v2c2, c2v2);
2636 
2637  size_t nVertices = 0;
2638 
2639  assert(faces.size() == faceMeshes.size());
2640 
2641  // Compute curvature of the faces
2642  int counter = 0;
2643  for(auto it = gm->firstFace(); it != gm->lastFace(); ++it) {
2644  HXTMesh *meshFace = faceMeshes[counter++];
2645 
2646  if(meshFace == nullptr) { Msg::Error("meshFace == NULL"); }
2647 
2648  GFace *gf = *it;
2649  std::map<MVertex *, int> nodeIndex;
2650  std::vector<SPoint3> nodes;
2651  std::vector<int> tris;
2652  std::vector<std::pair<SVector3, SVector3> > curv;
2653  for(std::size_t i = 0; i < gf->triangles.size(); i++) {
2654  MTriangle *t = gf->triangles[i];
2655  for(int j = 0; j < 3; j++) {
2656  MVertex *v = t->getVertex(j);
2657  if(nodeIndex.find(v) == nodeIndex.end()) {
2658  int idx = nodes.size();
2659  nodeIndex[v] = idx;
2660  nodes.push_back(v->point());
2661  tris.push_back(idx);
2662  }
2663  else {
2664  tris.push_back(nodeIndex[v]);
2665  }
2666  }
2667  }
2668 
2669  for(size_t i = 0; i < meshFace->vertices.num; ++i) {
2670  nodes[i] = SPoint3(meshFace->vertices.coord[(size_t)4 * i + 0],
2671  meshFace->vertices.coord[(size_t)4 * i + 1],
2672  meshFace->vertices.coord[(size_t)4 * i + 2]);
2673  // nodes[i] = gm->getMeshVertexByTag(i)
2674  }
2675 
2676  for(size_t i = 0; i < meshFace->triangles.num; ++i) {
2677  tris[3 * i + 0] = meshFace->triangles.node[3 * i + 0];
2678  tris[3 * i + 1] = meshFace->triangles.node[3 * i + 1];
2679  tris[3 * i + 2] = meshFace->triangles.node[3 * i + 2];
2680  }
2681 
2682  if(gf->triangles.empty()) {
2683  Msg::Info("Skipping curvature computation on face %d with 0 element",
2684  counter - 1);
2685  }
2686  else {
2687  // Compute curvature of the face
2688  CurvatureRusinkiewicz(tris, nodes, curv);
2689 
2690  // Assemble curvature vectors of the face in global nodalCurvature
2691  // structure
2692  for(uint64_t i = 0; i < meshFace->vertices.num; ++i) {
2693  uint64_t nodeGlobal = v2c[c2v2[nVertices + i]];
2694 
2695  debug = true;
2696  if(debug) { // Check the vertex of the face and the global vertex
2697  // are the same
2698  double x1, y1, z1, x2, y2, z2;
2699  x1 = meshFace->vertices.coord[(size_t)4 * i + 0];
2700  y1 = meshFace->vertices.coord[(size_t)4 * i + 1];
2701  z1 = meshFace->vertices.coord[(size_t)4 * i + 2];
2702  x2 = mesh->vertices.coord[(size_t)4 * nodeGlobal + 0];
2703  y2 = mesh->vertices.coord[(size_t)4 * nodeGlobal + 1];
2704  z2 = mesh->vertices.coord[(size_t)4 * nodeGlobal + 2];
2705  if(!isPoint(x1, y1, z1, x2, y2, z2, 1e-12)) {
2706  printf("Mismatch : (%10.12e,%10.12e,%10.12e) - "
2707  "(%10.12e,%10.12e,%10.12e)\n",
2708  x1, y1, z1, x2, y2, z2);
2709  }
2710  assert(isPoint(x1, y1, z1, x2, y2, z2, 1e-15));
2711  }
2712 
2713  nodalCurvature[6 * nodeGlobal + 0] = curv[i].first[0];
2714  nodalCurvature[6 * nodeGlobal + 1] = curv[i].first[1];
2715  nodalCurvature[6 * nodeGlobal + 2] = curv[i].first[2];
2716  nodalCurvature[6 * nodeGlobal + 3] = curv[i].second[0];
2717  nodalCurvature[6 * nodeGlobal + 4] = curv[i].second[1];
2718  nodalCurvature[6 * nodeGlobal + 5] = curv[i].second[2];
2719  }
2720 
2721  nVertices += meshFace->vertices.num;
2722  }
2723  } // for faces.size()
2724 
2725  debug = true;
2726  if(debug)
2727  writeNodalCurvature(nodalCurvature, mesh->vertices.num,
2728  "nodalCurvature.txt");
2729 
2730  // Compute Delaunay tetrahedrization of the (empty) surface mesh
2731  HXTDelaunayOptions delaunayOptions = {nullptr, nullptr, 0, 0, 0, 2, 1, 0};
2732  HXT_CHECK(hxtEmptyMesh(mesh, &delaunayOptions));
2733 
2734  // Compute normal vectors
2735  std::vector<int> tris(3 * mesh->triangles.num, 0);
2736  for(std::size_t i = 0; i < mesh->triangles.num; ++i) {
2737  tris[3 * i + 0] = mesh->triangles.node[3 * i + 0];
2738  tris[3 * i + 1] = mesh->triangles.node[3 * i + 1];
2739  tris[3 * i + 2] = mesh->triangles.node[3 * i + 2];
2740  }
2741 
2742  std::vector<SPoint3> nodes(mesh->vertices.num);
2743  for(size_t i = 0; i < mesh->vertices.num; ++i) {
2744  nodes[i] = SPoint3(mesh->vertices.coord[(size_t)4 * i + 0],
2745  mesh->vertices.coord[(size_t)4 * i + 1],
2746  mesh->vertices.coord[(size_t)4 * i + 2]);
2747  }
2748 
2749  std::vector<std::pair<SVector3, SVector3> > curv;
2750  nodeNormals.reserve(3 * mesh->vertices.num);
2751 
2752  // Same function but returns the normals
2753  CurvatureRusinkiewicz(tris, nodes, curv, nodeNormals);
2754 
2755  // Add bboxes of the surface mesh to rtree
2756  HXTBbox bbox_triangle;
2757  for(uint64_t i = 0; i < mesh->triangles.num; ++i) {
2758  hxtBboxInit(&bbox_triangle);
2759  for(uint64_t j = 0; j < 3; ++j) {
2760  double coord[3];
2761  uint32_t node = mesh->triangles.node[3 * i + j];
2762  for(uint32_t k = 0; k < 3; ++k) {
2763  coord[k] = mesh->vertices.coord[(size_t)4 * node + k];
2764  }
2765  hxtBboxAddOne(&bbox_triangle, coord);
2766  }
2767  SBoundingBox3d cube_bbox(bbox_triangle.min[0], bbox_triangle.min[1],
2768  bbox_triangle.min[2], bbox_triangle.max[0],
2769  bbox_triangle.max[1], bbox_triangle.max[2]);
2770  // cube_bbox.makeCube();
2771  triRTree.Insert((double *)(cube_bbox.min()),
2772  (double *)(cube_bbox.max()), i);
2773  }
2774 
2775  // Compute bbox of the mesh
2776  HXTBbox bbox_mesh;
2777  hxtBboxInit(&bbox_mesh);
2778  hxtBboxAdd(&bbox_mesh, mesh->vertices.coord, mesh->vertices.num);
2779  for(int i = 0; i < 3; ++i) {
2780  bbox_vertices[i] = bbox_mesh.min[i];
2781  bbox_vertices[i + 3] = bbox_mesh.max[i];
2782  }
2783 
2784  // Export RTree in .pos file
2785  bool exportRTree = false;
2786  if(exportRTree) {
2787  RTree<uint64_t, double, 3>::Iterator it;
2788  FILE *f = fopen("rtree.pos", "w");
2789  fprintf(f, "View \"sizeField\" {\n");
2790  int itIndex = 0;
2791  double s = 1.0;
2792  double x[8], y[8], z[8];
2793  // Using custom GetNext2 to display intermediary rectangles
2794  for(triRTree.GetFirst(it); !triRTree.IsNull(it);
2795  triRTree.GetNext2(it)) {
2796  int value = triRTree.GetAt(it);
2797  double boundsMin[3] = {0, 0, 0};
2798  double boundsMax[3] = {0, 0, 0};
2799  it.GetBounds(boundsMin, boundsMax);
2800  std::cout << "it[" << itIndex++ << "] " << value << " = ("
2801  << boundsMin[0] << "," << boundsMin[1] << ","
2802  << boundsMin[2] << "," << boundsMax[0] << ","
2803  << boundsMax[1] << "," << boundsMax[2] << ")\n";
2804  x[0] = x[3] = x[4] = x[7] = boundsMin[0];
2805  x[1] = x[2] = x[5] = x[6] = boundsMax[0];
2806  y[0] = y[1] = y[4] = y[5] = boundsMin[1];
2807  y[2] = y[3] = y[6] = y[7] = boundsMax[1];
2808  z[0] = z[1] = z[2] = z[3] = boundsMin[2];
2809  z[4] = z[5] = z[6] = z[7] = boundsMax[2];
2810  fprintf(f,
2811  "SH(%f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f, "
2812  "%f,%f,%f, %f,%f,%f, %f,%f,%f){%f,%f,%f,%f,%f,%f,%f,%f};\n",
2813  x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], x[3],
2814  y[3], z[3], x[4], y[4], z[4], x[5], y[5], z[5], x[6], y[6],
2815  z[6], x[7], y[7], z[7], s, s, s, s, s, s, s, s);
2816  }
2817  // Adding the root rectangle (the bbox)
2818  x[0] = x[3] = x[4] = x[7] = bbox_mesh.min[0];
2819  x[1] = x[2] = x[5] = x[6] = bbox_mesh.max[0];
2820  y[0] = y[1] = y[4] = y[5] = bbox_mesh.min[1];
2821  y[2] = y[3] = y[6] = y[7] = bbox_mesh.max[1];
2822  z[0] = z[1] = z[2] = z[3] = bbox_mesh.min[2];
2823  z[4] = z[5] = z[6] = z[7] = bbox_mesh.max[2];
2824  fprintf(f,
2825  "SH(%f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f, %f,%f,%f, "
2826  "%f,%f,%f, %f,%f,%f, %f,%f,%f){%f,%f,%f,%f,%f,%f,%f,%f};\n",
2827  x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], x[3],
2828  y[3], z[3], x[4], y[4], z[4], x[5], y[5], z[5], x[6], y[6],
2829  z[6], x[7], y[7], z[7], s, s, s, s, s, s, s, s);
2830  fprintf(f, "};");
2831  fclose(f);
2832  }
2833  }
2834 
2835  if(dim == 2) {
2836  Msg::Error("2D size field is not yet operational");
2837 
2838  SBoundingBox3d bbox = GModel::current()->bounds();
2839  for(int i = 0; i < 3; ++i) {
2840  bbox_vertices[i] = bbox.min()[i];
2841  bbox_vertices[i + 3] = bbox.max()[i];
2842  }
2843 
2844  mesh = nullptr;
2845 
2846  // for(auto it = GModel::current()->firstEdge(); it !=
2847  // GModel::current()->lastEdge(); ++it){
2848  // GEdge *e = *it;
2849  // curvFunctions.push_back([e](double par){ return e->curvature(par);
2850  // }); xFunctions.push_back([e](double par){ return e->point(par).x();
2851  // }); yFunctions.push_back([e](double par){ return e->point(par).y();
2852  // });
2853  // }
2854 
2855  // forestOptions->curvFunctions = &curvFunctions;
2856  // forestOptions->xFunctions = &xFunctions;
2857  // forestOptions->yFunctions = &yFunctions;
2858  }
2859 
2860  // Set the bulk size and the min size from the bbox
2861  if(_hbulk < 0 || _hmin < 0) {
2862  double L = -1.0;
2863  for(int i = 0; i < 3; ++i) {
2864  L = fmax(L, bbox_vertices[i + 3] - bbox_vertices[i]);
2865  }
2866  _hbulk < 0 ? _hbulk = L / 20. : _hbulk;
2867  _hmin < 0 ? _hmin = L / 1000. : _hmin;
2868  Msg::Info("Bulk size is set to %f", _hbulk);
2869  Msg::Info("Min size is set to %f", _hmin);
2870  }
2871 
2872  if(_hmax < 0) _hmax = _hbulk;
2873 
2874  std::vector<double> sizeAtVertices(mesh->vertices.num, DBL_MAX);
2875 
2876  forestOptions->dim = dim;
2877  forestOptions->hmax = _hmax;
2878  forestOptions->hmin = _hmin;
2879  forestOptions->hbulk = _hbulk;
2880  forestOptions->gradation = _gradation;
2881  forestOptions->nodePerTwoPi = _nPointsPerCircle;
2882  forestOptions->nodePerGap = _nPointsPerGap;
2883  forestOptions->bbox = bbox_vertices;
2884  forestOptions->sizeFunction = nullptr;
2885  forestOptions->nodalCurvature = nodalCurvature;
2886  forestOptions->nodeNormals = &nodeNormals[0];
2887  forestOptions->featureSizeAtVertices = &sizeAtVertices;
2888  forestOptions->triRTree = (dim == 3) ? &triRTree : nullptr;
2889  forestOptions->mesh = mesh;
2890 
2891  HXT_CHECK(forestCreate(0, nullptr, &forest, nullptr, forestOptions));
2892 
2893  if(dim == 3) {
2894  if(_nPointsPerGap > 0) {
2895  Msg::Info("Detecting features...");
2896  HXT_CHECK(featureSize(forest));
2897  }
2898 
2899  if(_nPointsPerCircle > 0) {
2900  Msg::Info("Refining octree...");
2901  HXT_CHECK(forestRefine(forest));
2902  }
2903 
2904  if(_smoothing) {
2905  Msg::Info("Smoothing size gradient...");
2906  HXT_CHECK(forestSizeSmoothing(forest));
2907  }
2908 
2909  double elemEstimation;
2910  HXT_CHECK(hxtOctreeElementEstimation(forest, &elemEstimation));
2911  Msg::Info("Estimated number of tetrahedra in the bounding box : %ld",
2912  (uint64_t)ceil(elemEstimation));
2913  }
2914  else {
2915  HXT_CHECK(forestRefine(forest));
2916  if(_smoothing) {
2917  Msg::Info("Smoothing size gradient...");
2918  HXT_CHECK(forestSizeSmoothing(forest));
2919  }
2920  }
2921 
2922  // forestInterpolateDirections(forest);
2923 
2924  // Save forest in .p4est file
2925  std::string forestFile = GModel::current()->getName() + ".p4est";
2926  std::string dataFile = GModel::current()->getName() + ".data";
2927  Msg::Info("Saving size field in %s", forestFile.c_str());
2928  HXT_CHECK(forestSave(forest, forestFile.c_str(), dataFile.c_str()));
2929 
2930  debug = false;
2931  if(debug) {
2932  // Export size field in .pos file
2933  forestFile = GModel::current()->getName() + ".pos";
2934  HXT_CHECK(forestExport(forest, forestFile.c_str()));
2935  }
2936 
2937  if(dim == 3) {
2938  if(nodalCurvature) HXT_CHECK(hxtFree(&nodalCurvature));
2939  // if(nodeNormals) HXT_CHECK(hxtFree(&nodeNormals));
2940  HXT_CHECK(hxtMeshDelete(&mesh));
2941  }
2942  }
2943 
2944  return HXT_STATUS_OK;
2945 }
2946 
2947 #endif
2948 
2949 void automaticMeshSizeField::update()
2950 {
2951 #if defined(HAVE_HXT) && defined(HAVE_P4EST)
2952  HXTStatus s = updateHXT();
2953  if(s != HXT_STATUS_OK)
2954  Msg::Error("Something went wrong when computing the octree");
2955 #else
2956  Msg::Error(
2957  "Gmsh has to be compiled with HXT and P4EST to use automaticMeshSizeField");
2958 #endif
2959 };
crossprod
SVector3 crossprod(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:150
D
#define D
Definition: DefaultOptions.h:24
size_point::z
double z
Definition: automaticMeshSizeField.h:91
automaticMeshSizeField::_forestFile
std::string _forestFile
Definition: automaticMeshSizeField.h:131
MEdge
Definition: MEdge.h:14
ForestOptions::nodePerGap
int nodePerGap
Definition: automaticMeshSizeField.h:46
ForestOptions::gradation
double gradation
Definition: automaticMeshSizeField.h:44
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
ForestOptions::file1
FILE * file1
Definition: automaticMeshSizeField.h:59
SMetric3
Definition: STensor3.h:17
ForestOptions::file3
FILE * file3
Definition: automaticMeshSizeField.h:61
GFace.h
ForestOptions::dim
int dim
Definition: automaticMeshSizeField.h:40
GModel::getNumMeshVertices
std::size_t getNumMeshVertices(int dim=-1) const
Definition: GModel.cpp:1529
RTree::GetNext2
void GetNext2(Iterator &a_it)
Definition: rtree.h:363
MTetrahedron
Definition: MTetrahedron.h:34
curvature.h
GFace
Definition: GFace.h:33
GFace::triangles
std::vector< MTriangle * > triangles
Definition: GFace.h:428
Msg::Info
static void Info(const char *fmt,...)
Definition: GmshMessage.cpp:587
GEdge::lines
std::vector< MLine * > lines
Definition: GEdge.h:42
angle
double angle(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:157
ForestOptions::featureSizeAtVertices
std::vector< double > * featureSizeAtVertices
Definition: automaticMeshSizeField.h:49
GEntity.h
interpolation_data
Definition: automaticMeshSizeField.h:98
RTree::GetAt
DATATYPE & GetAt(Iterator &a_it)
Get object at iterator position.
Definition: rtree.h:369
robustPredicates.h
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
MVertex
Definition: MVertex.h:24
GModel::getName
std::string getName() const
Definition: GModel.h:329
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
ForestOptions::hmin
double hmin
Definition: automaticMeshSizeField.h:42
SPoint3
Definition: SPoint3.h:14
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
size_point::x
double x
Definition: automaticMeshSizeField.h:89
robustPredicates::epsilon
static REAL epsilon
Definition: robustPredicates.cpp:371
MVertex::getNum
std::size_t getNum() const
Definition: MVertex.h:86
SBoundingBox3d::min
SPoint3 min() const
Definition: SBoundingBox3d.h:90
ForestOptions::file2
FILE * file2
Definition: automaticMeshSizeField.h:60
SWAP
#define SWAP(a, b)
Definition: linearSystemCSR.cpp:15
SVector3
Definition: SVector3.h:16
automaticMeshSizeField::_smoothing
bool _smoothing
Definition: automaticMeshSizeField.h:137
MTriangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MTriangle.h:62
ForestOptions::bbox
double * bbox
Definition: automaticMeshSizeField.h:47
MEdge::tangent
SVector3 tangent() const
Definition: MEdge.h:69
ForestOptions
Definition: automaticMeshSizeField.h:39
Forest
Definition: automaticMeshSizeField.h:65
MVertex::point
SPoint3 point() const
Definition: MVertex.h:67
ForestOptions::hbulk
double hbulk
Definition: automaticMeshSizeField.h:43
GmshMessage.h
edges
static int edges[6][2]
Definition: meshGRegionLocalMeshMod.cpp:23
GEntity
Definition: GEntity.h:31
SMetric3::determinant
double determinant() const
Definition: STensor3.cpp:74
size_point::parcourus
int parcourus
Definition: automaticMeshSizeField.h:95
automaticMeshSizeField::operator()
virtual double operator()(double X, double Y, double Z, GEntity *ge=nullptr)
Definition: automaticMeshSizeField.cpp:2462
ForestOptions::nodePerTwoPi
int nodePerTwoPi
Definition: automaticMeshSizeField.h:45
automaticMeshSizeField::update
void update()
Definition: automaticMeshSizeField.cpp:2949
GModel::lastFace
fiter lastFace()
Definition: GModel.h:359
RTree< uint64_t, double, 3 >
GRegion.h
GModel::bounds
SBoundingBox3d bounds(bool aroundVisible=false)
Definition: GModel.cpp:1043
SMetric3::print
void print(const char *) const
Definition: STensor3.cpp:112
MEdge::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MEdge.h:39
ForestOptions::nodalCurvature
double * nodalCurvature
Definition: automaticMeshSizeField.h:54
ForestOptions::nodeNormals
double * nodeNormals
Definition: automaticMeshSizeField.h:55
automaticMeshSizeField::_nPointsPerCircle
int _nPointsPerCircle
Definition: automaticMeshSizeField.h:132
SBoundingBox3d.h
MVertex.h
ForestOptions::triRTree
RTree< uint64_t, double, 3 > * triRTree
Definition: automaticMeshSizeField.h:50
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
GModel::getDim
int getDim() const
Definition: GModel.cpp:989
Numeric.h
size_data::t2
SVector3 t2
Definition: automaticMeshSizeField.h:83
GModel
Definition: GModel.h:44
automaticMeshSizeField.h
size_data::h
double h
Definition: automaticMeshSizeField.h:78
GModel::firstFace
fiter firstFace()
Definition: GModel.h:355
automaticMeshSizeField::_hbulk
double _hbulk
Definition: automaticMeshSizeField.h:135
automaticMeshSizeField::_gradation
double _gradation
Definition: automaticMeshSizeField.h:136
SVector3::unit
SVector3 unit() const
Definition: SVector3.h:48
Msg::Exit
static void Exit(int level)
Definition: GmshMessage.cpp:384
automaticMeshSizeField::_hmax
double _hmax
Definition: automaticMeshSizeField.h:134
GEntity::tag
int tag() const
Definition: GEntity.h:280
sort4
void sort4(T *t[4])
Definition: MQuadrangle.h:581
size_point::size
double size
Definition: automaticMeshSizeField.h:92
size_point
Definition: automaticMeshSizeField.h:88
RTree::Insert
void Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE &a_dataId)
MTriangle
Definition: MTriangle.h:26
length
double length(Quaternion &q)
Definition: Camera.cpp:346
automaticMeshSizeField::_hmin
double _hmin
Definition: automaticMeshSizeField.h:134
automaticMeshSizeField::_nPointsPerGap
int _nPointsPerGap
Definition: automaticMeshSizeField.h:133
b1
const double b1
Definition: GaussQuadratureHex.cpp:14
size_data::M
SMetric3 M
Definition: automaticMeshSizeField.h:79
Forest::forestOptions
ForestOptions * forestOptions
Definition: automaticMeshSizeField.h:69
ForestOptions::sizeFunction
double(* sizeFunction)(double, double, double, double)
Definition: automaticMeshSizeField.h:48
GModel::lastRegion
riter lastRegion()
Definition: GModel.h:358
MTetrahedron.h
RTree::Search
int Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], bool a_resultCallback(DATATYPE a_data, void *a_context), void *a_context)
z
const double z
Definition: GaussQuadratureQuad.cpp:56
automaticMeshSizeField::~automaticMeshSizeField
~automaticMeshSizeField()
Definition: automaticMeshSizeField.cpp:2493
GEdge
Definition: GEdge.h:26
size_data::n
SVector3 n
Definition: automaticMeshSizeField.h:83
MEdge::length
double length() const
Definition: MEdge.h:76
RTree::GetFirst
void GetFirst(Iterator &a_it)
Get 'first' for iteration.
Definition: rtree.h:339
robustPredicates::orient3d
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
Definition: robustPredicates.cpp:2321
size_data::t1
SVector3 t1
Definition: automaticMeshSizeField.h:83
size_point::m
SMetric3 m
Definition: automaticMeshSizeField.h:93
size_point::isFound
bool isFound
Definition: automaticMeshSizeField.h:94
GModel.h
MVertex::y
double y() const
Definition: MVertex.h:61
size_data::hasIntersection
bool hasIntersection
Definition: automaticMeshSizeField.h:82
RTree::IsNull
bool IsNull(Iterator &a_it)
Is iterator NULL, or at end?
Definition: rtree.h:366
ForestOptions::hmax
double hmax
Definition: automaticMeshSizeField.h:41
SBoundingBox3d::max
SPoint3 max() const
Definition: SBoundingBox3d.h:91
CurvatureRusinkiewicz
bool CurvatureRusinkiewicz(const std::vector< int > &triangles, const std::vector< SPoint3 > &nodes, std::vector< std::pair< SVector3, SVector3 > > &nodalCurvatures)
Definition: curvature.cpp:172
SBoundingBox3d
Definition: SBoundingBox3d.h:21
size_data::size
double size
Definition: automaticMeshSizeField.h:74
GModel::current
static GModel * current(int index=-1)
Definition: GModel.cpp:136
size_data
Definition: automaticMeshSizeField.h:73
MVertex::x
double x() const
Definition: MVertex.h:60
size_point::y
double y
Definition: automaticMeshSizeField.h:90
SVector3::normalize
double normalize()
Definition: SVector3.h:38
faces
static int faces[4][3]
Definition: meshGRegionDelaunayInsertion.cpp:165
Field::updateNeeded
bool updateNeeded
Definition: Field.h:129