src/Tools/MeshCut/MeshCut_DC.cxx

   1 // Copyright (C) 2006-2016  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, or (at your option) any later version.
   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 = 0;
  76   float yNormal = 0;
  77   float zNormal = 0;
  78   float xm = 0;
  79   float ym = 0;
  80   float zm = 0;
  81   float tolerance = 0;
  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 }