Salome HOME
36bcae4b0853d3aea47401bb2d916d310aa71342
[modules/smesh.git] / src / Tools / MeshCut / MeshCut_DC.cxx
1 // Copyright (C) 2006-2012  EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 // Classes et fonctions XMeshLab
21
22 #include "MeshCut_Utils.hxx"
23 #include "MeshCut_Maillage.hxx"
24
25 #include "MeshCut_Carre.hxx"
26 #include "MeshCut_Cube.hxx"
27
28 #include "MeshCut_Fonctions.hxx"
29 #include "MeshCut_Cas.hxx"
30
31 #include "MeshCut_Globals.hxx"
32
33 #include <iostream>
34 #include <cmath>
35 #include <cstdlib>
36 #include <cstring>
37
38 using namespace MESHCUT;
39 using namespace std;
40
41 // ==================================  DECLARATION DES VARIABLES GLOBALES  ==================================================
42
43 std::map<std::string, int> MESHCUT::intersections;
44
45 int MESHCUT::indexNouvellesMailles, MESHCUT::indexNouveauxNoeuds, MESHCUT::offsetMailles;
46 std::string MESHCUT::str_id_GMplus, MESHCUT::str_id_GMmoins;
47 Maillage *MESHCUT::MAILLAGE1, *MESHCUT::MAILLAGE2;
48
49 std::vector<float> MESHCUT::newXX, MESHCUT::newYY, MESHCUT::newZZ;
50 std::map<TYPE_MAILLE, std::vector<int> > MESHCUT::newCNX;
51 std::map<TYPE_MAILLE, int> MESHCUT::cptNouvellesMailles;
52 std::map<TYPE_MAILLE, std::vector<int> > MESHCUT::GMplus, MESHCUT::GMmoins;
53 std::vector<int> MESHCUT::cutTetras;
54
55 float *MESHCUT::DNP;
56 int *MESHCUT::POSN;
57
58 std::string MESHCUT::str_id_maillagenew;
59
60 float MESHCUT::normale[3], MESHCUT::pointPlan[3];
61 float MESHCUT::d;
62 float MESHCUT::epsilon;
63
64 bool MESHCUT::debug;
65 int MESHCUT::Naretes;
66
67 // ==================================   PROGRAMME PRINCIPAL  ==================================================
68
69 int main(int argc, char *argv[])
70 {
71
72   debug = false;
73   string ficMEDin;
74   string ficMEDout;
75   float xNormal;
76   float yNormal;
77   float zNormal;
78   float xm;
79   float ym;
80   float zm;
81   float tolerance;
82   try
83     {
84       if (argc != 13)
85         throw std::exception();
86       char *ficMEDin0 = argv[1];
87       ficMEDin = (string) ficMEDin0;
88       char *ficMEDout0 = argv[2];
89       ficMEDout = (string) ficMEDout0;
90       char *id_maillagenew = argv[3];
91       str_id_maillagenew = (string) id_maillagenew;
92
93       // Groupes créés
94       char *id_GMplus = argv[4];
95       str_id_GMplus = (string) id_GMplus;
96       char *id_GMmoins = argv[5];
97       str_id_GMmoins = (string) id_GMmoins;
98
99       // Vecteur normal au plan de coupe
100       char *charxn = argv[6];
101       xNormal = char2float(charxn);
102       char *charyn = argv[7];
103       yNormal = char2float(charyn);
104       char *charzn = argv[8];
105       zNormal = char2float(charzn);
106
107       // Point du plan de coupe
108       char *charxm = argv[9];
109       xm = char2float(charxm);
110       char *charym = argv[10];
111       ym = char2float(charym);
112       char *charzm = argv[11];
113       zm = char2float(charzm);
114
115       // Tolérance :  epsilon = tolérance * longueur arête moyenne - où epsilon est la tolérance absolue (distance)
116       char *chtolerance = argv[12];
117       tolerance = char2float(chtolerance);
118     }
119   catch (...)
120     {
121       cout << endl;
122       cout << "                 Cut a tetrahedron mesh by a plane" << endl;
123       cout << "                 ---------------------------------" << endl;
124       cout << "Syntax:" << endl << endl;
125       cout << argv[0] << " input.med output.med resuMeshName aboveGroup belowGroup nx ny nz px py pz T " << endl;
126       cout << endl << "where:" << endl;
127       cout << "  input.med    = name of the original mesh file in med format" << endl;
128       cout << "  output.med   = name of the result mesh file in med format" << endl;
129       cout << "  resuMeshName = name of the result mesh" << endl;
130       cout << "  aboveGroup   = name of the group of volumes above the cut plane" << endl;
131       cout << "  belowGroups  = name of the group of volumes below the cut plane" << endl;
132       cout << "  nx ny nz     = vector normal to the cut plane" << endl;
133       cout << "  px py pz     = a point of the cut plane" << endl;
134       cout << "  T            = 0 < T < 1 : vertices of a tetrahedron are considered as belonging" << endl;
135       cout << "                 the cut plane if their distance to the plane is inferior to L*T" << endl;
136       cout << "                 where L is the mean edge size of the tetrahedron" << endl;
137       ERREUR("--> check arguments!");
138     }
139
140   cout << "Cut by a plane :" << endl;
141   cout << "  source mesh: " << ficMEDin << endl;
142   cout << "  result mesh:  " << ficMEDout << endl;
143   cout << "  mesh name:  " << str_id_maillagenew << endl;
144   cout << "  group above plane:  " << str_id_GMplus << endl;
145   cout << "  group below plane: " << str_id_GMmoins << endl;
146   cout << "  vector normal to the cut plane: xn=" << xNormal << " yn=" << yNormal << " zn=" << zNormal << endl;
147   cout << "  point in the cut plane: xm=" << xm << " ym=" << ym << " zm=" << zm << endl;
148   cout << "  tolerance: " << tolerance << endl;
149   cout << endl;
150
151   if (tolerance <= 0.0)
152     ERREUR("Tolerance must not be negative or null");
153
154   // Il faut normer la normale
155   float normeNormal = sqrt(xNormal * xNormal + yNormal * yNormal + zNormal * zNormal);
156   if (normeNormal == 0.0)
157     ERREUR("null normal vector");
158   normale[0] = xNormal / normeNormal;
159   normale[1] = yNormal / normeNormal;
160   normale[2] = zNormal / normeNormal;
161
162   pointPlan[0] = xm;
163   pointPlan[1] = ym;
164   pointPlan[2] = zm;
165
166   // Calcul du coefficient d de l'équation du plan  xn x + yn y + zn n + d = 0
167   d = -normale[0] * xm - normale[1] * ym - normale[2] * zm;
168
169   intersections.clear();
170
171   // Initialisation des compteurs de nouvelles mailles
172   for (int itm = (int) POI1; itm <= (int) HEXA20; itm++)
173     {
174       TYPE_MAILLE tm = (TYPE_MAILLE) itm;
175       cptNouvellesMailles[tm] = 0;
176     }
177
178   int V[6];
179   int S[4]; // Signature du T4 courant
180   //int NG[4]; // Num. globaux des sommets
181
182   // Acquisition maillage initial
183   //cout << chrono() << " - Acquisition du maillage initial" << endl;
184   MAILLAGE1 = new Maillage((string) "TEMP");
185   MAILLAGE1->inputMED(ficMEDin);
186   cout << chrono() << " - End of mesh read" << endl;
187   indexNouveauxNoeuds = MAILLAGE1->nombreNoeudsMaillage;
188
189   // Le maillage ne contient aucun TETRA4 : on rend le maillage initial sans modification
190   if (!MAILLAGE1->EFFECTIFS_TYPES[TETRA4])
191     {
192       cout << "WARNING: mesh does not contain tetra4 elements, it will not be modified" << endl;
193       MAILLAGE1->ID = str_id_maillagenew;
194       MAILLAGE1->outputMED(ficMEDout);
195       cout << chrono() << " - Finished!" << endl << endl;
196       exit(0);
197     }
198   // A partir de cet instant le maillage contient forcément des TETRA4
199
200
201   // Chargement des distances noeud-plan DNP
202   DNP = (float*) malloc(sizeof(float) * MAILLAGE1->nombreNoeudsMaillage);
203   for (int k = 0; k < MAILLAGE1->nombreNoeudsMaillage; k++)
204     DNP[k] = distanceNoeudPlan(k + 1);
205   cout << chrono() << " - End of computation of distances between nodes and plane" << endl;
206
207   // Longueur d'arête moyenne des T4 intersectant le plan de coupe
208   float LONGUEURS = 0.0;
209   int cptLONGUEURS = 0;
210   for (int it4 = 0; it4 < MAILLAGE1->EFFECTIFS_TYPES[TETRA4]; it4++)
211     {
212       bool plus = false;
213       bool moins = false;
214       int *offset = MAILLAGE1->CNX[TETRA4] + 4 * it4;
215       for (int is = 0; is < 4; is++)
216         {
217           int ng = *(offset + is);
218           if (DNP[ng - 1] > 0.0)
219             plus = true;
220           else if (DNP[ng - 1] < 0.0)
221             moins = true;
222         }
223       if (plus && moins)
224         {
225           // Ce tetra est à cheval sur le plan de coupe: on calcule ses longueurs d'arêtes
226           LONGUEURS += longueurSegment(*(offset + 0), *(offset + 1));
227           cptLONGUEURS++;
228           LONGUEURS += longueurSegment(*(offset + 0), *(offset + 2));
229           cptLONGUEURS++;
230           LONGUEURS += longueurSegment(*(offset + 0), *(offset + 3));
231           cptLONGUEURS++;
232           LONGUEURS += longueurSegment(*(offset + 1), *(offset + 2));
233           cptLONGUEURS++;
234           LONGUEURS += longueurSegment(*(offset + 1), *(offset + 3));
235           cptLONGUEURS++;
236           LONGUEURS += longueurSegment(*(offset + 2), *(offset + 3));
237           cptLONGUEURS++;
238         }
239     }
240
241   // Aucun TETRA4 intercepté par le plan de coupe : on rend MAILLAGE1
242   if (cptLONGUEURS == 0)
243     {
244       cout
245           << "WARNING: the cut plane does not cut any tetra4 element, initial mesh will not be modified"
246           << endl;
247       MAILLAGE1->ID = str_id_maillagenew;
248       MAILLAGE1->outputMED(ficMEDout);
249       cout << chrono() << " - Finished!" << endl << endl;
250       exit(0);
251     }
252   // A partir de cet instant le maillage contient forcément des TETRA4 intersectant le plan de coupe
253
254
255   float longueurMoyenne = LONGUEURS / cptLONGUEURS;
256   epsilon = tolerance * longueurMoyenne;
257
258   int nT4coupe = cptLONGUEURS / 6;
259   cout << chrono() << " - End of computation of mean length of tetra4 edges near the cut plane" << endl;
260
261   cout << "Number of tetra4 to be cut = " << nT4coupe << endl;
262   cout << "Mean length = " << longueurMoyenne << endl;
263   cout << "Tolerance = " << tolerance << endl;
264   cout << "Epsilon = " << epsilon << endl;
265
266   // Détermination des positions de noeuds par rapport au plan de coupe - POSN
267   POSN = (int*) malloc(sizeof(int) * MAILLAGE1->nombreNoeudsMaillage);
268   for (int k = 0; k < MAILLAGE1->nombreNoeudsMaillage; k++)
269     {
270       if (DNP[k] > epsilon)
271         POSN[k] = 1;
272       else if (DNP[k] < -epsilon)
273         POSN[k] = -1;
274       else
275         POSN[k] = 0;
276     }
277   cout << chrono() << " - End of nodes qualification above or below the cut plane" << endl;
278   cout << "Start of iteration on tetra4" << endl;
279
280   for (int it4 = 0; it4 < MAILLAGE1->EFFECTIFS_TYPES[TETRA4]; it4++)
281     {
282
283       for (int is = 0; is < 4; is++)
284         {
285           int ng = *(MAILLAGE1->CNX[TETRA4] + 4 * it4 + is);
286           //NG[is] = ng;
287           S[is] = *(POSN + ng - 1);
288         }
289
290       // -------------------------------------------------------------------
291
292       if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == -1)
293         GMmoins[TETRA4].push_back(it4);
294
295       else if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == 0)
296         GMmoins[TETRA4].push_back(it4);
297
298       else if (S[0] == -1 && S[1] == -1 && S[2] == -1 && S[3] == 1)
299         { // Cas 3 - Arêtes 2 4 5
300           V[0] = -1;
301           V[1] = -1;
302           V[2] = intersectionSegmentPlan(it4, 2);
303           V[3] = -1;
304           V[4] = intersectionSegmentPlan(it4, 4);
305           V[5] = intersectionSegmentPlan(it4, 5);
306           cas3(V, it4);
307         }
308
309       // -------------------------------------------------------------------
310
311       else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == -1)
312         GMmoins[TETRA4].push_back(it4);
313
314       else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == 0)
315         GMmoins[TETRA4].push_back(it4);
316
317       else if (S[0] == -1 && S[1] == -1 && S[2] == 0 && S[3] == 1)
318         { // Cas 2, arêtes 2 4
319           V[0] = -1;
320           V[1] = -1;
321           V[2] = intersectionSegmentPlan(it4, 2);
322           V[3] = -1;
323           V[4] = intersectionSegmentPlan(it4, 4);
324           V[5] = -1;
325           cas2(V, it4);
326         }
327
328       // -------------------------------------------------------------------
329
330       else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == -1)
331         { // Cas 3, arêtes 1 3 5
332           V[0] = -1;
333           V[1] = intersectionSegmentPlan(it4, 1);
334           V[2] = -1;
335           V[3] = intersectionSegmentPlan(it4, 3);
336           V[4] = -1;
337           V[5] = intersectionSegmentPlan(it4, 5);
338           cas3(V, it4);
339         }
340
341       else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == 0)
342         { // Cas 2, arêtes 1 3
343           V[0] = -1;
344           V[1] = intersectionSegmentPlan(it4, 1);
345           V[2] = -1;
346           V[3] = intersectionSegmentPlan(it4, 3);
347           V[4] = -1;
348           V[5] = -1;
349           cas2(V, it4);
350         }
351
352       else if (S[0] == -1 && S[1] == -1 && S[2] == 1 && S[3] == 1)
353         { // Cas 4, arêtes 1 2 3 4
354           V[0] = -1;
355           V[1] = intersectionSegmentPlan(it4, 1);
356           V[2] = intersectionSegmentPlan(it4, 2);
357           V[3] = intersectionSegmentPlan(it4, 3);
358           V[4] = intersectionSegmentPlan(it4, 4);
359           V[5] = -1;
360           cas4(V, it4);
361         }
362
363       // -------------------------------------------------------------------
364
365       else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == -1)
366         GMmoins[TETRA4].push_back(it4);
367
368       else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == 0)
369         GMmoins[TETRA4].push_back(it4);
370
371       else if (S[0] == -1 && S[1] == 0 && S[2] == -1 && S[3] == 1)
372         { // Cas 2, arêtes 2 5
373           V[0] = -1;
374           V[1] = -1;
375           V[2] = intersectionSegmentPlan(it4, 2);
376           V[3] = -1;
377           V[4] = -1;
378           V[5] = intersectionSegmentPlan(it4, 5);
379           cas2(V, it4);
380         }
381
382       // -------------------------------------------------------------------
383
384       else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == -1)
385         GMmoins[TETRA4].push_back(it4);
386
387       else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == 0)
388         GMmoins[TETRA4].push_back(it4);
389
390       else if (S[0] == -1 && S[1] == 0 && S[2] == 0 && S[3] == 1)
391         { // Cas 1, arête 2
392           V[0] = -1;
393           V[1] = -1;
394           V[2] = intersectionSegmentPlan(it4, 2);
395           V[3] = -1;
396           V[4] = -1;
397           V[5] = -1;
398           cas1(V, it4);
399         }
400
401       // -------------------------------------------------------------------
402
403       else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == -1)
404         { // Cas 2, arêtes 1 5
405           V[0] = -1;
406           V[1] = intersectionSegmentPlan(it4, 1);
407           V[2] = -1;
408           V[3] = -1;
409           V[4] = -1;
410           V[5] = intersectionSegmentPlan(it4, 5);
411           cas2(V, it4);
412         }
413
414       else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == 0)
415         { // Cas 1, arête 1
416           V[0] = -1;
417           V[1] = intersectionSegmentPlan(it4, 1);
418           V[2] = -1;
419           V[3] = -1;
420           V[4] = -1;
421           V[5] = -1;
422           cas1(V, it4);
423         }
424
425       else if (S[0] == -1 && S[1] == 0 && S[2] == 1 && S[3] == 1)
426         { // Cas 2, arêtes 1 2
427           V[0] = -1;
428           V[1] = intersectionSegmentPlan(it4, 1);
429           V[2] = intersectionSegmentPlan(it4, 2);
430           V[3] = -1;
431           V[4] = -1;
432           V[5] = -1;
433           cas2(V, it4);
434         }
435
436       // -------------------------------------------------------------------
437
438       else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == -1)
439         { // Cas 3, arêtes 0 3 4
440           V[0] = intersectionSegmentPlan(it4, 0);
441           V[1] = -1;
442           V[2] = -1;
443           V[3] = intersectionSegmentPlan(it4, 3);
444           V[4] = intersectionSegmentPlan(it4, 4);
445           V[5] = -1;
446           cas3(V, it4);
447         }
448
449       else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == 0)
450         { // Cas 2, arêtes 0 3
451           V[0] = intersectionSegmentPlan(it4, 0);
452           V[1] = -1;
453           V[2] = -1;
454           V[3] = intersectionSegmentPlan(it4, 3);
455           V[4] = -1;
456           V[5] = -1;
457           cas2(V, it4);
458         }
459
460       else if (S[0] == -1 && S[1] == 1 && S[2] == -1 && S[3] == 1)
461         { // Cas 4, arêtes 0 2 3 5
462           V[0] = intersectionSegmentPlan(it4, 0);
463           V[1] = -1;
464           V[2] = intersectionSegmentPlan(it4, 2);
465           V[3] = intersectionSegmentPlan(it4, 3);
466           V[4] = -1;
467           V[5] = intersectionSegmentPlan(it4, 5);
468           cas4(V, it4);
469         }
470
471       // -------------------------------------------------------------------
472
473       else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == -1)
474         { // Cas 2, arêtes 0 4
475           V[0] = intersectionSegmentPlan(it4, 0);
476           V[1] = -1;
477           V[2] = -1;
478           V[3] = -1;
479           V[4] = intersectionSegmentPlan(it4, 4);
480           V[5] = -1;
481           cas2(V, it4);
482         }
483
484       else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == 0)
485         { // Cas 1, arête 0
486           V[0] = intersectionSegmentPlan(it4, 0);
487           V[1] = -1;
488           V[2] = -1;
489           V[3] = -1;
490           V[4] = -1;
491           V[5] = -1;
492           cas1(V, it4);
493         }
494
495       else if (S[0] == -1 && S[1] == 1 && S[2] == 0 && S[3] == 1)
496         { // Cas 2, arêtes 0 2
497           V[0] = intersectionSegmentPlan(it4, 0);
498           V[1] = -1;
499           V[2] = intersectionSegmentPlan(it4, 2);
500           V[3] = -1;
501           V[4] = -1;
502           V[5] = -1;
503           cas2(V, it4);
504         }
505
506       // -------------------------------------------------------------------
507
508       else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == -1)
509         { // Cas 4, arêtes 0 1 4 5
510           V[0] = intersectionSegmentPlan(it4, 0);
511           V[1] = intersectionSegmentPlan(it4, 1);
512           V[2] = -1;
513           V[3] = -1;
514           V[4] = intersectionSegmentPlan(it4, 4);
515           V[5] = intersectionSegmentPlan(it4, 5);
516           cas4(V, it4);
517         }
518
519       else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == 0)
520         { // Cas 2, arêtes 0 1
521           V[0] = intersectionSegmentPlan(it4, 0);
522           V[1] = intersectionSegmentPlan(it4, 1);
523           V[2] = -1;
524           V[3] = -1;
525           V[4] = -1;
526           V[5] = -1;
527           cas2(V, it4);
528         }
529
530       else if (S[0] == -1 && S[1] == 1 && S[2] == 1 && S[3] == 1)
531         { // Cas 3, arêtes 0 1 2
532           V[0] = intersectionSegmentPlan(it4, 0);
533           V[1] = intersectionSegmentPlan(it4, 1);
534           V[2] = intersectionSegmentPlan(it4, 2);
535           V[3] = -1;
536           V[4] = -1;
537           V[5] = -1;
538           cas3(V, it4);
539         }
540
541       // -------------------------------------------------------------------
542
543       else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == -1)
544         GMmoins[TETRA4].push_back(it4);
545
546       else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == 0)
547         GMmoins[TETRA4].push_back(it4);
548
549       else if (S[0] == 0 && S[1] == -1 && S[2] == -1 && S[3] == 1)
550         { // Cas 2, arêtes 4 5
551           V[0] = -1;
552           V[1] = -1;
553           V[2] = -1;
554           V[3] = -1;
555           V[4] = intersectionSegmentPlan(it4, 4);
556           V[5] = intersectionSegmentPlan(it4, 5);
557           cas2(V, it4);
558         }
559
560       // -------------------------------------------------------------------
561
562       else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == -1)
563         GMmoins[TETRA4].push_back(it4);
564
565       else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == 0)
566         GMmoins[TETRA4].push_back(it4);
567
568       else if (S[0] == 0 && S[1] == -1 && S[2] == 0 && S[3] == 1)
569         { // Cas 1, arête 4
570           V[0] = -1;
571           V[1] = -1;
572           V[2] = -1;
573           V[3] = -1;
574           V[4] = intersectionSegmentPlan(it4, 4);
575           V[5] = -1;
576           cas1(V, it4);
577         }
578
579       // -------------------------------------------------------------------
580
581       else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == -1)
582         { // Cas 2, arêtes 3 5
583           V[0] = -1;
584           V[1] = -1;
585           V[2] = -1;
586           V[3] = intersectionSegmentPlan(it4, 3);
587           V[4] = -1;
588           V[5] = intersectionSegmentPlan(it4, 5);
589           cas2(V, it4);
590         }
591
592       else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == 0)
593         { // Cas 1, arête 3
594           V[0] = -1;
595           V[1] = -1;
596           V[2] = -1;
597           V[3] = intersectionSegmentPlan(it4, 3);
598           V[4] = -1;
599           V[5] = -1;
600           cas1(V, it4);
601         }
602
603       else if (S[0] == 0 && S[1] == -1 && S[2] == 1 && S[3] == 1)
604         { // Cas 2, arêtes 3 4
605           V[0] = -1;
606           V[1] = -1;
607           V[2] = -1;
608           V[3] = intersectionSegmentPlan(it4, 3);
609           V[4] = intersectionSegmentPlan(it4, 4);
610           V[5] = -1;
611           cas2(V, it4);
612         }
613
614       // -------------------------------------------------------------------
615
616       else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == -1)
617         GMmoins[TETRA4].push_back(it4);
618
619       else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == 0)
620         GMmoins[TETRA4].push_back(it4);
621
622       else if (S[0] == 0 && S[1] == 0 && S[2] == -1 && S[3] == 1)
623         { // Cas 1, arête 5
624           V[0] = -1;
625           V[1] = -1;
626           V[2] = -1;
627           V[3] = -1;
628           V[4] = -1;
629           V[5] = intersectionSegmentPlan(it4, 5);
630           cas1(V, it4);
631         }
632
633       // -------------------------------------------------------------------
634
635       else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == -1)
636         GMmoins[TETRA4].push_back(it4);
637
638       else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == 0)
639         {
640           cout << "WARNING: TETRA4 number " << it4
641               << " entirely in the tolerance zone near the cut plane" << endl;
642           cout << " --> affected to group " << str_id_GMmoins << endl;
643           GMmoins[TETRA4].push_back(it4);
644         }
645
646       else if (S[0] == 0 && S[1] == 0 && S[2] == 0 && S[3] == 1)
647         GMplus[TETRA4].push_back(it4);
648
649       // -------------------------------------------------------------------
650
651       else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == -1)
652         { // Cas 1, arête 5
653           V[0] = -1;
654           V[1] = -1;
655           V[2] = -1;
656           V[3] = -1;
657           V[4] = -1;
658           V[5] = intersectionSegmentPlan(it4, 5);
659           cas1(V, it4);
660         }
661
662       else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == 0)
663         GMplus[TETRA4].push_back(it4);
664
665       else if (S[0] == 0 && S[1] == 0 && S[2] == 1 && S[3] == 1)
666         GMplus[TETRA4].push_back(it4);
667
668       // -------------------------------------------------------------------
669
670       else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == -1)
671         { // Cas 2, arêtes 3 4
672           V[0] = -1;
673           V[1] = -1;
674           V[2] = -1;
675           V[3] = intersectionSegmentPlan(it4, 3);
676           V[4] = intersectionSegmentPlan(it4, 4);
677           V[5] = -1;
678           cas2(V, it4);
679         }
680
681       else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == 0)
682         { // Cas 1, arête 3
683           V[0] = -1;
684           V[1] = -1;
685           V[2] = -1;
686           V[3] = intersectionSegmentPlan(it4, 3);
687           V[4] = -1;
688           V[5] = -1;
689           cas1(V, it4);
690         }
691
692       else if (S[0] == 0 && S[1] == 1 && S[2] == -1 && S[3] == 1)
693         { // Cas 2, arêtes 3 5
694           V[0] = -1;
695           V[1] = -1;
696           V[2] = -1;
697           V[3] = intersectionSegmentPlan(it4, 3);
698           V[4] = -1;
699           V[5] = intersectionSegmentPlan(it4, 5);
700           cas2(V, it4);
701         }
702
703       // -------------------------------------------------------------------
704
705       else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == -1)
706         { // Cas 1, arête 4
707           V[0] = -1;
708           V[1] = -1;
709           V[2] = -1;
710           V[3] = -1;
711           V[4] = intersectionSegmentPlan(it4, 4);
712           V[5] = -1;
713           cas1(V, it4);
714         }
715
716       else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == 0)
717         GMplus[TETRA4].push_back(it4);
718
719       else if (S[0] == 0 && S[1] == 1 && S[2] == 0 && S[3] == 1)
720         GMplus[TETRA4].push_back(it4);
721
722       // -------------------------------------------------------------------
723
724       else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == -1)
725         { // Cas 2, arêtes 4 5
726           V[0] = -1;
727           V[1] = -1;
728           V[2] = -1;
729           V[3] = -1;
730           V[4] = intersectionSegmentPlan(it4, 4);
731           V[5] = intersectionSegmentPlan(it4, 5);
732           cas2(V, it4);
733         }
734
735       else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == 0)
736         GMplus[TETRA4].push_back(it4);
737
738       else if (S[0] == 0 && S[1] == 1 && S[2] == 1 && S[3] == 1)
739         GMplus[TETRA4].push_back(it4);
740
741       // -------------------------------------------------------------------
742
743       else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == -1)
744         { // Cas 3, arêtes 0 1 2
745           V[0] = intersectionSegmentPlan(it4, 0);
746           V[1] = intersectionSegmentPlan(it4, 1);
747           V[2] = intersectionSegmentPlan(it4, 2);
748           V[3] = -1;
749           V[4] = -1;
750           V[5] = -1;
751           cas3(V, it4);
752         }
753
754       else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == 0)
755         { // Cas 2, arêtes 0 1
756           V[0] = intersectionSegmentPlan(it4, 0);
757           V[1] = intersectionSegmentPlan(it4, 1);
758           V[2] = -1;
759           V[3] = -1;
760           V[4] = -1;
761           V[5] = -1;
762           cas2(V, it4);
763         }
764
765       else if (S[0] == 1 && S[1] == -1 && S[2] == -1 && S[3] == 1)
766         { // Cas 4, arêtes 0 1 4 5
767           V[0] = intersectionSegmentPlan(it4, 0);
768           V[1] = intersectionSegmentPlan(it4, 1);
769           V[2] = -1;
770           V[3] = -1;
771           V[4] = intersectionSegmentPlan(it4, 4);
772           V[5] = intersectionSegmentPlan(it4, 5);
773           cas4(V, it4);
774         }
775
776       // -------------------------------------------------------------------
777
778       else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == -1)
779         { // Cas 2, arêtes 0 2
780           V[0] = intersectionSegmentPlan(it4, 0);
781           V[1] = -1;
782           V[2] = intersectionSegmentPlan(it4, 2);
783           V[3] = -1;
784           V[4] = -1;
785           V[5] = -1;
786           cas2(V, it4);
787         }
788
789       else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == 0)
790         { // Cas 1, arête 0
791           V[0] = intersectionSegmentPlan(it4, 0);
792           V[1] = -1;
793           V[2] = -1;
794           V[3] = -1;
795           V[4] = -1;
796           V[5] = -1;
797           cas1(V, it4);
798         }
799
800       else if (S[0] == 1 && S[1] == -1 && S[2] == 0 && S[3] == 1)
801         { // Cas 2, arêtes 0 4
802           V[0] = intersectionSegmentPlan(it4, 0);
803           V[1] = -1;
804           V[2] = -1;
805           V[3] = -1;
806           V[4] = intersectionSegmentPlan(it4, 4);
807           V[5] = -1;
808           cas2(V, it4);
809         }
810
811       // -------------------------------------------------------------------
812
813       else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == -1)
814         { // Cas 4, arêtes 0 2 3 5
815           V[0] = intersectionSegmentPlan(it4, 0);
816           V[1] = -1;
817           V[2] = intersectionSegmentPlan(it4, 2);
818           V[3] = intersectionSegmentPlan(it4, 3);
819           V[4] = -1;
820           V[5] = intersectionSegmentPlan(it4, 5);
821           cas4(V, it4);
822         }
823
824       else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == 0)
825         { // Cas 2, arêtes 0 3
826           V[0] = intersectionSegmentPlan(it4, 0);
827           V[1] = -1;
828           V[2] = -1;
829           V[3] = intersectionSegmentPlan(it4, 3);
830           V[4] = -1;
831           V[5] = -1;
832           cas2(V, it4);
833         }
834
835       else if (S[0] == 1 && S[1] == -1 && S[2] == 1 && S[3] == 1)
836         { // Cas 3, arêtes 0 3 4
837           V[0] = intersectionSegmentPlan(it4, 0);
838           V[1] = -1;
839           V[2] = -1;
840           V[3] = intersectionSegmentPlan(it4, 3);
841           V[4] = intersectionSegmentPlan(it4, 4);
842           V[5] = -1;
843           cas3(V, it4);
844         }
845
846       // -------------------------------------------------------------------
847
848       else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == -1)
849         { // Cas 2, arêtes 1 2
850           V[0] = -1;
851           V[1] = intersectionSegmentPlan(it4, 1);
852           V[2] = intersectionSegmentPlan(it4, 2);
853           V[3] = -1;
854           V[4] = -1;
855           V[5] = -1;
856           cas2(V, it4);
857         }
858
859       else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == 0)
860         { // Cas 1, arête 1
861           V[0] = -1;
862           V[1] = intersectionSegmentPlan(it4, 1);
863           V[2] = -1;
864           V[3] = -1;
865           V[4] = -1;
866           V[5] = -1;
867           cas1(V, it4);
868         }
869
870       else if (S[0] == 1 && S[1] == 0 && S[2] == -1 && S[3] == 1)
871         { // Cas 2, arêtes 1 5
872           V[0] = -1;
873           V[1] = intersectionSegmentPlan(it4, 1);
874           V[2] = -1;
875           V[3] = -1;
876           V[4] = -1;
877           V[5] = intersectionSegmentPlan(it4, 5);
878           cas2(V, it4);
879         }
880
881       // -------------------------------------------------------------------
882
883       else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == -1)
884         { // Cas 1, arête 2
885           V[0] = -1;
886           V[1] = -1;
887           V[2] = intersectionSegmentPlan(it4, 2);
888           V[3] = -1;
889           V[4] = -1;
890           V[5] = -1;
891           cas1(V, it4);
892         }
893
894       else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == 0)
895         GMplus[TETRA4].push_back(it4);
896
897       else if (S[0] == 1 && S[1] == 0 && S[2] == 0 && S[3] == 1)
898         GMplus[TETRA4].push_back(it4);
899
900       // -------------------------------------------------------------------
901
902       else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == -1)
903         { // Cas 2, arêtes 2 5
904           V[0] = -1;
905           V[1] = -1;
906           V[2] = intersectionSegmentPlan(it4, 2);
907           V[3] = -1;
908           V[4] = -1;
909           V[5] = intersectionSegmentPlan(it4, 5);
910           cas2(V, it4);
911         }
912
913       else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == 0)
914         GMplus[TETRA4].push_back(it4);
915
916       else if (S[0] == 1 && S[1] == 0 && S[2] == 1 && S[3] == 1)
917         GMplus[TETRA4].push_back(it4);
918
919       // -------------------------------------------------------------------
920
921       else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == -1)
922         { // Cas 4, arêtes 1 2 3 4
923           V[0] = -1;
924           V[1] = intersectionSegmentPlan(it4, 1);
925           V[2] = intersectionSegmentPlan(it4, 2);
926           V[3] = intersectionSegmentPlan(it4, 3);
927           V[4] = intersectionSegmentPlan(it4, 4);
928           V[5] = -1;
929           cas4(V, it4);
930         }
931
932       else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == 0)
933         { // Cas 2, arêtes 1 3
934           V[0] = -1;
935           V[1] = intersectionSegmentPlan(it4, 1);
936           V[2] = -1;
937           V[3] = intersectionSegmentPlan(it4, 3);
938           V[4] = -1;
939           V[5] = -1;
940           cas2(V, it4);
941         }
942
943       else if (S[0] == 1 && S[1] == 1 && S[2] == -1 && S[3] == 1)
944         { // Cas 3, arêtes 1 3 5
945           V[0] = -1;
946           V[1] = intersectionSegmentPlan(it4, 1);
947           V[2] = -1;
948           V[3] = intersectionSegmentPlan(it4, 3);
949           V[4] = -1;
950           V[5] = intersectionSegmentPlan(it4, 5);
951           cas3(V, it4);
952         }
953
954       // -------------------------------------------------------------------
955
956       else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == -1)
957         { // Cas 2, arêtes 2 4
958           V[0] = -1;
959           V[1] = -1;
960           V[2] = intersectionSegmentPlan(it4, 2);
961           V[3] = -1;
962           V[4] = intersectionSegmentPlan(it4, 4);
963           V[5] = -1;
964           cas2(V, it4);
965         }
966
967       else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == 0)
968         GMplus[TETRA4].push_back(it4);
969
970       else if (S[0] == 1 && S[1] == 1 && S[2] == 0 && S[3] == 1)
971         GMplus[TETRA4].push_back(it4);
972
973       // -------------------------------------------------------------------
974
975       else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == -1)
976         { // Cas 3, arêtes 2 4 5
977           V[0] = -1;
978           V[1] = -1;
979           V[2] = intersectionSegmentPlan(it4, 2);
980           V[3] = -1;
981           V[4] = intersectionSegmentPlan(it4, 4);
982           V[5] = intersectionSegmentPlan(it4, 5);
983           cas3(V, it4);
984         }
985
986       else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == 0)
987         GMplus[TETRA4].push_back(it4);
988
989       else if (S[0] == 1 && S[1] == 1 && S[2] == 1 && S[3] == 1)
990         GMplus[TETRA4].push_back(it4);
991
992       else
993         ERREUR("Case not taken into account");
994
995     }
996   cout << chrono() << " - End of iteration on tetra4" << endl;
997
998   // cout << "indexNouveauxNoeuds = " << indexNouveauxNoeuds << endl;
999   newXX.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
1000   newYY.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
1001   newZZ.resize(indexNouveauxNoeuds - MAILLAGE1->nombreNoeudsMaillage);
1002
1003   if (cptNouvellesMailles[TETRA4])
1004     newCNX[TETRA4].resize(4 * cptNouvellesMailles[TETRA4]);
1005   if (cptNouvellesMailles[PYRAM5])
1006     newCNX[PYRAM5].resize(5 * cptNouvellesMailles[PYRAM5]);
1007   if (cptNouvellesMailles[PENTA6])
1008     newCNX[PENTA6].resize(6 * cptNouvellesMailles[PENTA6]);
1009
1010   // =========================================================================================
1011   //                          2. Constitution du maillage final
1012   // =========================================================================================
1013
1014   cout << chrono() << " - Constitution of final mesh" << endl;
1015
1016   MAILLAGE2 = new Maillage(str_id_maillagenew);
1017   MAILLAGE2->dimensionMaillage = MAILLAGE1->dimensionMaillage;
1018   MAILLAGE2->dimensionEspace = MAILLAGE1->dimensionEspace;
1019   strcpy(MAILLAGE2->axisname, MAILLAGE1->axisname);
1020   strcpy(MAILLAGE2->unitname, MAILLAGE1->unitname);
1021   MAILLAGE2->nombreNoeudsMaillage = indexNouveauxNoeuds;
1022   MAILLAGE2->nombreMaillesMaillage = MAILLAGE1->nombreMaillesMaillage + cptNouvellesMailles[TETRA4]
1023       + cptNouvellesMailles[PYRAM5] + cptNouvellesMailles[PENTA6];
1024
1025   // ---------- Coordonnées
1026   // Optimisation de la mémoire au détriment du temps
1027
1028   // Héritage des coordonnées MAILLAGE1
1029   MAILLAGE2->XX = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1030   for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1031     *(MAILLAGE2->XX + i) = *(MAILLAGE1->XX + i);
1032   free(MAILLAGE1->XX);
1033   MAILLAGE2->YY = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1034   for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1035     *(MAILLAGE2->YY + i) = *(MAILLAGE1->YY + i);
1036   free(MAILLAGE1->YY);
1037   MAILLAGE2->ZZ = (float*) malloc(sizeof(float) * MAILLAGE2->nombreNoeudsMaillage);
1038   for (int i = 0; i < MAILLAGE1->nombreNoeudsMaillage; i++)
1039     *(MAILLAGE2->ZZ + i) = *(MAILLAGE1->ZZ + i);
1040   free(MAILLAGE1->ZZ);
1041
1042   // Coordonnées des noeuds créés
1043   for (int i = 0; i < MAILLAGE2->nombreNoeudsMaillage - MAILLAGE1->nombreNoeudsMaillage; i++)
1044     {
1045       *(MAILLAGE2->XX + MAILLAGE1->nombreNoeudsMaillage + i) = newXX[i];
1046       *(MAILLAGE2->YY + MAILLAGE1->nombreNoeudsMaillage + i) = newYY[i];
1047       *(MAILLAGE2->ZZ + MAILLAGE1->nombreNoeudsMaillage + i) = newZZ[i];
1048       // cout << "Nouveaux noeuds, indice " << i << " : " << newXX[i] << " " << newYY[i] << " " << newZZ[i] << " " << endl;
1049     }
1050
1051   // Legacy mailles maillage 1 (volumes seulement)
1052   for (int itm = (int) TETRA4; itm <= (int) HEXA20; itm++)
1053     {
1054       TYPE_MAILLE tm = (TYPE_MAILLE) itm;
1055       if (tm != TETRA4 && tm != PYRAM5 && tm != PENTA6)
1056         {
1057           // Pour les types autres que TETRA4 PYRAM5 PENTA6 on fait seulement pointer CNX2 vers CNX1
1058           if (MAILLAGE1->EFFECTIFS_TYPES[tm])
1059             MAILLAGE2->CNX[tm] = MAILLAGE1->CNX[tm];
1060           MAILLAGE2->EFFECTIFS_TYPES[tm] = MAILLAGE1->EFFECTIFS_TYPES[tm];
1061         }
1062       else
1063         {
1064           // Pour les types TETRA4 PYRAM5 PENTA6 on recopie CNX1 et on ajoute à la suite les newCNX
1065           // cout << "Legacy " << tm << " effectif " << MAILLAGE1->EFFECTIFS_TYPES[tm] << endl;
1066           int tailleType = Nnoeuds(tm);
1067
1068           MAILLAGE2->CNX[tm] = (int*) malloc(sizeof(int) * tailleType * (MAILLAGE1->EFFECTIFS_TYPES[tm]
1069               + cptNouvellesMailles[tm]));
1070           for (int i = 0; i < MAILLAGE1->EFFECTIFS_TYPES[tm]; i++)
1071             for (int j = 0; j < tailleType; j++)
1072               *(MAILLAGE2->CNX[tm] + tailleType * i + j) = *(MAILLAGE1->CNX[tm] + tailleType * i + j);
1073
1074           for (int i = 0; i < cptNouvellesMailles[tm]; i++)
1075             for (int j = 0; j < tailleType; j++)
1076               *(MAILLAGE2->CNX[tm] + tailleType * (MAILLAGE1->EFFECTIFS_TYPES[tm] + i) + j) = newCNX[tm][i * tailleType
1077                   + j];
1078
1079           MAILLAGE2->EFFECTIFS_TYPES[tm] = MAILLAGE1->EFFECTIFS_TYPES[tm] + cptNouvellesMailles[tm];
1080         }
1081     }
1082
1083   // Restit CNX
1084
1085   //   cout << "Maillage 2 - CNX TETRA4 : " << endl;
1086   //  ;
1087   //  for (int i = 0; i < MAILLAGE2->EFFECTIFS_TYPES[TETRA4]; i++)
1088   //    {
1089   //      cout << "Maille " << i << " : ";
1090   //      for (int j = 0; j < 4; j++)
1091   //        cout << MAILLAGE2->CNX[TETRA4][i * 4 + j] << " ";
1092   //      cout << endl;
1093   //    }
1094   //  cout << endl;
1095   //  cout << "Maillage 2 - CNX PENTA6 : " << endl;
1096   //  ;
1097   //  for (int i = 0; i < MAILLAGE2->EFFECTIFS_TYPES[PENTA6]; i++)
1098   //    {
1099   //      cout << "Maille " << i << " : ";
1100   //      for (int j = 0; j < 6; j++)
1101   //        cout << MAILLAGE2->CNX[PENTA6][i * 6 + j] << " ";
1102   //      cout << endl;
1103   //    }
1104   //  cout << endl;
1105
1106   // Groupes de mailles
1107   // MAILLAGE2->GM = MAILLAGE1->GM;
1108   MAILLAGE2->GN.clear();
1109   MAILLAGE2->GM.clear();
1110   MAILLAGE2->GM[str_id_GMplus] = GMplus;
1111   MAILLAGE2->GM[str_id_GMmoins] = GMmoins;
1112
1113   // MAILLAGE2->GN = MAILLAGE1->GN;
1114
1115   MAILLAGE2->eliminationMailles(TETRA4, cutTetras);
1116
1117   cout << chrono() << " - MED file writing" << endl;
1118
1119   MAILLAGE2->outputMED(ficMEDout);
1120   cout << chrono() << " - Finished!" << endl << endl;
1121
1122   return 0;
1123
1124 }