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