src/Tools/MeshCut/MeshCut_Utils.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 #include "MeshCut_Utils.hxx"
  21
  22 #include <iostream>
  23 #include <string>
  24 #include <sstream>
  25 #include <cstdlib>
  26 #include <cstring>
  27 #include <ctime>
  28
  29 using namespace std;
  30 using namespace MESHCUT;
  31
  32 //string pathUsers = (string) "/var/www/XMeshLab/users/";
  33
  34 bool MESHCUT::estUnTypeMaille(std::string S)
  35 {
  36   if (S == (string) "POI1" || S == (string) "SEG2" || S == (string) "SEG3" || S == (string) "TRIA3" || S
  37       == (string) "TRIA6" || S == (string) "QUAD4" || S == (string) "QUAD8" || S == (string) "QUAD9" || S
  38       == (string) "TETRA4" || S == (string) "TETRA10" || S == (string) "PYRAM5" || S == (string) "PYRAM13" || S
  39       == (string) "PENTA6" || S == (string) "PENTA15" || S == (string) "HEXA8" || S == (string) "HEXA20" || S
  40       == (string) "HEXA27")
  41     return true;
  42   else
  43     return false;
  44 }
  45
  46 void MESHCUT::ERREUR(const char* msg)
  47 {
  48   cout << endl << "====== ERROR ====== " << msg << endl << endl;
  49   exit(-1);
  50 }
  51
  52 char* MESHCUT::string2char(std::string str)
  53 {
  54   // créer le buffer pour copier la chaîne
  55   size_t size = str.size() + 1;
  56   char* buffer = new char[size];
  57   // copier la chaîne
  58   strncpy(buffer, str.c_str(), size);
  59
  60   // libérer la mémoire
  61   //delete [] buffer;
  62
  63   return buffer;
  64 }
  65
  66 std::string MESHCUT::int2string(int k)
  67 {
  68   std::stringstream oss;
  69   oss << k;
  70   return oss.str(); //  oss.seekp (ios_base::beg);
  71 }
  72
  73 float MESHCUT::char2float(const char* ch)
  74 {
  75   return atof(ch);
  76 }
  77
  78 std::string MESHCUT::float2string(float f)
  79 {
  80   stringstream buf;
  81   buf << fixed << f;
  82   string s = buf.str();
  83   return s;
  84 }
  85
  86 bool MESHCUT::appartient(std::string e, std::string tableau[], int taille)
  87 {
  88   for (int i = 0; i < taille; i++)
  89     if (tableau[i] == e)
  90       return true;
  91   return false;
  92 }
  93
  94 float MESHCUT::arrondi(float x)
  95 {
  96   if (x > 0 && x < 1.0e-5)
  97     return 0;
  98   else if (x < 0 && x > -1.0e-5)
  99     return 0;
 100   else
 101     return x;
 102 }
 103
 104 int MESHCUT::numNoeudPointe(std::string b1, std::string b2, std::string b3)
 105 {
 106   if (b1 == "1" && b2 == "1" && b3 == "2")
 107     return 2;
 108   else if (b1 == "1" && b2 == "2" && b3 == "1")
 109     return 1;
 110   else if (b1 == "1" && b2 == "2" && b3 == "2")
 111     return 2;
 112   else if (b1 == "2" && b2 == "1" && b3 == "1")
 113     return 3;
 114   else if (b1 == "2" && b2 == "1" && b3 == "2")
 115     return 3;
 116   else if (b1 == "2" && b2 == "2" && b3 == "1")
 117     return 1;
 118   else
 119     return -1;
 120 }
 121
 122 std::string MESHCUT::strip(std::string S)
 123 {
 124   if (S.empty())
 125     return S;
 126   int startIndex = S.find_first_not_of(" ");
 127   int endIndex = S.find_last_not_of(" ");
 128   return S.substr(startIndex, (endIndex - startIndex + 1));
 129 }
 130
 131 std::string MESHCUT::entierSur10_g(int i)
 132 {
 133   if (i > 999999999)
 134     ERREUR("trying to write a number superior to 999999999 on more than 10 chars");
 135   if (i < 10)
 136     return int2string(i) + (string) "         ";
 137   else if (i < 100)
 138     return int2string(i) + (string) "        ";
 139   else if (i < 1000)
 140     return int2string(i) + (string) "       ";
 141   else if (i < 10000)
 142     return int2string(i) + (string) "      ";
 143   else if (i < 100000)
 144     return int2string(i) + (string) "     ";
 145   else if (i < 1000000)
 146     return int2string(i) + (string) "    ";
 147   else if (i < 10000000)
 148     return int2string(i) + (string) "   ";
 149   else if (i < 100000000)
 150     return int2string(i) + (string) "  ";
 151   else if (i < 1000000000)
 152     return int2string(i) + (string) " ";
 153   else
 154     return int2string(i);
 155 }
 156
 157 std::string MESHCUT::entierSur10_d(int i)
 158 {
 159   if (i > 999999999)
 160     ERREUR("trying to write a number superior to 999999999 on more than 10 chars");
 161   if (i < 10)
 162     return (string) "         " + int2string(i);
 163   else if (i < 100)
 164     return (string) "        " + int2string(i);
 165   else if (i < 1000)
 166     return (string) "       " + int2string(i);
 167   else if (i < 10000)
 168     return (string) "      " + int2string(i);
 169   else if (i < 100000)
 170     return (string) "     " + int2string(i);
 171   else if (i < 1000000)
 172     return (string) "    " + int2string(i);
 173   else if (i < 10000000)
 174     return (string) "   " + int2string(i);
 175   else if (i < 100000000)
 176     return (string) "  " + int2string(i);
 177   else if (i < 1000000000)
 178     return (string) " " + int2string(i);
 179   else
 180     return int2string(i);
 181 }
 182
 183 std::string MESHCUT::typeEnsight(std::string type)
 184 {
 185   if (type == (string) "POI1")
 186     return (string) "point";
 187   else if (type == (string) "SEG2")
 188     return (string) "bar2";
 189   else if (type == (string) "SEG3")
 190     return (string) "bar2";// ATTENTION, triche!
 191   else if (type == (string) "TRIA3")
 192     return (string) "tria3";
 193   else if (type == (string) "TRIA6")
 194     return (string) "tria3";// ATTENTION, triche!
 195   else if (type == (string) "QUAD4")
 196     return (string) "quad4";
 197   else if (type == (string) "QUAD8")
 198     return (string) "quad4"; // ATTENTION, triche!
 199   else if (type == (string) "QUAD9")
 200     ERREUR("Type QUAD9 not supported by Ensight");
 201   else if (type == (string) "TETRA4")
 202     return (string) "tetra4";
 203   else if (type == (string) "TETRA10")
 204     return (string) "tetra4"; // ATTENTION, triche!
 205   else if (type == (string) "PYRAM5")
 206     return (string) "pyramid5";
 207   else if (type == (string) "PYRAM13")
 208     return (string) "pyramid5"; // ATTENTION, triche!
 209   else if (type == (string) "PENTA6")
 210     return (string) "penta6";
 211   else if (type == (string) "PENTA15")
 212     return (string) "penta6"; // ATTENTION, triche!
 213   else if (type == (string) "HEXA8")
 214     return (string) "hexa8";
 215   else if (type == (string) "HEXA20")
 216     return (string) "hexa8"; // ATTENTION, triche!
 217   else if (type == (string) "HEXA27")
 218     ERREUR("Type HEXA27 not supported by Ensight");
 219   else
 220     ERREUR("Type of element not accepted (method \"typeEnsight\"");
 221   return (string) "";
 222 }
 223
 224 int MESHCUT::Nnoeuds(TYPE_MAILLE type)
 225 {
 226   switch (type)
 227   {
 228     case POI1:
 229       {
 230         return 1;
 231         break;
 232       }
 233     case SEG2:
 234       {
 235         return 2;
 236         break;
 237       }
 238     case SEG3:
 239       {
 240         return 3;
 241         break;
 242       }
 243     case TRIA3:
 244       {
 245         return 3;
 246         break;
 247       }
 248     case TRIA6:
 249       {
 250         return 6;
 251         break;
 252       }
 253     case QUAD4:
 254       {
 255         return 4;
 256         break;
 257       }
 258     case QUAD8:
 259       {
 260         return 8;
 261         break;
 262       }
 263       //case QUAD9:                   { return 9; break; }
 264     case TETRA4:
 265       {
 266         return 4;
 267         break;
 268       }
 269     case TETRA10:
 270       {
 271         return 10;
 272         break;
 273       }
 274     case PYRAM5:
 275       {
 276         return 5;
 277         break;
 278       }
 279     case PYRAM13:
 280       {
 281         return 13;
 282         break;
 283       }
 284     case PENTA6:
 285       {
 286         return 6;
 287         break;
 288       }
 289     case PENTA15:
 290       {
 291         return 15;
 292         break;
 293       }
 294     case HEXA8:
 295       {
 296         return 8;
 297         break;
 298       }
 299     case HEXA20:
 300       {
 301         return 20;
 302         break;
 303       }
 304       //case HEXA27:                      { return 27;    break; }
 305     default:
 306       ERREUR("Type of elem not accepted (method Nnoeuds)");
 307   }
 308   return 0;
 309 }
 310
 311 int MESHCUT::NnoeudsGeom(TYPE_MAILLE type)
 312 {
 313   switch (type)
 314   {
 315     case POI1:
 316       {
 317         return 1;
 318         break;
 319       }
 320     case SEG2:
 321       {
 322         return 2;
 323         break;
 324       }
 325     case SEG3:
 326       {
 327         return 2;
 328         break;
 329       }
 330     case TRIA3:
 331       {
 332         return 3;
 333         break;
 334       }
 335     case TRIA6:
 336       {
 337         return 3;
 338         break;
 339       }
 340     case QUAD4:
 341       {
 342         return 4;
 343         break;
 344       }
 345     case QUAD8:
 346       {
 347         return 4;
 348         break;
 349       }
 350       //case QUAD9:                   { return 9; break; }
 351     case TETRA4:
 352       {
 353         return 4;
 354         break;
 355       }
 356     case TETRA10:
 357       {
 358         return 4;
 359         break;
 360       }
 361     case PYRAM5:
 362       {
 363         return 5;
 364         break;
 365       }
 366     case PYRAM13:
 367       {
 368         return 5;
 369         break;
 370       }
 371     case PENTA6:
 372       {
 373         return 6;
 374         break;
 375       }
 376     case PENTA15:
 377       {
 378         return 6;
 379         break;
 380       }
 381     case HEXA8:
 382       {
 383         return 8;
 384         break;
 385       }
 386     case HEXA20:
 387       {
 388         return 8;
 389         break;
 390       }
 391       //case HEXA27:                      { return 27;    break; }
 392     default:
 393       ERREUR("Type of elem not accepted (method NnoeudsGeom)");
 394   }
 395   return 0;
 396 }
 397
 398 int MESHCUT::codeGMSH(std::string type)
 399 {
 400   if (type == (string) "POI1")
 401     ERREUR("POI1 not taken into account by GMSH");
 402   else if (type == (string) "SEG2")
 403     return 1;
 404   else if (type == (string) "SEG3")
 405     return 8;
 406   else if (type == (string) "TRIA3")
 407     return 2;
 408   else if (type == (string) "TRIA6")
 409     return 9;
 410   else if (type == (string) "QUAD4")
 411     return 3;
 412   else if (type == (string) "QUAD8")
 413     return 16;
 414   else if (type == (string) "QUAD9")
 415     return 10;
 416   else if (type == (string) "TETRA4")
 417     return 4;
 418   else if (type == (string) "TETRA10")
 419     return 11;
 420   else if (type == (string) "PYRAM5")
 421     return 7;
 422   else if (type == (string) "PENTA6")
 423     return 6;
 424   else if (type == (string) "PENTA15")
 425     return 18;
 426   else if (type == (string) "HEXA8")
 427     return 5;
 428   else if (type == (string) "HEXA20")
 429     return 17;
 430   else if (type == (string) "HEXA27")
 431     return 12;
 432   else
 433     ERREUR("Type of elem not accepted (method codeGMSH)");
 434   return 0;
 435 }
 436
 437 std::string MESHCUT::floatEnsight(float x)
 438 {
 439   char buf[12];
 440   string s;
 441   if (x < 0.0)
 442     sprintf(buf, "%1.5E", x);
 443   else
 444     sprintf(buf, " %1.5E", x);
 445   s = (string) buf;
 446   s.erase(10, 1);
 447   return s;
 448 }
 449
 450 bool MESHCUT::typeComplexe(std::string type)
 451 {
 452   if (type == (string) "SEG3")
 453     return true;
 454   else if (type == (string) "TRIA6")
 455     return true;
 456   else if (type == (string) "QUAD8")
 457     return true;
 458   else if (type == (string) "QUAD9")
 459     return true;
 460   else if (type == (string) "TETRA10")
 461     return true;
 462   else if (type == (string) "PYRAM13")
 463     return true;
 464   else if (type == (string) "PENTA15")
 465     return true;
 466   else if (type == (string) "HEXA20")
 467     return true;
 468   else if (type == (string) "HEXA27")
 469     return true;
 470   else
 471     return false;
 472 }
 473
 474 std::string MESHCUT::ASTER8(std::string s)
 475 {
 476   if (s.size() == 0)
 477     return (s + (string) "        ");
 478   else if (s.size() == 1)
 479     return (s + (string) "       ");
 480   else if (s.size() == 2)
 481     return (s + (string) "      ");
 482   else if (s.size() == 3)
 483     return (s + (string) "     ");
 484   else if (s.size() == 4)
 485     return (s + (string) "    ");
 486   else if (s.size() == 5)
 487     return (s + (string) "   ");
 488   else if (s.size() == 6)
 489     return (s + (string) "  ");
 490   else if (s.size() == 7)
 491     return (s + (string) " ");
 492   else if (s.size() == 8)
 493     return (s);
 494   else
 495     ERREUR("More than 8 char for an ASTER string");
 496   return (s);
 497 }
 498
 499 /*!
 500  *  Distance à laquelle doit se tenir l'observateur sur un axe
 501  *  pour voir sous 90° un objet centré de dimensions a et b selon les deux autres axes.
 502  *  Si on ne tient pas compte de la dimension de l'objet selon l'axe choisi,
 503  *  la formule d_obs=max(a,b)/2 donne la cote
 504  *  qui permet de voir l'objet plat dans un angle de 90°.
 505  *  A cela il faut ajouter la dimension de l'objet selon l'axe d'observation = c.
 506  *
 507  *  @param a dimensions de l'objet selon un des axes normal à l'axe d'observation
 508  *  @param b dimensions de l'objet selon l'autre axe normal à l'axe d'observation
 509  *  @param c est la dimension de l'objet selon l'axe d'observation
 510  */
 511 float MESHCUT::dObservateur(float a, float b, float c)
 512 {
 513   return (max(a, b) / 2.0 + c);
 514 }
 515
 516 int MESHCUT::copieFichier(std::string source, std::string cible)
 517 {
 518   FILE *fsource, *fcible;
 519   char buffer[512];
 520   int NbLu;
 521   if ((fsource = fopen(string2char(source), "rb")) == NULL)
 522     return -1;
 523   if ((fcible = fopen(string2char(cible), "wb")) == NULL)
 524     {
 525       fclose(fsource);
 526       return -2;
 527     }
 528   while ((NbLu = fread(buffer, 1, 512, fsource)) != 0)
 529     fwrite(buffer, 1, NbLu, fcible);
 530   fclose(fcible);
 531   fclose(fsource);
 532   return 0;
 533 }
 534
 535 med_geometry_type MESHCUT::InstanceMGE(TYPE_MAILLE TYPE)
 536 {
 537   med_geometry_type typeBanaliseMED;
 538
 539   switch (TYPE)
 540   {
 541     case POI1:
 542       typeBanaliseMED = MED_POINT1;
 543       break; // Attention, piège !
 544     case SEG2:
 545       typeBanaliseMED = MED_SEG2;
 546       break;
 547     case SEG3:
 548       typeBanaliseMED = MED_SEG3;
 549       break;
 550     case TRIA3:
 551       typeBanaliseMED = MED_TRIA3;
 552       break;
 553     case TRIA6:
 554       typeBanaliseMED = MED_TRIA6;
 555       break;
 556     case QUAD4:
 557       typeBanaliseMED = MED_QUAD4;
 558       break;
 559     case QUAD8:
 560       typeBanaliseMED = MED_QUAD8;
 561       break;
 562     case TETRA4:
 563       typeBanaliseMED = MED_TETRA4;
 564       break;
 565     case TETRA10:
 566       typeBanaliseMED = MED_TETRA10;
 567       break;
 568     case PYRAM5:
 569       typeBanaliseMED = MED_PYRA5;
 570       break; // Attention, piège !
 571     case PYRAM13:
 572       typeBanaliseMED = MED_PYRA13;
 573       break; // Attention, piège !
 574     case PENTA6:
 575       typeBanaliseMED = MED_PENTA6;
 576       break;
 577     case PENTA15:
 578       typeBanaliseMED = MED_PENTA15;
 579       break;
 580     case HEXA8:
 581       typeBanaliseMED = MED_HEXA8;
 582       break;
 583     case HEXA20:
 584       typeBanaliseMED = MED_HEXA20;
 585       break;
 586     default:
 587       ERREUR("Method InstanceMGE, unknown type ");
 588   }
 589   return typeBanaliseMED;
 590 }
 591
 592 int MESHCUT::chrono()
 593 {
 594   return clock() / CLOCKS_PER_SEC;
 595 }
 596
 597 TYPE_MAILLE MESHCUT::typeMaille(std::string type)
 598 {
 599   if (type == (string) "POI1")
 600     return POI1;
 601   else if (type == (string) "SEG2")
 602     return SEG2;
 603   else if (type == (string) "SEG3")
 604     return SEG3;
 605   else if (type == (string) "TRIA3")
 606     return TRIA3;
 607   else if (type == (string) "TRIA6")
 608     return TRIA6;
 609   else if (type == (string) "QUAD4")
 610     return QUAD4;
 611   else if (type == (string) "QUAD8")
 612     return QUAD8;
 613   else if (type == (string) "TETRA4")
 614     return TETRA4;
 615   else if (type == (string) "TETRA10")
 616     return TETRA10;
 617   else if (type == (string) "PYRAM5")
 618     return PYRAM5;
 619   else if (type == (string) "PYRAM13")
 620     return PYRAM13;
 621   else if (type == (string) "PENTA6")
 622     return PENTA6;
 623   else if (type == (string) "PENTA15")
 624     return PENTA15;
 625   else if (type == (string) "HEXA8")
 626     return HEXA8;
 627   else if (type == (string) "HEXA20")
 628     return HEXA20;
 629   else
 630     ERREUR("ERROR method typeMaille, unknown type");
 631   return POI1;
 632 }
 633
 634 std::string MESHCUT::MGE2string(med_geometry_type MGE)
 635 {
 636   if (MGE == MED_NONE)
 637     return (string) "NOEUD";
 638   else if (MGE == MED_POINT1)
 639     return (string) "POI1";
 640   else if (MGE == MED_SEG2)
 641     return (string) "SEG2";
 642   else if (MGE == MED_SEG3)
 643     return (string) "SEG3";
 644   else if (MGE == MED_TRIA3)
 645     return (string) "TRIA3";
 646   else if (MGE == MED_TRIA6)
 647     return (string) "TRIA6";
 648   else if (MGE == MED_QUAD4)
 649     return (string) "QUAD4";
 650   else if (MGE == MED_QUAD8)
 651     return (string) "QUAD8";
 652   else if (MGE == MED_TETRA4)
 653     return (string) "TETRA4";
 654   else if (MGE == MED_TETRA10)
 655     return (string) "TETRA10";
 656   else if (MGE == MED_PYRA5)
 657     return (string) "PYRAM5";
 658   else if (MGE == MED_PYRA13)
 659     return (string) "PYRAM13";
 660   else if (MGE == MED_PENTA6)
 661     return (string) "PENTA6";
 662   else if (MGE == MED_PENTA15)
 663     return (string) "PENTA15";
 664   else if (MGE == MED_HEXA8)
 665     return (string) "HEXA8";
 666   else if (MGE == MED_HEXA20)
 667     return (string) "HEXA20";
 668   else
 669     ERREUR("ERROR method MGE2string, unknown type");
 670   return (string) "NOEUD";
 671 }
 672
 673 std::string MESHCUT::TM2string(TYPE_MAILLE MGE)
 674 {
 675   if (MGE == POI1)
 676     return (string) "POI1";
 677   else if (MGE == SEG2)
 678     return (string) "SEG2";
 679   else if (MGE == SEG3)
 680     return (string) "SEG3";
 681   else if (MGE == TRIA3)
 682     return (string) "TRIA3";
 683   else if (MGE == TRIA6)
 684     return (string) "TRIA6";
 685   else if (MGE == QUAD4)
 686     return (string) "QUAD4";
 687   else if (MGE == QUAD8)
 688     return (string) "QUAD8";
 689   else if (MGE == TETRA4)
 690     return (string) "TETRA4";
 691   else if (MGE == TETRA10)
 692     return (string) "TETRA10";
 693   else if (MGE == PYRAM5)
 694     return (string) "PYRAM5";
 695   else if (MGE == PYRAM13)
 696     return (string) "PYRAM13";
 697   else if (MGE == PENTA6)
 698     return (string) "PENTA6";
 699   else if (MGE == PENTA15)
 700     return (string) "PENTA15";
 701   else if (MGE == HEXA8)
 702     return (string) "HEXA8";
 703   else if (MGE == HEXA20)
 704     return (string) "HEXA20";
 705   else
 706     ERREUR("ERROR method TM2string, unknown type");
 707   return (string) "POI1";
 708 }
 709
 710 TYPE_MAILLE MESHCUT::string2TM(std::string stm)
 711 {
 712   if (stm == (string) "POI1")
 713     return POI1;
 714   else if (stm == (string) "SEG2")
 715     return SEG2;
 716   else if (stm == (string) "SEG3")
 717     return SEG3;
 718   else if (stm == (string) "TRIA3")
 719     return TRIA3;
 720   else if (stm == (string) "TRIA6")
 721     return TRIA6;
 722   else if (stm == (string) "QUAD4")
 723     return QUAD4;
 724   else if (stm == (string) "QUAD8")
 725     return QUAD8;
 726   else if (stm == (string) "TETRA4")
 727     return TETRA4;
 728   else if (stm == (string) "TETRA10")
 729     return TETRA10;
 730   else if (stm == (string) "PYRAM5")
 731     return PYRAM5;
 732   else if (stm == (string) "PYRAM13")
 733     return PYRAM13;
 734   else if (stm == (string) "PENTA6")
 735     return PENTA6;
 736   else if (stm == (string) "PENTA15")
 737     return PENTA15;
 738   else if (stm == (string) "HEXA8")
 739     return HEXA8;
 740   else if (stm == (string) "HEXA20")
 741     return HEXA20;
 742   else
 743     ERREUR("ERROR method string2TM, unknown type");
 744   return POI1;
 745 }
 746
 747 std::string MESHCUT::coordIndex_ILS(TYPE_MAILLE tm)
 748 {
 749   if (tm == SEG2)
 750     return (string) " 0,1 ";
 751   else if (tm == SEG3)
 752     return (string) " 0,1 "; // Idem SEG2
 753   else if (tm == TRIA3)
 754     return (string) " 0,1,2,0 ";
 755   else if (tm == TRIA6)
 756     return (string) " 0,1,2,0 ";
 757   else if (tm == QUAD4)
 758     return (string) " 0,1,2,3,0 ";
 759   else if (tm == QUAD8)
 760     return (string) " 0,1,2,3,0 ";
 761   else if (tm == TETRA4)
 762     return (string) " 0,1,2,0,-1, 0,3,-1, 1,3,-1, 2,3,-1 ";
 763   else if (tm == TETRA10)
 764     return (string) " 0,1,2,0,-1, 0,3,-1, 1,3,-1, 2,3,-1 ";
 765   else if (tm == PYRAM5)
 766     return (string) " 0,1,2,3,0,-1, 0,4,-1, 1,4,-1, 2,4,-1, 3,4,-1 ";
 767   else if (tm == PYRAM13)
 768     return (string) " 0,1,2,3,0,-1, 0,4,-1, 1,4,-1, 2,4,-1, 3,4,-1 ";
 769   else if (tm == PENTA6)
 770     return (string) " 0,1,2,0,-1, 3,4,5,3,-1, 0,3,-1, 1,4,-1, 2,5,-1 ";
 771   else if (tm == PENTA15)
 772     return (string) " 0,1,2,0,-1, 3,4,5,3,-1, 0,3,-1, 1,4,-1, 2,5,-1 ";
 773   else if (tm == HEXA8)
 774     return (string) " 0,1,2,3,0,-1, 4,5,6,7,4,-1, 0,4,-1, 1,5,-1, 2,6,-1, 3,7,-1 ";
 775   else if (tm == HEXA20)
 776     return (string) " 0,1,2,3,0,-1, 4,5,6,7,4,-1, 0,4,-1, 1,5,-1, 2,6,-1, 3,7,-1 ";
 777   else
 778     return (string) "";
 779 }
 780
 781 std::string MESHCUT::coordIndex_IFS(TYPE_MAILLE tm)
 782 {
 783   if (tm == SEG2)
 784     return (string) "  ";
 785   else if (tm == SEG3)
 786     return (string) "  "; // Idem SEG2
 787   else if (tm == TRIA3)
 788     return (string) " 0,1,2,0,-1, 0,2,1,0,-1 ";
 789   else if (tm == TRIA6)
 790     return (string) " 0,1,2,0,-1, 0,2,1,0,-1 ";
 791   else if (tm == QUAD4)
 792     return (string) " 0,1,2,3,0,-1, 0,3,2,1,0,-1 ";
 793   else if (tm == QUAD8)
 794     return (string) " 0,1,2,3,0,-1, 0,3,2,1,0,-1 ";
 795   else if (tm == TETRA4)
 796     return (string) " 0,1,2,0,-1, 0,2,1,0,-1, 0,3,1,0,-1, 0,1,3,0,-1, 1,3,2,1,-1, 1,2,3,1,-1, 0,2,3,0,-1, 0,3,2,0,-1 ";
 797   else if (tm == TETRA10)
 798     return (string) " 0,1,2,0,-1, 0,2,1,0,-1, 0,3,1,0,-1, 0,1,3,0,-1, 1,3,2,1,-1, 1,2,3,1,-1, 0,2,3,0,-1, 0,3,2,0,-1 ";
 799   else if (tm == PYRAM5)
 800     return (string) " 0,1,2,3,0,-1, 0,3,2,1,0,-1, 0,1,4,0,-1, 0,4,1,0,-1, 1,2,4,1,-1, 1,4,2,1,-1, 2,4,3,2,-1, 2,3,4,2,-1, 3,4,0,3,-1, 3,0,4,3,-1 ";
 801   else if (tm == PYRAM13)
 802     return (string) " 0,1,2,3,0,-1, 0,3,2,1,0,-1, 0,1,4,0,-1, 0,4,1,0,-1, 1,2,4,1,-1, 1,4,2,1,-1, 2,4,3,2,-1, 2,3,4,2,-1, 3,4,0,3,-1, 3,0,4,3,-1 ";
 803   else if (tm == PENTA6)
 804     return (string) " 0,1,2,0,-1, 0,2,1,0,-1, 3,4,5,3,-1, 3,5,4,3,-1, 0,1,4,3,0,-1, 0,3,4,1,0,-1, 1,4,5,2,1,-1, 1,2,5,4,1,-1, 0,3,5,2,0,-1, 0,2,5,3,0,-1 ";
 805   else if (tm == PENTA15)
 806     return (string) " 0,1,2,0,-1, 0,2,1,0,-1, 3,4,5,3,-1, 3,5,4,3,-1, 0,1,4,3,0,-1, 0,3,4,1,0,-1, 1,4,5,2,1,-1, 1,2,5,4,1,-1, 0,3,5,2,0,-1, 0,2,5,3,0,-1 ";
 807   else if (tm == HEXA8)
 808     return (string) " 0,1,2,3,0,-1, 0,3,2,1,0,-1, 1,5,6,2,1,-1, 1,2,6,5,1,-1, 5,4,7,6,5,-1, 5,6,7,4,5,-1, 4,0,3,7,4,-1, 4,7,3,0,4,-1, 0,4,5,1,0,-1, 0,1,5,4,0,-1, 3,7,6,2,3,-1, 3,2,6,7,3,-1 ";
 809   else if (tm == HEXA20)
 810     return (string) " 0,1,2,3,0,-1, 0,3,2,1,0,-1, 1,5,6,2,1,-1, 1,2,6,5,1,-1, 5,4,7,6,5,-1, 5,6,7,4,5,-1, 4,0,3,7,4,-1, 4,7,3,0,4,-1, 0,4,5,1,0,-1, 0,1,5,4,0,-1, 3,7,6,2,3,-1, 3,2,6,7,3,-1 ";
 811   else
 812     return (string) "";
 813 }
 814
 815 std::string MESHCUT::SIGNE(double x)
 816 {
 817   if (x < 0)
 818     return "-";
 819   else if (x > 0)
 820     return "+";
 821   else
 822     return "0";
 823 }
 824
 825 void MESHCUT::champType(std::string type, med_entity_type MEM, med_geometry_type MGE, med_idt fid, med_idt fidout,
 826                         char *maa, char *nomChamp, char *nomChampMoy, med_field_type typeChamp, char *compChamp,
 827                         char *unitChamp, med_int nCompChamp, std::map<std::string, int> REFGAUSS, int ichamp)
 828 {
 829
 830   bool debug = true;
 831   int ipt, nmailles, ngauss, imaille, igauss, icomp;
 832   //  int ival, ngpdt;
 833   med_int nval, numdt, numo, nPasTemps;
 834   char dtunit[MED_SNAME_SIZE + 1] = "";
 835   char locname[MED_NAME_SIZE + 1] = "";
 836   char nomprofil[MED_NAME_SIZE + 1] = "";
 837   med_float dt = 0.0;
 838   med_float *valr = NULL;
 839   med_float *valr2 = NULL;
 840   med_bool local;
 841   //  med_int nbrefmaa;
 842   med_field_type fieldType;
 843
 844   if (MEDfieldInfo(fid, ichamp, nomChamp, maa, &local, &fieldType, compChamp, unitChamp, dtunit, &nPasTemps) < 0)
 845     ERREUR("Error MEDfieldInfo");
 846   cout << type << " : " << (int) nPasTemps << " timestep  " << endl;
 847
 848   for (ipt = 1; ipt <= nPasTemps; ipt++)
 849     {
 850       //for (ipt=1; ipt<=min(nPasTemps,1); ipt++) {
 851       if (debug)
 852         cout << endl;
 853       if (debug)
 854         cout << "************************************************************" << endl;
 855       if (debug)
 856         cout << "                    FIELD " << ichamp << endl;
 857       if (debug)
 858         cout << "          " << nomChamp << endl;
 859       if (debug)
 860         cout << "          " << type << "   ---   Timestep " << ipt << endl;
 861       if (debug)
 862         cout << "************************************************************" << endl;
 863       if (debug)
 864         cout << endl;
 865
 866       if (MEDfieldComputingStepInfo(fid, nomChamp, ipt, &numdt, &numo, &dt) < 0)
 867         {
 868           cout << endl;
 869           cout << endl << "####################################################################" << endl;
 870           cout << "                   ERROR MEDpasdetempsInfo                         " << endl;
 871           cout << endl << "####################################################################" << endl;
 872           cout << "                  Field: " << (string) nomChamp << endl;
 873           cout << "                  Geometrie: " << MGE2string(MGE) << endl;
 874           cout << "                  Timestep " << ipt << " ignored" << endl;
 875
 876           continue;
 877         }
 878
 879       med_int profilesize, nintegrationpoint;
 880       nval = MEDfieldnValueWithProfile(fid, nomChamp, numdt, numo, MEM, MGE, ipt, MED_COMPACT_PFLMODE, nomprofil,
 881                                        &profilesize, locname, &nintegrationpoint);
 882       if (debug)
 883         cout << "     Number of values in this timestep: " << (int) nval << endl;
 884
 885       if (typeChamp == MED_FLOAT64)
 886         valr = (med_float*) calloc(nCompChamp * nval, sizeof(med_float));
 887       else
 888         ERREUR("Type of field not taken into account");
 889
 890       if (MEDfieldValueWithProfileRd(fid, maa, numdt, numo, MEM, MGE, MED_COMPACT_PFLMODE, nomprofil,
 891                                      MED_FULL_INTERLACE, MED_ALL_CONSTITUENT, (unsigned char*) valr) < 0)
 892         {
 893           cout << endl;
 894           cout << endl << "####################################################################" << endl;
 895           cout << "                         ERROR MEDchampLire                        " << endl;
 896           cout << endl << "####################################################################" << endl;
 897           cout << endl;
 898           cout << "   Field: " << (string) nomChamp << endl;
 899           cout << "   Geometry: " << MGE2string(MGE) << endl;
 900           cout << "   Timestep " << ipt << " ignored" << endl;
 901           cout << endl << endl;
 902           continue;
 903         }
 904
 905       if (debug)
 906         cout << "       profile  = " << (string) nomprofil << endl;
 907       // Localisation du champ aux points de Gauss
 908       if (debug)
 909         cout << "       locname = " << (string) locname << endl;
 910
 911       if (REFGAUSS[(string) locname])
 912         {
 913           ngauss = REFGAUSS[(string) locname];
 914           if (debug)
 915             cout << "       " << ngauss << " Gauss points by element)" << endl;
 916         }
 917       else
 918         ngauss = 1;
 919
 920       nmailles = nval / ngauss;
 921       if (debug)
 922         cout << "      Nbre de mailles: " << nmailles << endl;
 923
 924       if (debug)
 925         {
 926           cout << endl << "       Liste des valeurs du champ brut aux 3 premiers éléments:" << endl;
 927           for (imaille = 0; imaille < min(nmailles, 3); imaille++)
 928             {
 929               cout << "         Maille " << imaille << endl;
 930               for (igauss = 0; igauss < ngauss; igauss++)
 931                 {
 932                   cout << "             PG " << igauss << " : ";
 933                   for (icomp = 0; icomp < nCompChamp; icomp++)
 934                     cout << " " << *(valr + imaille * ngauss * nCompChamp + igauss * nCompChamp + icomp);
 935                   cout << endl;
 936                 }
 937               cout << endl;
 938             }
 939           cout << endl;
 940         }
 941
 942       if (ngauss > 1)
 943         {
 944
 945           valr2 = (med_float*) calloc(nCompChamp * nmailles, sizeof(med_float));
 946
 947           if (debug)
 948             cout << endl << "       Moyenne sur les PG des mailles" << endl;
 949           for (imaille = 0; imaille < nmailles; imaille++)
 950             {
 951               for (icomp = 0; icomp < nCompChamp; icomp++)
 952                 {
 953                   float valCompMaille = 0.0;
 954                   for (igauss = 0; igauss < ngauss; igauss++)
 955                     valCompMaille += *(valr + imaille * ngauss * nCompChamp + igauss * nCompChamp + icomp);
 956                   *(valr2 + imaille * nCompChamp + icomp) = valCompMaille / ngauss;
 957
 958                 }
 959             }
 960
 961           //cout << endl << "Nom champ moy = " <<  (string)nomChampMoy << endl;
 962           //cout << endl << "Type champ = " <<  typeChamp << endl;
 963           //cout << endl << "Comp champ = " <<  (string)compChamp << endl;
 964           //cout << endl << "Unit champ = " <<  (string)unitChamp << endl;
 965           //cout << endl << "N comp champ = " <<  nCompChamp << endl;
 966
 967           if (MEDfieldValueWithProfileWr(fidout, nomChampMoy, numdt, numo, dt, MEM, MGE, MED_COMPACT_PFLMODE,
 968                                          nomprofil, MED_NO_LOCALIZATION, MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
 969                                          (med_int) nmailles, (unsigned char*) valr2) < 0)
 970             {
 971               cout << endl;
 972               cout << endl << "********************************************************************" << endl;
 973               cout << "********************                         ***********************" << endl;
 974               cout << "********************   ERROR MEDchampEcr     ***********************" << endl;
 975               cout << "********************                         ***********************" << endl;
 976               cout << "********************************************************************" << endl;
 977               cout << endl;
 978               cout << "   Champ: " << (string) nomChampMoy << endl;
 979               cout << "   Géométrie: " << MGE2string(MGE) << endl;
 980               cout << "   Pas de temps " << ipt << " ignoré" << endl;
 981               cout << endl << endl;
 982               continue;
 983             }
 984
 985           if (debug)
 986             cout << "    Writing mean values in new field: OK " << endl;
 987
 988           // Restitution du champ moyenné
 989           if (debug)
 990             {
 991               cout << endl << "       Liste des valeurs du champ moyenné aux 3 premiers éléments:" << endl;
 992               for (imaille = 0; imaille < min(nmailles, 3); imaille++)
 993                 {
 994                   cout << "         Maille " << imaille << endl;
 995                   for (icomp = 0; icomp < nCompChamp; icomp++)
 996                     cout << " " << *(valr2 + imaille * nCompChamp + icomp);
 997                   cout << endl;
 998                 }
 999               cout << endl;
1000             }
1001
1002         }
1003
1004       free(valr);
1005       free(valr2);
1006
1007     } // boucle sur les pas de temps
1008
1009   cout << endl;
1010 }
1011
1012 std::string MESHCUT::nomMaille(TYPE_MAILLE tm, int nl)
1013 {
1014   return (TM2string(tm) + (string) "_" + int2string(nl));
1015 }
1016
1017 bool MESHCUT::appartientVN(int n, std::vector<int> V)
1018 {
1019   bool app = false;
1020   for (unsigned int i = 0; i < V.size(); i++)
1021     if (n == V[i])
1022       {
1023         app = true;
1024         break;
1025       }
1026   return app;
1027 }
1028
1029 float MESHCUT::distance2(float x1, float y1, float z1, float x2, float y2, float z2)
1030 {
1031   return (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1);
1032 }
1033
1034 /*!
1035  *  Conversion HL-MED d'une table de connectivités
1036  */
1037 void MESHCUT::conversionCNX(int *CNXtm, TYPE_MAILLE tm, int N)
1038 {
1039
1040   int n = Nnoeuds(tm);
1041
1042   if (tm == TETRA4)
1043     {
1044       for (int i = 0; i < N; i++)
1045         {
1046           int i1 = CNXtm[i * n + 1];
1047           int i2 = CNXtm[i * n + 2];
1048           CNXtm[i * n + 1] = i2;
1049           CNXtm[i * n + 2] = i1;
1050         }
1051     }
1052   else if (tm == PYRAM5)
1053     {
1054       for (int i = 0; i < N; i++)
1055         {
1056           int i1 = CNXtm[i * n + 1];
1057           int i3 = CNXtm[i * n + 3];
1058           CNXtm[i * n + 1] = i3;
1059           CNXtm[i * n + 3] = i1;
1060         }
1061     }
1062   else if (tm == PENTA6)
1063     {
1064       for (int i = 0; i < N; i++)
1065         {
1066           int i0 = CNXtm[i * n + 0];
1067           int i1 = CNXtm[i * n + 1];
1068           int i2 = CNXtm[i * n + 2];
1069           int i3 = CNXtm[i * n + 3];
1070           int i4 = CNXtm[i * n + 4];
1071           int i5 = CNXtm[i * n + 5];
1072           CNXtm[i * n + 0] = i3;
1073           CNXtm[i * n + 1] = i4;
1074           CNXtm[i * n + 2] = i5;
1075           CNXtm[i * n + 3] = i0;
1076           CNXtm[i * n + 4] = i1;
1077           CNXtm[i * n + 5] = i2;
1078         }
1079     }
1080
1081   else if (tm == HEXA8)
1082     {
1083       for (int i = 0; i < N; i++)
1084         {
1085           int i0 = CNXtm[i * n + 0];
1086           int i1 = CNXtm[i * n + 1];
1087           int i2 = CNXtm[i * n + 2];
1088           int i3 = CNXtm[i * n + 3];
1089           int i4 = CNXtm[i * n + 4];
1090           int i5 = CNXtm[i * n + 5];
1091           int i6 = CNXtm[i * n + 6];
1092           int i7 = CNXtm[i * n + 7];
1093           CNXtm[i * n + 0] = i4;
1094           CNXtm[i * n + 1] = i5;
1095           CNXtm[i * n + 2] = i6;
1096           CNXtm[i * n + 3] = i7;
1097           CNXtm[i * n + 4] = i0;
1098           CNXtm[i * n + 5] = i1;
1099           CNXtm[i * n + 6] = i2;
1100           CNXtm[i * n + 7] = i3;
1101         }
1102     }
1103 }
1104