gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
OctreePost.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 #include "Octree.h"
7 #include "OctreePost.h"
8 #include "PView.h"
9 #include "PViewData.h"
10 #include "PViewDataList.h"
11 #include "PViewDataGModel.h"
12 #include "Numeric.h"
13 #include "GmshMessage.h"
14 #include "shapeFunctions.h"
15 #include "GModel.h"
16 #include "MElement.h"
17 #include "Context.h"
18 #include "SBoundingBox3d.h"
19 
20 // helper routines for list-based views
21 
22 static void minmax(int n, double *X, double *Y, double *Z, double *min,
23  double *max)
24 {
25  min[0] = X[0];
26  min[1] = Y[0];
27  min[2] = Z[0];
28  max[0] = X[0];
29  max[1] = Y[0];
30  max[2] = Z[0];
31  for(int i = 1; i < n; i++) {
32  min[0] = (X[i] < min[0]) ? X[i] : min[0];
33  min[1] = (Y[i] < min[1]) ? Y[i] : min[1];
34  min[2] = (Z[i] < min[2]) ? Z[i] : min[2];
35  max[0] = (X[i] > max[0]) ? X[i] : max[0];
36  max[1] = (Y[i] > max[1]) ? Y[i] : max[1];
37  max[2] = (Z[i] > max[2]) ? Z[i] : max[2];
38  }
39 
40  SBoundingBox3d bb(min[0], min[1], min[2], max[0], max[1], max[2]);
41 
42  if(CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS ||
43  CTX::instance()->mesh.algo2d == ALGO_2D_QUAD_QUASI_STRUCT) {
44  /* TODO FIXME: ugly temporary fix, but we need larger bbox for the
45  * guiding field sampling on curved surfaces, thicken is not sufficient */
46  bb *= 1.1;
47  bb.makeCube();
48  }
49  else {
50  bb.thicken(0.01); // make 1% thicker
51  }
52 
53  max[0] = bb.max().x();
54  max[1] = bb.max().y();
55  max[2] = bb.max().z();
56  min[0] = bb.min().x();
57  min[1] = bb.min().y();
58  min[2] = bb.min().z();
59 }
60 
61 static void centroid(int n, double *X, double *Y, double *Z, double *c)
62 {
63  const double oc = 1. / (double)n;
64  c[0] = X[0];
65  c[1] = Y[0];
66  c[2] = Z[0];
67  for(int i = 1; i < n; i++) {
68  c[0] += X[i];
69  c[1] += Y[i];
70  c[2] += Z[i];
71  }
72  c[0] *= oc;
73  c[1] *= oc;
74  c[2] *= oc;
75 }
76 
77 static void pntBB(void *a, double *min, double *max)
78 {
79  double *X = (double *)a, *Y = &X[1], *Z = &X[2];
80  minmax(1, X, Y, Z, min, max);
81 }
82 
83 static void linBB(void *a, double *min, double *max)
84 {
85  double *X = (double *)a, *Y = &X[2], *Z = &X[4];
86  minmax(2, X, Y, Z, min, max);
87 }
88 
89 static void triBB(void *a, double *min, double *max)
90 {
91  double *X = (double *)a, *Y = &X[3], *Z = &X[6];
92  minmax(3, X, Y, Z, min, max);
93 }
94 
95 static void quaBB(void *a, double *min, double *max)
96 {
97  double *X = (double *)a, *Y = &X[4], *Z = &X[8];
98  minmax(4, X, Y, Z, min, max);
99 }
100 
101 static void tetBB(void *a, double *min, double *max)
102 {
103  double *X = (double *)a, *Y = &X[4], *Z = &X[8];
104  minmax(4, X, Y, Z, min, max);
105 }
106 
107 static void hexBB(void *a, double *min, double *max)
108 {
109  double *X = (double *)a, *Y = &X[8], *Z = &X[16];
110  minmax(8, X, Y, Z, min, max);
111 }
112 
113 static void priBB(void *a, double *min, double *max)
114 {
115  double *X = (double *)a, *Y = &X[6], *Z = &X[12];
116  minmax(6, X, Y, Z, min, max);
117 }
118 
119 static void pyrBB(void *a, double *min, double *max)
120 {
121  double *X = (double *)a, *Y = &X[5], *Z = &X[10];
122  minmax(5, X, Y, Z, min, max);
123 }
124 
125 static int pntInEle(void *a, double *x) { return 1; }
126 
127 static int linInEle(void *a, double *x)
128 {
129  double *X = (double *)a, *Y = &X[2], *Z = &X[4], uvw[3];
130  line lin(X, Y, Z);
131  lin.xyz2uvw(x, uvw);
132  return lin.isInside(uvw[0], uvw[1], uvw[2]);
133 }
134 
135 static int triInEle(void *a, double *x)
136 {
137  double *X = (double *)a, *Y = &X[3], *Z = &X[6], uvw[3];
138  triangle tri(X, Y, Z);
139  tri.xyz2uvw(x, uvw);
140  return tri.isInside(uvw[0], uvw[1], uvw[2]);
141 }
142 
143 static int quaInEle(void *a, double *x)
144 {
145  double *X = (double *)a, *Y = &X[4], *Z = &X[8], uvw[3];
146  quadrangle qua(X, Y, Z);
147  qua.xyz2uvw(x, uvw);
148  return qua.isInside(uvw[0], uvw[1], uvw[2]);
149 }
150 
151 static int tetInEle(void *a, double *x)
152 {
153  double *X = (double *)a, *Y = &X[4], *Z = &X[8], uvw[3];
154  tetrahedron tet(X, Y, Z);
155  tet.xyz2uvw(x, uvw);
156  return tet.isInside(uvw[0], uvw[1], uvw[2]);
157 }
158 
159 static int hexInEle(void *a, double *x)
160 {
161  double *X = (double *)a, *Y = &X[8], *Z = &X[16], uvw[3];
162  hexahedron hex(X, Y, Z);
163  hex.xyz2uvw(x, uvw);
164  return hex.isInside(uvw[0], uvw[1], uvw[2]);
165 }
166 
167 static int priInEle(void *a, double *x)
168 {
169  double *X = (double *)a, *Y = &X[6], *Z = &X[12], uvw[3];
170  prism pri(X, Y, Z);
171  pri.xyz2uvw(x, uvw);
172  return pri.isInside(uvw[0], uvw[1], uvw[2]);
173 }
174 
175 static int pyrInEle(void *a, double *x)
176 {
177  double *X = (double *)a, *Y = &X[5], *Z = &X[10], uvw[3];
178  pyramid pyr(X, Y, Z);
179  pyr.xyz2uvw(x, uvw);
180  return pyr.isInside(uvw[0], uvw[1], uvw[2]);
181 }
182 
183 static void pntCentroid(void *a, double *x)
184 {
185  double *X = (double *)a, *Y = &X[1], *Z = &X[2];
186  centroid(1, X, Y, Z, x);
187 }
188 
189 static void linCentroid(void *a, double *x)
190 {
191  double *X = (double *)a, *Y = &X[2], *Z = &X[4];
192  centroid(2, X, Y, Z, x);
193 }
194 
195 static void triCentroid(void *a, double *x)
196 {
197  double *X = (double *)a, *Y = &X[3], *Z = &X[6];
198  centroid(3, X, Y, Z, x);
199 }
200 
201 static void quaCentroid(void *a, double *x)
202 {
203  double *X = (double *)a, *Y = &X[4], *Z = &X[8];
204  centroid(4, X, Y, Z, x);
205 }
206 
207 static void tetCentroid(void *a, double *x)
208 {
209  double *X = (double *)a, *Y = &X[4], *Z = &X[8];
210  centroid(4, X, Y, Z, x);
211 }
212 
213 static void hexCentroid(void *a, double *x)
214 {
215  double *X = (double *)a, *Y = &X[8], *Z = &X[16];
216  centroid(8, X, Y, Z, x);
217 }
218 
219 static void priCentroid(void *a, double *x)
220 {
221  double *X = (double *)a, *Y = &X[6], *Z = &X[12];
222  centroid(6, X, Y, Z, x);
223 }
224 
225 static void pyrCentroid(void *a, double *x)
226 {
227  double *X = (double *)a, *Y = &X[5], *Z = &X[10];
228  centroid(5, X, Y, Z, x);
229 }
230 
231 static void addListOfStuff(Octree *o, std::vector<double> &l, int nbelm)
232 {
233  for(std::size_t i = 0; i < l.size(); i += nbelm) Octree_Insert(&l[i], o);
234 }
235 
236 // OctreePost implementation
237 
238 OctreePost::~OctreePost()
239 {
240  Octree_Delete(_sp);
241  Octree_Delete(_vp);
242  Octree_Delete(_tp);
243  Octree_Delete(_sl);
244  Octree_Delete(_vl);
245  Octree_Delete(_tl);
246  Octree_Delete(_st);
247  Octree_Delete(_vt);
248  Octree_Delete(_tt);
249  Octree_Delete(_sq);
250  Octree_Delete(_vq);
251  Octree_Delete(_tq);
252  Octree_Delete(_ss);
253  Octree_Delete(_vs);
254  Octree_Delete(_ts);
255  Octree_Delete(_sh);
256  Octree_Delete(_vh);
257  Octree_Delete(_th);
258  Octree_Delete(_si);
259  Octree_Delete(_vi);
260  Octree_Delete(_ti);
261  Octree_Delete(_sy);
262  Octree_Delete(_vy);
263  Octree_Delete(_ty);
264 }
265 
266 OctreePost::OctreePost(PView *v)
267 {
268  _create(v->getData(true)); // use adaptive data if available
269 }
270 
271 OctreePost::OctreePost(PViewData *data) { _create(data); }
272 
273 void OctreePost::_create(PViewData *data)
274 {
275  _sp = _vp = _tp = _sl = _vl = _tl = _st = _vt = _tt = nullptr;
276  _sq = _vq = _tq = _ss = _vs = _ts = _sh = _vh = _th = nullptr;
277  _si = _vi = _ti = _sy = _vy = _ty = nullptr;
278  _theViewDataList = nullptr;
279  _theViewDataGModel = nullptr;
280 
281  _theViewDataGModel = dynamic_cast<PViewDataGModel *>(data);
282 
283  if(_theViewDataGModel) return; // the octree is already available in the model
284 
285  _theViewDataList = dynamic_cast<PViewDataList *>(data);
286 
287  if(_theViewDataList) {
288  PViewDataList *l = _theViewDataList;
289 
290  if(l->haveInterpolationMatrices() && !l->isAdapted()) {
291  Msg::Error("Cannot create octree for non-adapted high-order list-based "
292  "view: you need");
293  Msg::Error("to select 'Adapt visualization grid' first");
294  return;
295  }
296 
297  SBoundingBox3d bb = l->getBoundingBox();
298 
299  if(CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS ||
300  CTX::instance()->mesh.algo2d == ALGO_2D_QUAD_QUASI_STRUCT) {
301  /* TODO FIXME: ugly temporary fix, but we need larger bbox for the
302  * guiding field sampling on curved surfaces, thicken is not sufficient */
303  bb *= 1.1;
304  }
305  else {
306  bb.thicken(0.01); // make 1% thicker
307  }
308 
309  SPoint3 bbmin = bb.min(), bbmax = bb.max();
310  double min[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
311  double size[3] = {bbmax.x() - bbmin.x(), bbmax.y() - bbmin.y(),
312  bbmax.z() - bbmin.z()};
313  const int maxElePerBucket = 100; // memory vs. speed trade-off
314 
315  _sp =
316  Octree_Create(maxElePerBucket, min, size, pntBB, pntCentroid, pntInEle);
317  addListOfStuff(_sp, l->SP, 3 + 1 * l->getNumTimeSteps());
318  Octree_Arrange(_sp);
319  _vp =
320  Octree_Create(maxElePerBucket, min, size, pntBB, pntCentroid, pntInEle);
321  addListOfStuff(_vp, l->VP, 3 + 3 * l->getNumTimeSteps());
322  Octree_Arrange(_vp);
323  _tp =
324  Octree_Create(maxElePerBucket, min, size, pntBB, pntCentroid, pntInEle);
325  addListOfStuff(_tp, l->TP, 3 + 9 * l->getNumTimeSteps());
326  Octree_Arrange(_tp);
327 
328  _sl =
329  Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
330  addListOfStuff(_sl, l->SL, 6 + 2 * l->getNumTimeSteps());
331  Octree_Arrange(_sl);
332  _vl =
333  Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
334  addListOfStuff(_vl, l->VL, 6 + 6 * l->getNumTimeSteps());
335  Octree_Arrange(_vl);
336  _tl =
337  Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
338  addListOfStuff(_tl, l->TL, 6 + 18 * l->getNumTimeSteps());
339  Octree_Arrange(_tl);
340 
341  _st =
342  Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
343  addListOfStuff(_st, l->ST, 9 + 3 * l->getNumTimeSteps());
344  Octree_Arrange(_st);
345  _vt =
346  Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
347  addListOfStuff(_vt, l->VT, 9 + 9 * l->getNumTimeSteps());
348  Octree_Arrange(_vt);
349  _tt =
350  Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
351  addListOfStuff(_tt, l->TT, 9 + 27 * l->getNumTimeSteps());
352  Octree_Arrange(_tt);
353 
354  _sq =
355  Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
356  addListOfStuff(_sq, l->SQ, 12 + 4 * l->getNumTimeSteps());
357  Octree_Arrange(_sq);
358  _vq =
359  Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
360  addListOfStuff(_vq, l->VQ, 12 + 12 * l->getNumTimeSteps());
361  Octree_Arrange(_vq);
362  _tq =
363  Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
364  addListOfStuff(_tq, l->TQ, 12 + 36 * l->getNumTimeSteps());
365  Octree_Arrange(_tq);
366 
367  _ss =
368  Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
369  addListOfStuff(_ss, l->SS, 12 + 4 * l->getNumTimeSteps());
370  Octree_Arrange(_ss);
371  _vs =
372  Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
373  addListOfStuff(_vs, l->VS, 12 + 12 * l->getNumTimeSteps());
374  Octree_Arrange(_vs);
375  _ts =
376  Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
377  addListOfStuff(_ts, l->TS, 12 + 36 * l->getNumTimeSteps());
378  Octree_Arrange(_ts);
379 
380  _sh =
381  Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
382  addListOfStuff(_sh, l->SH, 24 + 8 * l->getNumTimeSteps());
383  Octree_Arrange(_sh);
384  _vh =
385  Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
386  addListOfStuff(_vh, l->VH, 24 + 24 * l->getNumTimeSteps());
387  Octree_Arrange(_vh);
388  _th =
389  Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
390  addListOfStuff(_th, l->TH, 24 + 72 * l->getNumTimeSteps());
391  Octree_Arrange(_th);
392 
393  _si =
394  Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
395  addListOfStuff(_si, l->SI, 18 + 6 * l->getNumTimeSteps());
396  Octree_Arrange(_si);
397  _vi =
398  Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
399  addListOfStuff(_vi, l->VI, 18 + 18 * l->getNumTimeSteps());
400  Octree_Arrange(_vi);
401  _ti =
402  Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
403  addListOfStuff(_ti, l->TI, 18 + 54 * l->getNumTimeSteps());
404  Octree_Arrange(_ti);
405 
406  _sy =
407  Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
408  addListOfStuff(_sy, l->SY, 15 + 5 * l->getNumTimeSteps());
409  Octree_Arrange(_sy);
410  _vy =
411  Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
412  addListOfStuff(_vy, l->VY, 15 + 15 * l->getNumTimeSteps());
413  Octree_Arrange(_vy);
414  _ty =
415  Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
416  addListOfStuff(_ty, l->TY, 15 + 45 * l->getNumTimeSteps());
417  Octree_Arrange(_ty);
418  }
419 }
420 
421 static void *getElement(double P[3], Octree *octree, int nbNod, int qn,
422  double *qx, double *qy, double *qz)
423 {
424  if(qn && qx && qy && qz) {
425  std::vector<void *> v;
426  Octree_SearchAll(P, octree, &v);
427  if(nbNod == qn) {
428  // try to use the value from the same geometrical element as the one
429  // provided in qx/y/z
430  double eps = CTX::instance()->geom.tolerance;
431  for(std::size_t i = 0; i < v.size(); i++) {
432  double *X = (double *)v[i], *Y = &X[qn], *Z = &X[2 * qn];
433  bool ok = true;
434  for(int j = 0; j < qn; j++) {
435  ok &= (fabs(X[j] - qx[j]) < eps && fabs(Y[j] - qy[j]) < eps &&
436  fabs(Z[j] - qz[j]) < eps);
437  }
438  if(ok) return v[i];
439  }
440  }
441  if(v.size()) return v[0];
442  }
443  else {
444  return Octree_Search(P, octree);
445  }
446  return nullptr;
447 }
448 
449 static MElement *getElement(double P[3], GModel *m, int qn, double *qx,
450  double *qy, double *qz, int dim)
451 {
452  SPoint3 pt(P);
453  if(qn && qx && qy && qz) {
454  // try to use the value from the same geometrical element as the one
455  // provided in qx/y/z
456  double eps = CTX::instance()->geom.tolerance;
457  std::vector<MElement *> elements = m->getMeshElementsByCoord(pt, dim);
458  for(std::size_t i = 0; i < elements.size(); i++) {
459  if(qn == static_cast<int>(elements[i]->getNumVertices())) {
460  bool ok = true;
461  for(int j = 0; j < qn; j++) {
462  MVertex *v = elements[i]->getVertex(j);
463  ok &=
464  (std::abs(v->x() - qx[j]) < eps && std::abs(v->y() - qy[j]) < eps &&
465  std::abs(v->z() - qz[j]) < eps);
466  }
467  if(ok) return elements[i];
468  }
469  }
470  if(elements.size()) return elements[0];
471  }
472  else {
473  SPoint3 uvw;
474  return m->getMeshElementByCoord(pt, uvw, dim);
475  }
476  return nullptr;
477 }
478 
479 bool OctreePost::_getValue(void *in, int dim, int nbNod, int nbComp,
480  double P[3], int step, double *values,
481  double *elementSize, bool grad)
482 {
483  if(!in) return false;
484 
485  double *X = (double *)in, *Y = &X[nbNod], *Z = &X[2 * nbNod],
486  *V = &X[3 * nbNod], U[3];
487 
488  elementFactory factory;
489  element *e = factory.create(nbNod, dim, X, Y, Z);
490  if(!e) return false;
491 
492  e->xyz2uvw(P, U);
493  if(step < 0) {
494  for(int i = 0; i < _theViewDataList->getNumTimeSteps(); i++) {
495  for(int j = 0; j < nbComp; j++) {
496  if(!grad) {
497  values[nbComp * i + j] = e->interpolate(&V[nbNod * nbComp * i + j],
498  U[0], U[1], U[2], nbComp);
499  }
500  else {
501  e->interpolateGrad(&V[nbNod * nbComp * i + j], U[0], U[1], U[2],
502  &values[3 * (nbComp * i + j)], nbComp);
503  }
504  }
505  }
506  }
507  else {
508  for(int j = 0; j < nbComp; j++) {
509  if(!grad) {
510  values[j] = e->interpolate(&V[nbNod * nbComp * step + j], U[0], U[1],
511  U[2], nbComp);
512  }
513  else
514  e->interpolateGrad(&V[nbNod * nbComp * step + j], U[0], U[1], U[2],
515  &values[3 * j], nbComp);
516  }
517  }
518 
519  if(elementSize) *elementSize = e->maxEdgeLength();
520 
521  delete e;
522  return true;
523 }
524 
525 bool OctreePost::_getValue(void *in, int nbComp, double P[3], int timestep,
526  double *values, double *elementSize, bool grad)
527 {
528  if(!in) return false;
529 
530  if(_theViewDataGModel->getNumComponents(0, 0, 0) != nbComp) return false;
531 
532  MElement *e = (MElement *)in;
533 
534  std::vector<int> dataIndex(e->getNumVertices());
535  if(_theViewDataGModel->getType() == PViewDataGModel::NodeData)
536  for(std::size_t i = 0; i < e->getNumVertices(); i++)
537  dataIndex[i] = e->getVertex(i)->getNum();
538  else
539  for(std::size_t i = 0; i < e->getNumVertices(); i++)
540  dataIndex[i] = e->getNum();
541 
542  double U[3];
543  e->xyz2uvw(P, U);
544 
545  std::vector<double> nodeval(e->getNumVertices() * 9);
546  for(int step = 0; step < _theViewDataGModel->getNumTimeSteps(); step++) {
547  if(!_theViewDataGModel->hasTimeStep(step)) continue;
548  if(timestep < 0 || step == timestep) {
549  for(std::size_t nod = 0; nod < e->getNumVertices(); nod++) {
550  for(int comp = 0; comp < nbComp; comp++)
551  _theViewDataGModel->getValueByIndex(step, dataIndex[nod], nod, comp,
552  nodeval[nod * nbComp + comp]);
553  }
554  for(int comp = 0; comp < nbComp; comp++) {
555  if(!grad) {
556  double val = e->interpolate(&nodeval[comp], U[0], U[1], U[2], nbComp);
557  if(timestep < 0)
558  values[nbComp * step + comp] = val;
559  else
560  values[comp] = val;
561  }
562  else {
563  if(timestep < 0)
564  e->interpolateGrad(&nodeval[comp], U[0], U[1], U[2],
565  &values[3 * (nbComp * step + comp)], nbComp);
566  else
567  e->interpolateGrad(&nodeval[comp], U[0], U[1], U[2],
568  &values[3 * comp], nbComp);
569  }
570  }
571  }
572  }
573 
574  if(elementSize) *elementSize = e->maxEdge();
575  return true;
576 }
577 
578 bool OctreePost::searchScalar(double x, double y, double z, double *values,
579  int step, double *size, int qn, double *qx,
580  double *qy, double *qz, bool grad, int dim)
581 {
582  double P[3] = {x, y, z};
583  int mult = grad ? 3 : 1;
584 
585  if(step < 0) {
586  int numSteps = 1;
587  if(_theViewDataList)
588  numSteps = _theViewDataList->getNumTimeSteps();
589  else if(_theViewDataGModel)
590  numSteps = _theViewDataGModel->getNumTimeSteps();
591  for(int i = 0; i < numSteps * mult; i++) { values[i] = 0.; }
592  }
593  else {
594  for(int i = 0; i < mult; i++) values[i] = 0.;
595  }
596 
597  if(_theViewDataList) {
598  if((dim < 0 || dim == 3) && _getValue(getElement(P, _ss, 4, qn, qx, qy, qz),
599  3, 4, 1, P, step, values, size, grad))
600  return true;
601  if((dim < 0 || dim == 3) && _getValue(getElement(P, _sh, 8, qn, qx, qy, qz),
602  3, 8, 1, P, step, values, size, grad))
603  return true;
604  if((dim < 0 || dim == 3) && _getValue(getElement(P, _si, 6, qn, qx, qy, qz),
605  3, 6, 1, P, step, values, size, grad))
606  return true;
607  if((dim < 0 || dim == 3) && _getValue(getElement(P, _sy, 5, qn, qx, qy, qz),
608  3, 5, 1, P, step, values, size, grad))
609  return true;
610  if((dim < 0 || dim == 2) && _getValue(getElement(P, _st, 3, qn, qx, qy, qz),
611  2, 3, 1, P, step, values, size, grad))
612  return true;
613  if((dim < 0 || dim == 2) && _getValue(getElement(P, _sq, 4, qn, qx, qy, qz),
614  2, 4, 1, P, step, values, size, grad))
615  return true;
616  if((dim < 0 || dim == 1) && _getValue(getElement(P, _sl, 2, qn, qx, qy, qz),
617  1, 2, 1, P, step, values, size, grad))
618  return true;
619  if((dim < 0 || dim == 0) && _getValue(getElement(P, _sp, 1, qn, qx, qy, qz),
620  0, 1, 1, P, step, values, size, grad))
621  return true;
622  }
623  else if(_theViewDataGModel) {
624  GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
625  if(m) {
626  MElement *e = getElement(P, m, qn, qx, qy, qz, dim);
627  if(_getValue(e, 1, P, step, values, size, grad)) { return true; }
628  }
629  }
630 
631  return false;
632 }
633 
634 bool OctreePost::searchVector(double x, double y, double z, double *values,
635  int step, double *size, int qn, double *qx,
636  double *qy, double *qz, bool grad, int dim)
637 {
638  double P[3] = {x, y, z};
639  int mult = grad ? 3 : 1;
640 
641  if(step < 0) {
642  int numSteps = 1;
643  if(_theViewDataList)
644  numSteps = _theViewDataList->getNumTimeSteps();
645  else if(_theViewDataGModel)
646  numSteps = _theViewDataGModel->getNumTimeSteps();
647  for(int i = 0; i < 3 * numSteps * mult; i++) values[i] = 0.;
648  }
649  else {
650  for(int i = 0; i < 3 * mult; i++) values[i] = 0.;
651  }
652 
653  if(_theViewDataList) {
654  if((dim < 0 || dim == 3) && _getValue(getElement(P, _vs, 4, qn, qx, qy, qz),
655  3, 4, 3, P, step, values, size, grad))
656  return true;
657  if((dim < 0 || dim == 3) && _getValue(getElement(P, _vh, 8, qn, qx, qy, qz),
658  3, 8, 3, P, step, values, size, grad))
659  return true;
660  if((dim < 0 || dim == 3) && _getValue(getElement(P, _vi, 6, qn, qx, qy, qz),
661  3, 6, 3, P, step, values, size, grad))
662  return true;
663  if((dim < 0 || dim == 3) && _getValue(getElement(P, _vy, 5, qn, qx, qy, qz),
664  3, 5, 3, P, step, values, size, grad))
665  return true;
666  if((dim < 0 || dim == 2) && _getValue(getElement(P, _vt, 3, qn, qx, qy, qz),
667  2, 3, 3, P, step, values, size, grad))
668  return true;
669  if((dim < 0 || dim == 2) && _getValue(getElement(P, _vq, 4, qn, qx, qy, qz),
670  2, 4, 3, P, step, values, size, grad))
671  return true;
672  if((dim < 0 || dim == 1) && _getValue(getElement(P, _vl, 2, qn, qx, qy, qz),
673  1, 2, 3, P, step, values, size, grad))
674  return true;
675  if((dim < 0 || dim == 0) && _getValue(getElement(P, _vp, 1, qn, qx, qy, qz),
676  0, 1, 3, P, step, values, size, grad))
677  return true;
678  }
679  else if(_theViewDataGModel) {
680  GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
681  if(m) {
682  MElement *e = getElement(P, m, qn, qx, qy, qz, dim);
683  if(_getValue(e, 3, P, step, values, size, grad)) { return true; }
684  }
685  }
686 
687  return false;
688 }
689 
690 bool OctreePost::searchTensor(double x, double y, double z, double *values,
691  int step, double *size, int qn, double *qx,
692  double *qy, double *qz, bool grad, int dim)
693 {
694  double P[3] = {x, y, z};
695  int mult = grad ? 3 : 1;
696 
697  if(step < 0) {
698  int numSteps = 1;
699  if(_theViewDataList)
700  numSteps = _theViewDataList->getNumTimeSteps();
701  else if(_theViewDataGModel)
702  numSteps = _theViewDataGModel->getNumTimeSteps();
703  for(int i = 0; i < 9 * numSteps * mult; i++) values[i] = 0.;
704  }
705  else {
706  for(int i = 0; i < 9 * mult; i++) values[i] = 0.;
707  }
708 
709  if(_theViewDataList) {
710  if((dim < 0 || dim == 3) && _getValue(getElement(P, _ts, 4, qn, qx, qy, qz),
711  3, 4, 9, P, step, values, size, grad))
712  return true;
713  if((dim < 0 || dim == 3) && _getValue(getElement(P, _th, 8, qn, qx, qy, qz),
714  3, 8, 9, P, step, values, size, grad))
715  return true;
716  if((dim < 0 || dim == 3) && _getValue(getElement(P, _ti, 6, qn, qx, qy, qz),
717  3, 6, 9, P, step, values, size, grad))
718  return true;
719  if((dim < 0 || dim == 3) && _getValue(getElement(P, _ty, 5, qn, qx, qy, qz),
720  3, 5, 9, P, step, values, size, grad))
721  return true;
722  if((dim < 0 || dim == 2) && _getValue(getElement(P, _tt, 3, qn, qx, qy, qz),
723  2, 3, 9, P, step, values, size, grad))
724  return true;
725  if((dim < 0 || dim == 2) && _getValue(getElement(P, _tq, 4, qn, qx, qy, qz),
726  2, 4, 9, P, step, values, size, grad))
727  return true;
728  if((dim < 0 || dim == 1) && _getValue(getElement(P, _tl, 2, qn, qx, qy, qz),
729  1, 2, 9, P, step, values, size, grad))
730  return true;
731  if((dim < 0 || dim == 0) && _getValue(getElement(P, _tp, 1, qn, qx, qy, qz),
732  0, 1, 9, P, step, values, size, grad))
733  return true;
734  }
735  else if(_theViewDataGModel) {
736  GModel *m = _theViewDataGModel->getModel((step < 0) ? 0 : step);
737  if(m) {
738  MElement *e = getElement(P, m, qn, qx, qy, qz, dim);
739  if(_getValue(e, 9, P, step, values, size, grad)) { return true; }
740  }
741  }
742 
743  return false;
744 }
MElement::getNum
virtual std::size_t getNum() const
Definition: MElement.h:68
OctreePost::_th
Octree * _th
Definition: OctreePost.h:23
quaBB
static void quaBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:95
PViewDataList::VQ
std::vector< double > VQ
Definition: PViewDataList.h:33
PViewDataList::TS
std::vector< double > TS
Definition: PViewDataList.h:37
MElement::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MElement.cpp:1128
OctreePost::~OctreePost
~OctreePost()
Definition: OctreePost.cpp:238
OctreePost::_theViewDataList
PViewDataList * _theViewDataList
Definition: OctreePost.h:26
triangle
Definition: shapeFunctions.h:414
contextGeometryOptions::tolerance
double tolerance
Definition: Context.h:99
PView
Definition: PView.h:27
OctreePost::_si
Octree * _si
Definition: OctreePost.h:24
PViewDataList::TQ
std::vector< double > TQ
Definition: PViewDataList.h:33
OctreePost::_sq
Octree * _sq
Definition: OctreePost.h:21
element::maxEdgeLength
double maxEdgeLength()
Definition: shapeFunctions.h:290
tetInEle
static int tetInEle(void *a, double *x)
Definition: OctreePost.cpp:151
OctreePost::_ts
Octree * _ts
Definition: OctreePost.h:22
PViewDataGModel
Definition: PViewDataGModel.h:191
OctreePost::searchScalar
bool searchScalar(double x, double y, double z, double *values, int step=-1, double *size=nullptr, int qn=0, double *qx=nullptr, double *qy=nullptr, double *qz=nullptr, bool grad=false, int dim=-1)
Definition: OctreePost.cpp:578
pyrCentroid
static void pyrCentroid(void *a, double *x)
Definition: OctreePost.cpp:225
PViewDataGModel::NodeData
@ NodeData
Definition: PViewDataGModel.h:194
PViewDataList::SS
std::vector< double > SS
Definition: PViewDataList.h:37
quadrangle::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:672
pyrInEle
static int pyrInEle(void *a, double *x)
Definition: OctreePost.cpp:175
triBB
static void triBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:89
PViewDataList::TL
std::vector< double > TL
Definition: PViewDataList.h:29
SBoundingBox3d::thicken
void thicken(double factor)
Definition: SBoundingBox3d.h:103
triCentroid
static void triCentroid(void *a, double *x)
Definition: OctreePost.cpp:195
OctreePost.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
PViewDataList
Definition: PViewDataList.h:17
Octree.h
MVertex
Definition: MVertex.h:24
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
OctreePost::_ti
Octree * _ti
Definition: OctreePost.h:24
PViewDataGModel::getNumComponents
int getNumComponents(int step, int ent, int ele)
Definition: PViewDataGModel.cpp:597
SPoint3
Definition: SPoint3.h:14
PViewDataList::TI
std::vector< double > TI
Definition: PViewDataList.h:41
MVertex::getNum
std::size_t getNum() const
Definition: MVertex.h:86
OctreePost::_ty
Octree * _ty
Definition: OctreePost.h:25
ALGO_2D_PACK_PRLGRMS
#define ALGO_2D_PACK_PRLGRMS
Definition: GmshDefines.h:245
SBoundingBox3d::min
SPoint3 min() const
Definition: SBoundingBox3d.h:90
Octree_Insert
void Octree_Insert(void *element, Octree *myOctree)
Definition: Octree.cpp:55
PViewDataList::SY
std::vector< double > SY
Definition: PViewDataList.h:43
line::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:404
PViewDataList::VL
std::vector< double > VL
Definition: PViewDataList.h:29
PView.h
OctreePost::_sh
Octree * _sh
Definition: OctreePost.h:23
Octree_Delete
void Octree_Delete(Octree *myOctree)
Definition: Octree.cpp:46
Octree_SearchAll
void Octree_SearchAll(double *pt, Octree *myOctree, std::vector< void * > *output)
Definition: Octree.cpp:88
PViewDataList::VY
std::vector< double > VY
Definition: PViewDataList.h:43
GmshMessage.h
OctreePost::_vh
Octree * _vh
Definition: OctreePost.h:23
PViewData.h
PViewData::haveInterpolationMatrices
bool haveInterpolationMatrices(int type=0)
Definition: PViewData.cpp:186
centroid
static void centroid(int n, double *X, double *Y, double *Z, double *c)
Definition: OctreePost.cpp:61
MElement::getNumVertices
virtual std::size_t getNumVertices() const =0
prism::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:1191
PViewDataList::VH
std::vector< double > VH
Definition: PViewDataList.h:39
minmax
static void minmax(int n, double *X, double *Y, double *Z, double *min, double *max)
Definition: OctreePost.cpp:22
OctreePost::_tp
Octree * _tp
Definition: OctreePost.h:18
elementFactory
Definition: shapeFunctions.h:1402
priBB
static void priBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:113
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
MElement::getVertex
virtual const MVertex * getVertex(int num) const =0
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
PViewDataList::TY
std::vector< double > TY
Definition: PViewDataList.h:43
OctreePost::_getValue
bool _getValue(void *in, int dim, int nbNod, int nbComp, double P[3], int step, double *values, double *elementSize, bool grad)
Definition: OctreePost.cpp:479
pntInEle
static int pntInEle(void *a, double *x)
Definition: OctreePost.cpp:125
Octree_Search
void * Octree_Search(double *pt, Octree *myOctree)
Definition: Octree.cpp:81
triInEle
static int triInEle(void *a, double *x)
Definition: OctreePost.cpp:135
Octree_Arrange
void Octree_Arrange(Octree *myOctree)
Definition: Octree.cpp:67
PViewDataList::TP
std::vector< double > TP
Definition: PViewDataList.h:27
PViewDataList::getNumTimeSteps
int getNumTimeSteps()
Definition: PViewDataList.h:79
OctreePost::_st
Octree * _st
Definition: OctreePost.h:20
MElement::interpolateGrad
void interpolateGrad(double val[], double u, double v, double w, double f[], int stride=1, double invjac[3][3]=nullptr, int order=-1)
Definition: MElement.cpp:1230
priInEle
static int priInEle(void *a, double *x)
Definition: OctreePost.cpp:167
PViewDataList::VP
std::vector< double > VP
Definition: PViewDataList.h:27
OctreePost::_vs
Octree * _vs
Definition: OctreePost.h:22
PViewDataList::SP
std::vector< double > SP
Definition: PViewDataList.h:27
SBoundingBox3d.h
prism
Definition: shapeFunctions.h:1031
Numeric.h
OctreePost::_vp
Octree * _vp
Definition: OctreePost.h:18
GModel
Definition: GModel.h:44
PView::getData
PViewData * getData(bool useAdaptiveIfAvailable=false)
Definition: PView.cpp:233
linBB
static void linBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:83
PViewDataList::SL
std::vector< double > SL
Definition: PViewDataList.h:29
PViewDataList::getBoundingBox
SBoundingBox3d getBoundingBox(int step=-1)
Definition: PViewDataList.h:87
PViewDataList.h
OctreePost::_sl
Octree * _sl
Definition: OctreePost.h:19
tetCentroid
static void tetCentroid(void *a, double *x)
Definition: OctreePost.cpp:207
hexBB
static void hexBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:107
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
Octree_Create
Octree * Octree_Create(int maxElements, double origin[3], double size[3], void(*BB)(void *, double *, double *), void(*Centroid)(void *, double *), int(*InEle)(void *, double *))
Definition: Octree.cpp:11
hexahedron::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:1021
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
MElement
Definition: MElement.h:30
PViewDataGModel.h
PViewDataList::isAdapted
bool isAdapted()
Definition: PViewDataList.h:76
quaInEle
static int quaInEle(void *a, double *x)
Definition: OctreePost.cpp:143
mult
Quaternion mult(const Quaternion &A, const Quaternion &B)
Definition: Camera.cpp:459
PViewDataGModel::hasTimeStep
bool hasTimeStep(int step)
Definition: PViewDataGModel.cpp:830
elementFactory::create
element * create(int numNodes, int dimension, double *x, double *y, double *z, bool copy=false)
Definition: shapeFunctions.h:1404
PViewDataList::ST
std::vector< double > ST
Definition: PViewDataList.h:31
PViewDataGModel::getNumTimeSteps
int getNumTimeSteps()
Definition: PViewDataGModel.cpp:270
PViewDataList::TT
std::vector< double > TT
Definition: PViewDataList.h:31
SBoundingBox3d::makeCube
void makeCube()
Definition: SBoundingBox3d.h:94
element
Definition: shapeFunctions.h:12
element::interpolateGrad
void interpolateGrad(double val[], double u, double v, double w, double f[3], int stride=1, double invjac[3][3]=nullptr)
Definition: shapeFunctions.h:167
PViewDataGModel::getType
int getType(int step, int ent, int ele)
Definition: PViewDataGModel.cpp:700
CTX::geom
contextGeometryOptions geom
Definition: Context.h:311
addListOfStuff
static void addListOfStuff(Octree *o, std::vector< double > &l, int nbelm)
Definition: OctreePost.cpp:231
linInEle
static int linInEle(void *a, double *x)
Definition: OctreePost.cpp:127
Octree
Definition: OctreeInternals.h:51
OctreePost::_vt
Octree * _vt
Definition: OctreePost.h:20
OctreePost::_tt
Octree * _tt
Definition: OctreePost.h:20
PViewData
Definition: PViewData.h:29
GModel::getMeshElementsByCoord
std::vector< MElement * > getMeshElementsByCoord(SPoint3 &p, int dim=-1, bool strict=true)
Definition: GModel.cpp:1865
triangle::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:545
tetrahedron::xyz2uvw
void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:804
PViewDataList::TH
std::vector< double > TH
Definition: PViewDataList.h:39
OctreePost::_sp
Octree * _sp
Definition: OctreePost.h:18
PViewDataList::SH
std::vector< double > SH
Definition: PViewDataList.h:39
OctreePost::_ss
Octree * _ss
Definition: OctreePost.h:22
Context.h
priCentroid
static void priCentroid(void *a, double *x)
Definition: OctreePost.cpp:219
element::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:248
OctreePost::searchVector
bool searchVector(double x, double y, double z, double *values, int step=-1, double *size=nullptr, int qn=0, double *qx=nullptr, double *qy=nullptr, double *qz=nullptr, bool grad=false, int dim=-1)
Definition: OctreePost.cpp:634
OctreePost::_theViewDataGModel
PViewDataGModel * _theViewDataGModel
Definition: OctreePost.h:27
z
const double z
Definition: GaussQuadratureQuad.cpp:56
OctreePost::_vi
Octree * _vi
Definition: OctreePost.h:24
hexCentroid
static void hexCentroid(void *a, double *x)
Definition: OctreePost.cpp:213
contextMeshOptions::algo2d
int algo2d
Definition: Context.h:29
OctreePost::_vq
Octree * _vq
Definition: OctreePost.h:21
MElement.h
tetBB
static void tetBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:101
tetrahedron
Definition: shapeFunctions.h:682
PViewDataList::VS
std::vector< double > VS
Definition: PViewDataList.h:37
MElement::maxEdge
virtual double maxEdge()
Definition: MElement.cpp:246
OctreePost::_create
void _create(PViewData *data)
Definition: OctreePost.cpp:273
quadrangle
Definition: shapeFunctions.h:554
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
OctreePost::_tl
Octree * _tl
Definition: OctreePost.h:19
PViewDataList::SI
std::vector< double > SI
Definition: PViewDataList.h:41
GModel.h
OctreePost::_tq
Octree * _tq
Definition: OctreePost.h:21
MVertex::y
double y() const
Definition: MVertex.h:61
ElementType::getNumVertices
int getNumVertices(int type)
Definition: ElementType.cpp:456
pyramid
Definition: shapeFunctions.h:1201
OctreePost::_sy
Octree * _sy
Definition: OctreePost.h:25
OctreePost::_vl
Octree * _vl
Definition: OctreePost.h:19
pntCentroid
static void pntCentroid(void *a, double *x)
Definition: OctreePost.cpp:183
shapeFunctions.h
line
Definition: shapeFunctions.h:342
OctreePost::OctreePost
OctreePost(PView *v)
Definition: OctreePost.cpp:266
pyrBB
static void pyrBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:119
linCentroid
static void linCentroid(void *a, double *x)
Definition: OctreePost.cpp:189
PViewDataList::VI
std::vector< double > VI
Definition: PViewDataList.h:41
element::interpolate
double interpolate(double val[], double u, double v, double w, int stride=1)
Definition: shapeFunctions.h:155
PViewDataList::VT
std::vector< double > VT
Definition: PViewDataList.h:31
quaCentroid
static void quaCentroid(void *a, double *x)
Definition: OctreePost.cpp:201
SBoundingBox3d::max
SPoint3 max() const
Definition: SBoundingBox3d.h:91
tetrahedron::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:822
PViewDataList::SQ
std::vector< double > SQ
Definition: PViewDataList.h:33
pyramid::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:1392
MElement::interpolate
double interpolate(double val[], double u, double v, double w, int stride=1, int order=-1)
Definition: MElement.cpp:1216
ALGO_2D_QUAD_QUASI_STRUCT
#define ALGO_2D_QUAD_QUASI_STRUCT
Definition: GmshDefines.h:247
SBoundingBox3d
Definition: SBoundingBox3d.h:21
getElement
static void * getElement(double P[3], Octree *octree, int nbNod, int qn, double *qx, double *qy, double *qz)
Definition: OctreePost.cpp:421
PViewDataGModel::getModel
GModel * getModel(int step)
Definition: PViewDataGModel.h:272
PViewDataGModel::getValueByIndex
bool getValueByIndex(int step, int dataIndex, int node, int comp, double &val)
Definition: PViewDataGModel.cpp:863
pntBB
static void pntBB(void *a, double *min, double *max)
Definition: OctreePost.cpp:77
triangle::xyz2uvw
void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:519
OctreePost::_vy
Octree * _vy
Definition: OctreePost.h:25
hexahedron
Definition: shapeFunctions.h:831
hexInEle
static int hexInEle(void *a, double *x)
Definition: OctreePost.cpp:159
MVertex::x
double x() const
Definition: MVertex.h:60
OctreePost::searchTensor
bool searchTensor(double x, double y, double z, double *values, int step=-1, double *size=nullptr, int qn=0, double *qx=nullptr, double *qy=nullptr, double *qz=nullptr, bool grad=false, int dim=-1)
Definition: OctreePost.cpp:690
GModel::getMeshElementByCoord
MElement * getMeshElementByCoord(SPoint3 &p, SPoint3 &param, int dim=-1, bool strict=true)
Definition: GModel.cpp:1842