1 // Gmsh - Copyright (C) 1997-2021 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 "Levelset.h"
7 #include "MakeSimplex.h"
8 #include "Numeric.h"
9 #include "Iso.h"
10 #include "adaptiveData.h"
11 #include "GmshDefines.h"
12 #include "PViewOptions.h"
13
14 static const int exn[13][12][2] = {
15 {{0, 0}}, // point
16 {{0, 1}}, // line
17 {{}}, // -
18 {{0, 1}, {0, 2}, {1, 2}}, // triangle
19 {{0, 1}, {0, 3}, {1, 2}, {2, 3}}, // quad
20 {{}}, // -
21 {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}}, // tetra
22 {{}}, // -
23 {{0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 4}, {2, 3}, {2, 4}, {3, 4}}, // pyramid
24 {{0, 1},
25 {0, 2},
26 {0, 3},
27 {1, 2},
28 {1, 4},
29 {2, 5},
30 {3, 4},
31 {3, 5},
32 {4, 5}}, // prism
33 {{}},
34 {{}}, // -
35 {{0, 1},
36 {0, 3},
37 {0, 4},
38 {1, 2},
39 {1, 5},
40 {2, 3},
41 {2, 6},
42 {3, 7},
43 {4, 5},
44 {4, 7},
45 {5, 6},
46 {6, 7}} // hexa
47 };
48
numSimplexDec(int type)49 static int numSimplexDec(int type)
50 {
51 switch(type) {
52 case TYPE_QUA: return 2;
53 case TYPE_HEX: return 6;
54 case TYPE_PRI: return 3;
55 case TYPE_PYR: return 2;
56 default: return 1;
57 }
58 }
59
getSimplexDec(int numNodes,int numEdges,int type,int i,int & n0,int & n1,int & n2,int & n3,int & nn,int & ne)60 static void getSimplexDec(int numNodes, int numEdges, int type, int i, int &n0,
61 int &n1, int &n2, int &n3, int &nn, int &ne)
62 {
63 static const int qua[2][3] = {{0, 1, 2}, {0, 2, 3}};
64 static const int hex[6][4] = {{0, 1, 3, 7}, {0, 4, 1, 7}, {1, 4, 5, 7},
65 {1, 2, 3, 7}, {1, 6, 2, 7}, {1, 5, 6, 7}};
66 static const int pri[3][4] = {{0, 1, 2, 4}, {0, 2, 4, 5}, {0, 3, 4, 5}};
67 static const int pyr[2][4] = {{0, 1, 3, 4}, {1, 2, 3, 4}};
68 switch(type) {
69 case TYPE_QUA:
70 n0 = qua[i][0];
71 n1 = qua[i][1];
72 n2 = qua[i][2];
73 nn = 3;
74 ne = 3;
75 break;
76 case TYPE_HEX:
77 n0 = hex[i][0];
78 n1 = hex[i][1];
79 n2 = hex[i][2];
80 n3 = hex[i][3];
81 nn = 4;
82 ne = 6;
83 break;
84 case TYPE_PRI:
85 n0 = pri[i][0];
86 n1 = pri[i][1];
87 n2 = pri[i][2];
88 n3 = pri[i][3];
89 nn = 4;
90 ne = 6;
91 break;
92 case TYPE_PYR:
93 n0 = pyr[i][0];
94 n1 = pyr[i][1];
95 n2 = pyr[i][2];
96 n3 = pyr[i][3];
97 nn = 4;
98 ne = 6;
99 break;
100 default:
101 n0 = 0;
102 n1 = 1;
103 n2 = 2;
104 n3 = 3;
105 nn = numNodes;
106 ne = numEdges;
107 break;
108 }
109 }
110
affect(double * xpi,double * ypi,double * zpi,double valpi[12][9],int epi[12],int i,double * xp,double * yp,double * zp,double valp[12][9],int ep[12],int j,int nb)111 static void affect(double *xpi, double *ypi, double *zpi, double valpi[12][9],
112 int epi[12], int i, double *xp, double *yp, double *zp,
113 double valp[12][9], int ep[12], int j, int nb)
114 {
115 xpi[i] = xp[j];
116 ypi[i] = yp[j];
117 zpi[i] = zp[j];
118 for(int k = 0; k < nb; k++) valpi[i][k] = valp[j][k];
119 epi[i] = ep[j];
120 }
121
removeIdenticalNodes(int * np,int numComp,double xp[12],double yp[12],double zp[12],double valp[12][9],int ep[12])122 static void removeIdenticalNodes(int *np, int numComp, double xp[12],
123 double yp[12], double zp[12],
124 double valp[12][9], int ep[12])
125 {
126 double xpi[12], ypi[12], zpi[12], valpi[12][9];
127 int epi[12];
128
129 affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, numComp);
130 int npi = 1;
131 for(int j = 1; j < *np; j++) {
132 for(int i = 0; i < npi; i++) {
133 if(fabs(xp[j] - xpi[i]) < 1.e-12 && fabs(yp[j] - ypi[i]) < 1.e-12 &&
134 fabs(zp[j] - zpi[i]) < 1.e-12) {
135 break;
136 }
137 if(i == npi - 1) {
138 affect(xpi, ypi, zpi, valpi, epi, npi, xp, yp, zp, valp, ep, j,
139 numComp);
140 npi++;
141 break;
142 }
143 }
144 }
145 for(int i = 0; i < npi; i++)
146 affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, numComp);
147 *np = npi;
148 }
149
reorderQuad(int numComp,double xp[12],double yp[12],double zp[12],double valp[12][9],int ep[12])150 static void reorderQuad(int numComp, double xp[12], double yp[12],
151 double zp[12], double valp[12][9], int ep[12])
152 {
153 double xpi[1], ypi[1], zpi[1], valpi[1][9];
154 int epi[12];
155 affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, numComp);
156 affect(xp, yp, zp, valp, ep, 3, xp, yp, zp, valp, ep, 2, numComp);
157 affect(xp, yp, zp, valp, ep, 2, xpi, ypi, zpi, valpi, epi, 0, numComp);
158 }
159
reorderPrism(int numComp,double xp[12],double yp[12],double zp[12],double valp[12][9],int ep[12],int nbCut)160 static void reorderPrism(int numComp, double xp[12], double yp[12],
161 double zp[12], double valp[12][9], int ep[12],
162 int nbCut)
163 {
164 double xpi[6], ypi[6], zpi[6], valpi[6][9];
165 int epi[12];
166
167 if(nbCut == 3) {
168 // 3 first nodes come from zero levelset intersection, next 3 are
169 // endpoints of relative edges
170 affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, numComp);
171 affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, numComp);
172 affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 5, numComp);
173 for(int i = 0; i < 3; i++) {
174 int edgecut = ep[i] - 1;
175 for(int j = 0; j < 3; j++) {
176 int p = -epi[j] - 1;
177 if(exn[9][edgecut][0] == p || exn[9][edgecut][1] == p)
178 affect(xp, yp, zp, valp, ep, 3 + i, xpi, ypi, zpi, valpi, epi, j,
179 numComp);
180 }
181 }
182 }
183 else if(nbCut == 4) {
184 // 4 first nodes come from zero levelset intersection
185 affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, numComp);
186 int edgecut = ep[0] - 1;
187 int p0 = -ep[4] - 1;
188
189 if(exn[9][edgecut][0] == p0 || exn[9][edgecut][1] == p0) {
190 affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, numComp);
191 if(exn[9][ep[1] - 1][0] == p0 || exn[9][ep[1] - 1][1] == p0) {
192 affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, numComp);
193 affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, numComp);
194 affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, numComp);
195 affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
196 }
197 else {
198 affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, numComp);
199 affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, numComp);
200 affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, numComp);
201 affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
202 }
203 }
204 else {
205 affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 5, numComp);
206 if(exn[9][ep[1] - 1][0] == p0 || exn[9][ep[1] - 1][1] == p0) {
207 affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, numComp);
208 affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, numComp);
209 affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, numComp);
210 affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
211 }
212 else {
213 affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, numComp);
214 affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, numComp);
215 affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, numComp);
216 affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
217 }
218 }
219 for(int i = 0; i < 6; i++)
220 affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, numComp);
221 }
222 }
223
GMSH_LevelsetPlugin()224 GMSH_LevelsetPlugin::GMSH_LevelsetPlugin()
225 {
226 _invert = 0.;
227 _ref[0] = _ref[1] = _ref[2] = 0.;
228 _valueIndependent = 0; // "moving" levelset
229 _valueView = -1; // use same view for levelset and field data
230 _valueTimeStep = -1; // use same time step in levelset and field data views
231 _recurLevel = 4;
232 _targetError = 0.;
233 _extractVolume =
234 0; // to create isovolumes (keep all elements < or > levelset)
235 _orientation = GMSH_LevelsetPlugin::NONE;
236 }
237
_addElement(int np,int numEdges,int numComp,double xp[12],double yp[12],double zp[12],double valp[12][9],PViewDataList * out,bool firstStep)238 void GMSH_LevelsetPlugin::_addElement(int np, int numEdges, int numComp,
239 double xp[12], double yp[12],
240 double zp[12], double valp[12][9],
241 PViewDataList *out, bool firstStep)
242 {
243 std::vector<double> *list;
244 int *nbPtr;
245 switch(np) {
246 case 1:
247 if(numComp == 1) {
248 list = &out->SP;
249 nbPtr = &out->NbSP;
250 }
251 else if(numComp == 3) {
252 list = &out->VP;
253 nbPtr = &out->NbVP;
254 }
255 else {
256 list = &out->TP;
257 nbPtr = &out->NbTP;
258 }
259 break;
260 case 2:
261 if(numComp == 1) {
262 list = &out->SL;
263 nbPtr = &out->NbSL;
264 }
265 else if(numComp == 3) {
266 list = &out->VL;
267 nbPtr = &out->NbVL;
268 }
269 else {
270 list = &out->TL;
271 nbPtr = &out->NbTL;
272 }
273 break;
274 case 3:
275 if(numComp == 1) {
276 list = &out->ST;
277 nbPtr = &out->NbST;
278 }
279 else if(numComp == 3) {
280 list = &out->VT;
281 nbPtr = &out->NbVT;
282 }
283 else {
284 list = &out->TT;
285 nbPtr = &out->NbTT;
286 }
287 break;
288 case 4:
289 if(!_extractVolume || numEdges <= 4) {
290 if(numComp == 1) {
291 list = &out->SQ;
292 nbPtr = &out->NbSQ;
293 }
294 else if(numComp == 3) {
295 list = &out->VQ;
296 nbPtr = &out->NbVQ;
297 }
298 else {
299 list = &out->TQ;
300 nbPtr = &out->NbTQ;
301 }
302 }
303 else {
304 if(numComp == 1) {
305 list = &out->SS;
306 nbPtr = &out->NbSS;
307 }
308 else if(numComp == 3) {
309 list = &out->VS;
310 nbPtr = &out->NbVS;
311 }
312 else {
313 list = &out->TS;
314 nbPtr = &out->NbTS;
315 }
316 }
317 break;
318 case 5:
319 if(numComp == 1) {
320 list = &out->SY;
321 nbPtr = &out->NbSY;
322 }
323 else if(numComp == 3) {
324 list = &out->VY;
325 nbPtr = &out->NbVY;
326 }
327 else {
328 list = &out->TY;
329 nbPtr = &out->NbTY;
330 }
331 break;
332 case 6:
333 if(numComp == 1) {
334 list = &out->SI;
335 nbPtr = &out->NbSI;
336 }
337 else if(numComp == 3) {
338 list = &out->VI;
339 nbPtr = &out->NbVI;
340 }
341 else {
342 list = &out->TI;
343 nbPtr = &out->NbTI;
344 }
345 break;
346 case 8:
347 if(numComp == 1) {
348 list = &out->SH;
349 nbPtr = &out->NbSH;
350 }
351 else if(numComp == 3) {
352 list = &out->VH;
353 nbPtr = &out->NbVH;
354 }
355 else {
356 list = &out->TH;
357 nbPtr = &out->NbTH;
358 }
359 break;
360 default: return;
361 }
362
363 // copy the elements in the output data
364 if(firstStep || !_valueIndependent) {
365 for(int k = 0; k < np; k++) list->push_back(xp[k]);
366 for(int k = 0; k < np; k++) list->push_back(yp[k]);
367 for(int k = 0; k < np; k++) list->push_back(zp[k]);
368 (*nbPtr)++;
369 }
370 for(int k = 0; k < np; k++)
371 for(int l = 0; l < numComp; l++) list->push_back(valp[k][l]);
372 }
373
_cutAndAddElements(PViewData * vdata,PViewData * wdata,int ent,int ele,int vstep,int wstep,double x[8],double y[8],double z[8],double levels[8],double scalarValues[8],PViewDataList * out)374 void GMSH_LevelsetPlugin::_cutAndAddElements(
375 PViewData *vdata, PViewData *wdata, int ent, int ele, int vstep, int wstep,
376 double x[8], double y[8], double z[8], double levels[8],
377 double scalarValues[8], PViewDataList *out)
378 {
379 int stepmin = vstep, stepmax = vstep + 1, otherstep = wstep;
380 if(stepmin < 0) {
381 stepmin = vdata->getFirstNonEmptyTimeStep();
382 stepmax = vdata->getNumTimeSteps();
383 }
384 if(wstep < 0) otherstep = wdata->getFirstNonEmptyTimeStep();
385
386 int numNodes = vdata->getNumNodes(stepmin, ent, ele);
387 int numEdges = vdata->getNumEdges(stepmin, ent, ele);
388 int numComp = wdata->getNumComponents(otherstep, ent, ele);
389 int type = vdata->getType(stepmin, ent, ele);
390
391 // decompose the element into simplices
392 for(int simplex = 0; simplex < numSimplexDec(type); simplex++) {
393 int n[4], ep[12], nsn, nse;
394 getSimplexDec(numNodes, numEdges, type, simplex, n[0], n[1], n[2], n[3],
395 nsn, nse);
396
397 // loop over time steps
398 for(int step = stepmin; step < stepmax; step++) {
399 // check which edges cut the iso and interpolate the value
400 if(wstep < 0) otherstep = step;
401
402 if(!wdata->hasTimeStep(otherstep)) continue;
403
404 int np = 0;
405 double xp[12], yp[12], zp[12], valp[12][9];
406 for(int i = 0; i < nse; i++) {
407 int n0 = exn[nse][i][0], n1 = exn[nse][i][1];
408 if(levels[n[n0]] * levels[n[n1]] <= 0.) {
409 double c = InterpolateIso(x, y, z, levels, 0., n[n0], n[n1], &xp[np],
410 &yp[np], &zp[np]);
411 for(int comp = 0; comp < numComp; comp++) {
412 double v0, v1;
413 wdata->getValue(otherstep, ent, ele, n[n0], comp, v0);
414 wdata->getValue(otherstep, ent, ele, n[n1], comp, v1);
415 valp[np][comp] = v0 + c * (v1 - v0);
416 }
417 ep[np++] = i + 1;
418 }
419 }
420
421 // remove identical nodes (this can happen if an edge actually
422 // belongs to the zero levelset, i.e., if levels[] * levels[] == 0)
423 if(np > 1) removeIdenticalNodes(&np, numComp, xp, yp, zp, valp, ep);
424
425 // if there are no cuts and we extract the volume, save the full
426 // element if it is on the correct side of the levelset
427 if(np <= 1 && _extractVolume) {
428 bool add = true;
429 for(int nod = 0; nod < nsn; nod++) {
430 if((_extractVolume < 0. && levels[n[nod]] > 0.) ||
431 (_extractVolume > 0. && levels[n[nod]] < 0.)) {
432 add = false;
433 break;
434 }
435 }
436 if(add) {
437 for(int nod = 0; nod < nsn; nod++) {
438 xp[nod] = x[n[nod]];
439 yp[nod] = y[n[nod]];
440 zp[nod] = z[n[nod]];
441 for(int comp = 0; comp < numComp; comp++)
442 wdata->getValue(otherstep, ent, ele, n[nod], comp,
443 valp[nod][comp]);
444 }
445 _addElement(nsn, nse, numComp, xp, yp, zp, valp, out,
446 step == stepmin);
447 }
448 continue;
449 }
450
451 // discard invalid cases
452 if(numEdges > 4 && np < 3) // 3D input should only lead to 2D output
453 continue;
454 else if(numEdges > 1 && np < 2) // 2D input should only lead to 1D output
455 continue;
456 else if(np < 1 || np > 4) // can't deal with this
457 continue;
458
459 // avoid "butterflies"
460 if(np == 4) reorderQuad(numComp, xp, yp, zp, valp, ep);
461
462 // orient the triangles and the quads to get the normals right
463 if(!_extractVolume && (np == 3 || np == 4)) {
464 // compute invertion test only once for spatially-fixed views
465 if(step == stepmin || !_valueIndependent) {
466 double v1[3] = {xp[2] - xp[0], yp[2] - yp[0], zp[2] - zp[0]};
467 double v2[3] = {xp[1] - xp[0], yp[1] - yp[0], zp[1] - zp[0]};
468 double gr[3], normal[3];
469 prodve(v1, v2, normal);
470 switch(_orientation) {
471 case MAP:
472 gradSimplex(x, y, z, scalarValues, gr);
473 _invert = prosca(gr, normal);
474 break;
475 case PLANE: _invert = prosca(normal, _ref); break;
476 case SPHERE:
477 gr[0] = xp[0] - _ref[0];
478 gr[1] = yp[0] - _ref[1];
479 gr[2] = zp[0] - _ref[2];
480 _invert = prosca(gr, normal);
481 case NONE:
482 default: break;
483 }
484 }
485 if(_invert > 0.) {
486 double xpi[12], ypi[12], zpi[12], valpi[12][9];
487 int epi[12];
488 for(int k = 0; k < np; k++)
489 affect(xpi, ypi, zpi, valpi, epi, k, xp, yp, zp, valp, ep, k,
490 numComp);
491 for(int k = 0; k < np; k++)
492 affect(xp, yp, zp, valp, ep, k, xpi, ypi, zpi, valpi, epi,
493 np - k - 1, numComp);
494 }
495 }
496
497 // if we extract volumes, add the nodes on the chosen side
498 // (FIXME: when cutting 2D views, the elts can have the wrong
499 // orientation)
500 if(_extractVolume) {
501 int nbCut = np;
502 for(int nod = 0; nod < nsn; nod++) {
503 if((_extractVolume < 0. && levels[n[nod]] < 0.) ||
504 (_extractVolume > 0. && levels[n[nod]] > 0.)) {
505 xp[np] = x[n[nod]];
506 yp[np] = y[n[nod]];
507 zp[np] = z[n[nod]];
508 for(int comp = 0; comp < numComp; comp++)
509 wdata->getValue(otherstep, ent, ele, n[nod], comp,
510 valp[np][comp]);
511 ep[np] = -(nod + 1); // store node num!
512 np++;
513 }
514 }
515 removeIdenticalNodes(&np, numComp, xp, yp, zp, valp, ep);
516 if(np == 4 && numEdges <= 4) reorderQuad(numComp, xp, yp, zp, valp, ep);
517 if(np == 6) reorderPrism(numComp, xp, yp, zp, valp, ep, nbCut);
518 if(np > 8) // can't deal with this
519 continue;
520 }
521
522 // finally, add the new element
523 _addElement(np, numEdges, numComp, xp, yp, zp, valp, out,
524 step == stepmin);
525 }
526
527 }
528
529 if(vstep < 0 && (stepmax - stepmin) > (int)out->Time.size()) {
530 out->Time.clear();
531 for(int i = stepmin; i < stepmax; i++) {
532 out->Time.push_back(vdata->getTime(i));
533 }
534 }
535 }
536
execute(PView * v)537 PView *GMSH_LevelsetPlugin::execute(PView *v)
538 {
539 // for adapted views we can only run the plugin on one step at a time
540 if(v->getData()->getAdaptiveData()) {
541 PViewOptions *opt = v->getOptions();
542 v->getData()->getAdaptiveData()->changeResolution(
543 opt->timeStep, _recurLevel, _targetError, this);
544 v->setChanged(true);
545 }
546
547 PViewData *vdata = getPossiblyAdaptiveData(v), *wdata;
548 if(_valueView < 0) { wdata = vdata; }
549 else if(_valueView > (int)PView::list.size() - 1) {
550 Msg::Error("View[%d] does not exist: reverting to View[%d]", _valueView,
551 v->getIndex());
552 wdata = vdata;
553 }
554 else {
555 wdata = getPossiblyAdaptiveData(PView::list[_valueView]);
556 }
557
558 // sanity checks
559 if(vdata->getNumEntities() != wdata->getNumEntities() ||
560 vdata->getNumElements() != wdata->getNumElements()) {
561 Msg::Error("Incompatible views");
562 return v;
563 }
564 if(_valueTimeStep >= wdata->getNumTimeSteps()) {
565 Msg::Error("Wrong time step %d in view", _valueTimeStep);
566 return v;
567 }
568
569 // Force creation of one view per time step if we have multi meshes
570 if(vdata->hasMultipleMeshes()) _valueIndependent = 0;
571
572 double x[8], y[8], z[8], levels[8];
573 double scalarValues[8] = {0., 0., 0., 0., 0., 0., 0., 0.};
574
575 PView *v2 = nullptr;
576 if(_valueIndependent) {
577 // create a single output view containing the (possibly multi-step) levelset
578 int firstNonEmptyStep = vdata->getFirstNonEmptyTimeStep();
579 v2 = new PView();
580 PViewDataList *out = getDataList(v2);
581 for(int ent = 0; ent < vdata->getNumEntities(firstNonEmptyStep); ent++) {
582 for(int ele = 0; ele < vdata->getNumElements(firstNonEmptyStep, ent);
583 ele++) {
584 if(vdata->skipElement(firstNonEmptyStep, ent, ele)) continue;
585 for(int nod = 0; nod < vdata->getNumNodes(firstNonEmptyStep, ent, ele);
586 nod++) {
587 vdata->getNode(firstNonEmptyStep, ent, ele, nod, x[nod], y[nod],
588 z[nod]);
589 levels[nod] = levelset(x[nod], y[nod], z[nod], 0.);
590 }
591 _cutAndAddElements(vdata, wdata, ent, ele, -1, _valueTimeStep, x, y, z,
592 levels, scalarValues, out);
593 }
594 }
595 out->setName(vdata->getName() + "_Levelset");
596 out->setFileName(vdata->getFileName() + "_Levelset.pos");
597 out->finalize();
598 }
599 else {
600 // create one view per timestep
601 for(int step = 0; step < vdata->getNumTimeSteps(); step++) {
602 if(!vdata->hasTimeStep(step)) continue;
603 v2 = new PView();
604 PViewDataList *out = getDataList(v2);
605 for(int ent = 0; ent < vdata->getNumEntities(step); ent++) {
606 for(int ele = 0; ele < vdata->getNumElements(step, ent); ele++) {
607 if(vdata->skipElement(step, ent, ele)) continue;
608 for(int nod = 0; nod < vdata->getNumNodes(step, ent, ele); nod++) {
609 vdata->getNode(step, ent, ele, nod, x[nod], y[nod], z[nod]);
610 vdata->getScalarValue(step, ent, ele, nod, scalarValues[nod]);
611 levels[nod] = levelset(x[nod], y[nod], z[nod], scalarValues[nod]);
612 }
613 int wstep = (_valueTimeStep < 0) ? step : _valueTimeStep;
614 _cutAndAddElements(vdata, wdata, ent, ele, step, wstep, x, y, z,
615 levels, scalarValues, out);
616 }
617 }
618 char tmp[246];
619 sprintf(tmp, "_Levelset_%d", step);
620 out->setName(vdata->getName() + tmp);
621 out->setFileName(vdata->getFileName() + tmp + ".pos");
622 out->finalize();
623 }
624 }
625
626 return v2;
627 }
628
629 // On high order maps, we draw only the elements that have a cut with
630 // the levelset, this is as accurate as it should be
631
recur_sign_change(adaptiveTriangle * t,const GMSH_LevelsetPlugin * plug)632 static bool recur_sign_change(adaptiveTriangle *t,
633 const GMSH_LevelsetPlugin *plug)
634 {
635 if(!t->e[0] || t->visible) {
636 double v1 =
637 plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
638 double v2 =
639 plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
640 double v3 =
641 plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
642 if(v1 * v2 > 0 && v1 * v3 > 0)
643 t->visible = false;
644 else
645 t->visible = true;
646 return t->visible;
647 }
648 else {
649 bool sc1 = recur_sign_change(t->e[0], plug);
650 bool sc2 = recur_sign_change(t->e[1], plug);
651 bool sc3 = recur_sign_change(t->e[2], plug);
652 bool sc4 = recur_sign_change(t->e[3], plug);
653 if(sc1 || sc2 || sc3 || sc4) {
654 if(!sc1) t->e[0]->visible = true;
655 if(!sc2) t->e[1]->visible = true;
656 if(!sc3) t->e[2]->visible = true;
657 if(!sc4) t->e[3]->visible = true;
658 return true;
659 }
660 t->visible = false;
661 return false;
662 }
663 }
664
recur_sign_change(adaptiveQuadrangle * q,const GMSH_LevelsetPlugin * plug)665 static bool recur_sign_change(adaptiveQuadrangle *q,
666 const GMSH_LevelsetPlugin *plug)
667 {
668 if(!q->e[0] || q->visible) {
669 double v1 =
670 plug->levelset(q->p[0]->X, q->p[0]->Y, q->p[0]->Z, q->p[0]->val);
671 double v2 =
672 plug->levelset(q->p[1]->X, q->p[1]->Y, q->p[1]->Z, q->p[1]->val);
673 double v3 =
674 plug->levelset(q->p[2]->X, q->p[2]->Y, q->p[2]->Z, q->p[2]->val);
675 double v4 =
676 plug->levelset(q->p[3]->X, q->p[3]->Y, q->p[3]->Z, q->p[3]->val);
677 if(v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0)
678 q->visible = false;
679 else
680 q->visible = true;
681 return q->visible;
682 }
683 else {
684 bool sc1 = recur_sign_change(q->e[0], plug);
685 bool sc2 = recur_sign_change(q->e[1], plug);
686 bool sc3 = recur_sign_change(q->e[2], plug);
687 bool sc4 = recur_sign_change(q->e[3], plug);
688 if(sc1 || sc2 || sc3 || sc4) {
689 if(!sc1) q->e[0]->visible = true;
690 if(!sc2) q->e[1]->visible = true;
691 if(!sc3) q->e[2]->visible = true;
692 if(!sc4) q->e[3]->visible = true;
693 return true;
694 }
695 q->visible = false;
696 return false;
697 }
698 }
699
recur_sign_change(adaptiveTetrahedron * t,const GMSH_LevelsetPlugin * plug)700 static bool recur_sign_change(adaptiveTetrahedron *t,
701 const GMSH_LevelsetPlugin *plug)
702 {
703 if(!t->e[0] || t->visible) {
704 double v1 =
705 plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
706 double v2 =
707 plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
708 double v3 =
709 plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
710 double v4 =
711 plug->levelset(t->p[3]->X, t->p[3]->Y, t->p[3]->Z, t->p[3]->val);
712 if(v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0)
713 t->visible = false;
714 else
715 t->visible = true;
716 return t->visible;
717 }
718 else {
719 bool sc1 = recur_sign_change(t->e[0], plug);
720 bool sc2 = recur_sign_change(t->e[1], plug);
721 bool sc3 = recur_sign_change(t->e[2], plug);
722 bool sc4 = recur_sign_change(t->e[3], plug);
723 bool sc5 = recur_sign_change(t->e[4], plug);
724 bool sc6 = recur_sign_change(t->e[5], plug);
725 bool sc7 = recur_sign_change(t->e[6], plug);
726 bool sc8 = recur_sign_change(t->e[7], plug);
727 if(sc1 || sc2 || sc3 || sc4 || sc5 || sc6 || sc7 || sc8) {
728 if(!sc1) t->e[0]->visible = true;
729 if(!sc2) t->e[1]->visible = true;
730 if(!sc3) t->e[2]->visible = true;
731 if(!sc4) t->e[3]->visible = true;
732 if(!sc5) t->e[4]->visible = true;
733 if(!sc6) t->e[5]->visible = true;
734 if(!sc7) t->e[6]->visible = true;
735 if(!sc8) t->e[7]->visible = true;
736 return true;
737 }
738 t->visible = false;
739 return false;
740 }
741 }
742
recur_sign_change(adaptiveHexahedron * t,const GMSH_LevelsetPlugin * plug)743 static bool recur_sign_change(adaptiveHexahedron *t,
744 const GMSH_LevelsetPlugin *plug)
745 {
746 if(!t->e[0] || t->visible) {
747 double v1 =
748 plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
749 double v2 =
750 plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
751 double v3 =
752 plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
753 double v4 =
754 plug->levelset(t->p[3]->X, t->p[3]->Y, t->p[3]->Z, t->p[3]->val);
755 double v5 =
756 plug->levelset(t->p[4]->X, t->p[4]->Y, t->p[4]->Z, t->p[4]->val);
757 double v6 =
758 plug->levelset(t->p[5]->X, t->p[5]->Y, t->p[5]->Z, t->p[5]->val);
759 double v7 =
760 plug->levelset(t->p[6]->X, t->p[6]->Y, t->p[6]->Z, t->p[6]->val);
761 double v8 =
762 plug->levelset(t->p[7]->X, t->p[7]->Y, t->p[7]->Z, t->p[7]->val);
763 if(v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0 && v1 * v5 > 0 &&
764 v1 * v6 > 0 && v1 * v7 > 0 && v1 * v8 > 0)
765 t->visible = false;
766 else
767 t->visible = true;
768 return t->visible;
769 }
770 else {
771 bool sc1 = recur_sign_change(t->e[0], plug);
772 bool sc2 = recur_sign_change(t->e[1], plug);
773 bool sc3 = recur_sign_change(t->e[2], plug);
774 bool sc4 = recur_sign_change(t->e[3], plug);
775 bool sc5 = recur_sign_change(t->e[4], plug);
776 bool sc6 = recur_sign_change(t->e[5], plug);
777 bool sc7 = recur_sign_change(t->e[6], plug);
778 bool sc8 = recur_sign_change(t->e[7], plug);
779 if(sc1 || sc2 || sc3 || sc4 || sc5 || sc6 || sc7 || sc8) {
780 if(!sc1) t->e[0]->visible = true;
781 if(!sc2) t->e[1]->visible = true;
782 if(!sc3) t->e[2]->visible = true;
783 if(!sc4) t->e[3]->visible = true;
784 if(!sc5) t->e[4]->visible = true;
785 if(!sc6) t->e[5]->visible = true;
786 if(!sc7) t->e[6]->visible = true;
787 if(!sc8) t->e[7]->visible = true;
788 return true;
789 }
790 t->visible = false;
791 return false;
792 }
793 }
794
recur_sign_change(adaptivePrism * t,const GMSH_LevelsetPlugin * plug)795 static bool recur_sign_change(adaptivePrism *t, const GMSH_LevelsetPlugin *plug)
796 {
797 if(!t->e[0] || t->visible) {
798 double v1 =
799 plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
800 double v2 =
801 plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
802 double v3 =
803 plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
804 double v4 =
805 plug->levelset(t->p[3]->X, t->p[3]->Y, t->p[3]->Z, t->p[3]->val);
806 double v5 =
807 plug->levelset(t->p[4]->X, t->p[4]->Y, t->p[4]->Z, t->p[4]->val);
808 double v6 =
809 plug->levelset(t->p[5]->X, t->p[5]->Y, t->p[5]->Z, t->p[5]->val);
810 if(v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0 && v1 * v5 > 0 && v1 * v6 > 0)
811 t->visible = false;
812 else
813 t->visible = true;
814 return t->visible;
815 }
816 else {
817 bool sc1 = recur_sign_change(t->e[0], plug);
818 bool sc2 = recur_sign_change(t->e[1], plug);
819 bool sc3 = recur_sign_change(t->e[2], plug);
820 bool sc4 = recur_sign_change(t->e[3], plug);
821 bool sc5 = recur_sign_change(t->e[4], plug);
822 bool sc6 = recur_sign_change(t->e[5], plug);
823 bool sc7 = recur_sign_change(t->e[6], plug);
824 bool sc8 = recur_sign_change(t->e[7], plug);
825 if(sc1 || sc2 || sc3 || sc4 || sc5 || sc6 || sc7 || sc8) {
826 if(!sc1) t->e[0]->visible = true;
827 if(!sc2) t->e[1]->visible = true;
828 if(!sc3) t->e[2]->visible = true;
829 if(!sc4) t->e[3]->visible = true;
830 if(!sc5) t->e[4]->visible = true;
831 if(!sc6) t->e[5]->visible = true;
832 if(!sc7) t->e[6]->visible = true;
833 if(!sc8) t->e[7]->visible = true;
834 return true;
835 }
836 t->visible = false;
837 return false;
838 }
839 }
840
recur_sign_change(adaptivePyramid * t,const GMSH_LevelsetPlugin * plug)841 static bool recur_sign_change(adaptivePyramid *t,
842 const GMSH_LevelsetPlugin *plug)
843 {
844 if(!t->e[0] || t->visible) {
845 double v1 =
846 plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
847 double v2 =
848 plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
849 double v3 =
850 plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
851 double v4 =
852 plug->levelset(t->p[3]->X, t->p[3]->Y, t->p[3]->Z, t->p[3]->val);
853 double v5 =
854 plug->levelset(t->p[4]->X, t->p[4]->Y, t->p[4]->Z, t->p[4]->val);
855 if(v1 * v2 > 0 && v1 * v3 > 0 && v1 * v4 > 0 && v1 * v5 > 0)
856 t->visible = false;
857 else
858 t->visible = true;
859 return t->visible;
860 }
861 else {
862 bool sc1 = recur_sign_change(t->e[0], plug);
863 bool sc2 = recur_sign_change(t->e[1], plug);
864 bool sc3 = recur_sign_change(t->e[2], plug);
865 bool sc4 = recur_sign_change(t->e[3], plug);
866 bool sc5 = recur_sign_change(t->e[4], plug);
867 bool sc6 = recur_sign_change(t->e[5], plug);
868 bool sc7 = recur_sign_change(t->e[6], plug);
869 bool sc8 = recur_sign_change(t->e[7], plug);
870 bool sc9 = recur_sign_change(t->e[8], plug);
871 bool sc10 = recur_sign_change(t->e[9], plug);
872 if(sc1 || sc2 || sc3 || sc4 || sc5 || sc6 || sc7 || sc8 || sc9 || sc10) {
873 if(!sc1) t->e[0]->visible = true;
874 if(!sc2) t->e[1]->visible = true;
875 if(!sc3) t->e[2]->visible = true;
876 if(!sc4) t->e[3]->visible = true;
877 if(!sc5) t->e[4]->visible = true;
878 if(!sc6) t->e[5]->visible = true;
879 if(!sc7) t->e[6]->visible = true;
880 if(!sc8) t->e[7]->visible = true;
881 if(!sc9) t->e[8]->visible = true;
882 if(!sc10) t->e[9]->visible = true;
883 return true;
884 }
885 t->visible = false;
886 return false;
887 }
888 }
889
assignSpecificVisibility() const890 void GMSH_LevelsetPlugin::assignSpecificVisibility() const
891 {
892 if(adaptiveTriangle::all.size()) {
893 adaptiveTriangle *t = *adaptiveTriangle::all.begin();
894 if(!t->visible) t->visible = !recur_sign_change(t, this);
895 }
896 if(adaptiveQuadrangle::all.size()) {
897 adaptiveQuadrangle *q = *adaptiveQuadrangle::all.begin();
898 if(!q->visible) q->visible = !recur_sign_change(q, this);
899 }
900 if(adaptiveTetrahedron::all.size()) {
901 adaptiveTetrahedron *t = *adaptiveTetrahedron::all.begin();
902 if(!t->visible) t->visible = !recur_sign_change(t, this);
903 }
904 if(adaptiveHexahedron::all.size()) {
905 adaptiveHexahedron *h = *adaptiveHexahedron::all.begin();
906 if(!h->visible) h->visible = !recur_sign_change(h, this);
907 }
908 if(adaptivePrism::all.size()) {
909 adaptivePrism *p = *adaptivePrism::all.begin();
910 if(!p->visible) p->visible = !recur_sign_change(p, this);
911 }
912 if(adaptivePyramid::all.size()) {
913 adaptivePyramid *p = *adaptivePyramid::all.begin();
914 if(!p->visible) p->visible = !recur_sign_change(p, this);
915 }
916 }
917