src/Tools/MeshCut/MeshCut_Utils.cxx

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