src/MULTIPR/MULTIPR_Mesh.cxx

   1 // Project MULTIPR, IOLS WP1.2.1 - EDF/CS
   2 // Partitioning/decimation module for the SALOME v3.2 platform
   3
   4 /**
   5  * \file    MULTIPR_Mesh.cxx
   6  *
   7  * \brief   see MULTIPR_Mesh.hxx
   8  *
   9  * \author  Olivier LE ROUX - CS, Virtual Reality Dpt
  10  *
  11  * \date    01/2007
  12  */
  13
  14 //*****************************************************************************
  15 // Includes section
  16 //*****************************************************************************
  17
  18 #include "MULTIPR_Mesh.hxx"
  19 #include "MULTIPR_Nodes.hxx"
  20 #include "MULTIPR_Elements.hxx"
  21 #include "MULTIPR_Family.hxx"
  22 #include "MULTIPR_Profil.hxx"
  23 #include "MULTIPR_GaussLoc.hxx"
  24 #include "MULTIPR_Field.hxx"
  25 #include "MULTIPR_MeshDis.hxx"
  26 #include "MULTIPR_PointOfField.hxx"
  27 #include "MULTIPR_DecimationFilter.hxx"
  28 #include "MULTIPR_Utils.hxx"
  29 #include "MULTIPR_Exceptions.hxx"
  30 #include "MULTIPR_Globals.hxx"
  31 #include "MULTIPR_API.hxx"
  32 #include <stdio.h>
  33 #include <cmath>
  34
  35 using namespace std;
  36
  37 namespace multipr
  38 {
  39
  40 const med_geometrie_element CELL_TYPES[MED_NBR_GEOMETRIE_MAILLE] =
  41 {
  42     MED_POINT1,
  43     MED_SEG2,
  44     MED_SEG3,
  45     MED_TRIA3,
  46     MED_TRIA6,
  47     MED_QUAD4,
  48     MED_QUAD8,
  49     MED_TETRA4,
  50     MED_TETRA10,
  51     MED_HEXA8,
  52     MED_HEXA20,
  53     MED_PENTA6,
  54     MED_PENTA15,
  55     MED_PYRA5,
  56     MED_PYRA13
  57 };
  58
  59
  60 char CELL_NAMES[MED_NBR_GEOMETRIE_MAILLE][MED_TAILLE_NOM + 1] =
  61 {
  62     "MED_POINT1",
  63     "MED_SEG2",
  64     "MED_SEG3",
  65     "MED_TRIA3",
  66     "MED_TRIA6",
  67     "MED_QUAD4",
  68     "MED_QUAD8",
  69     "MED_TETRA4",
  70     "MED_TETRA10",
  71     "MED_HEXA8",
  72     "MED_HEXA20",
  73     "MED_PENTA6",
  74     "MED_PENTA15",
  75     "MED_PYRA5",
  76     "MED_PYRA13"
  77 };
  78
  79 // Number of points to consider when looking for significant nodes in a mesh.
  80 // ie the n first nodes.
  81 const int CELL_NB_NODE[MED_NBR_GEOMETRIE_MAILLE] =
  82 {
  83     1,  //MED_POINT1
  84     2,  //MED_SEG2
  85     2,  //MED_SEG3
  86     3,  //MED_TRIA3
  87     3,  //MED_TRIA6
  88     4,  //MED_QUAD4
  89     4,  //MED_QUAD8
  90     4,  //MED_TETRA4
  91     4,  //MED_TETRA10
  92     8,  //MED_HEXA8
  93     8,  //MED_HEXA20
  94     6,  //MED_PENTA6
  95     6,  //MED_PENTA15
  96     5,  //MED_PYRA5
  97     5   //MED_PYRA13
  98 };
  99
 100
 101
 102 //*****************************************************************************
 103 // Class Mesh implementation
 104 //*****************************************************************************
 105
 106 Mesh::Mesh()
 107 {
 108     mNodes    = NULL;
 109     for (int i = 0; i < eMaxMedMesh; i++)
 110     {
 111         mElements[i] = NULL;
 112     }
 113
 114     reset();
 115 }
 116
 117
 118 Mesh::~Mesh()
 119 {
 120     reset();
 121 }
 122
 123
 124 void Mesh::reset()
 125 {
 126     mMEDfilename[0] = '\0';
 127     mMEDfile        = 0;
 128
 129     mMeshName[0]    = '\0';
 130     mMeshUName[0]   = '\0';
 131     mMeshDesc[0]    = '\0';
 132     mMeshDim        = -1;
 133     mMeshType       = MED_NON_STRUCTURE;
 134
 135     for (int itDim = 0 ; itDim < 3 ; itDim++)
 136     {
 137         mMeshBBoxMin[itDim] = numeric_limits<med_float>::quiet_NaN();
 138         mMeshBBoxMax[itDim] = numeric_limits<med_float>::quiet_NaN();
 139     }
 140
 141     if (mNodes != NULL)    { delete mNodes;    mNodes = NULL; }
 142
 143     for (int i = 0; i < eMaxMedMesh; i++)
 144     {
 145         if (mElements[i] != NULL)
 146         {
 147             delete mElements[i];
 148             mElements[i] = NULL;
 149         }
 150     }
 151     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
 152     {
 153         delete mFamilies[itFam];
 154     }
 155     mFamilies.clear();
 156     mFamIdToFam.clear();
 157
 158     for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
 159     {
 160         delete mGroups[itGroup];
 161     }
 162     mGroups.clear();
 163     mGroupNameToGroup.clear();
 164
 165     for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 166     {
 167         delete mGaussLoc[itGaussLoc];
 168     }
 169     mGaussLoc.clear();
 170     mGaussLocNameToGaussLoc.clear();
 171
 172     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 173     {
 174         //delete mFields[itField];
 175     }
 176     mFields.clear();
 177
 178     for (unsigned itProfil = 0 ; itProfil < mProfils.size() ; itProfil++)
 179     {
 180         delete mProfils[itProfil];
 181     }
 182     mProfils.clear();
 183
 184     mFlagPrintAll = false;
 185
 186     mGaussIndex.clear();
 187 }
 188
 189
 190 vector<string> Mesh::getNameScalarFields() const
 191 {
 192     vector<string> res;
 193
 194     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 195     {
 196         Field* currentField = mFields[itField];
 197
 198         // only get scalar fields, not vectorial fields
 199         // (because, currently, decimation can only be performed on scalar fields)
 200         if (currentField->getNumberOfComponents() == 1)
 201         {
 202             res.push_back(currentField->getName());
 203         }
 204     }
 205
 206     return res;
 207 }
 208
 209
 210 int Mesh::getTimeStamps(const char* pFieldName) const
 211 {
 212     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 213     {
 214         Field* currentField = mFields[itField];
 215         if (strcmp(currentField->getName(), pFieldName) == 0)
 216         {
 217             return currentField->getNumberOfTimeSteps();
 218         }
 219     }
 220
 221     return 0;
 222 }
 223
 224 Field* Mesh::getFieldByName(const char* pFieldName, eMeshType pGeomType) const
 225 {
 226     if (pFieldName == NULL) throw NullArgumentException("pFieldName should not be NULL", __FILE__, __LINE__);
 227
 228     Field* retField = NULL;
 229
 230     // for each field
 231     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 232     {
 233         Field* currentField = mFields[itField];
 234         // Check if this is the field we want.
 235         if (strncmp(pFieldName, currentField->getName(), MED_TAILLE_NOM) == 0 &&
 236             (currentField->getGeomIdx() == pGeomType || currentField->isFieldOnNodes()))
 237         {
 238             // field found!
 239             retField = currentField;
 240             break;
 241         }
 242     }
 243
 244     return retField;
 245 }
 246
 247 void Mesh::getFieldMinMax(const char* pFieldName, float& pMin, float& pMax) const
 248 {
 249     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 250     {
 251         Field* currentField = mFields[itField];
 252         // Check if this is the field we want.
 253         if (strncmp(pFieldName, currentField->getName(), MED_TAILLE_NOM) == 0)
 254         {
 255             currentField->getMinMax(pMin, pMax);
 256             break;
 257         }
 258     }
 259 }
 260
 261 GaussLoc* Mesh::getGaussLocByName(const char* pGaussLocName) const
 262 {
 263     if (pGaussLocName == NULL) throw NullArgumentException("pGaussLocName should not be NULL", __FILE__, __LINE__);
 264
 265     map<string, GaussLoc*>::const_iterator itGaussLoc = mGaussLocNameToGaussLoc.find(pGaussLocName);
 266     GaussLoc* retGaussLoc = NULL;
 267
 268     if (itGaussLoc != mGaussLocNameToGaussLoc.end())
 269     {
 270         retGaussLoc = (*itGaussLoc).second;
 271     }
 272
 273     return retGaussLoc;
 274 }
 275
 276
 277 int Mesh::getNumberOfElements(eMeshType pGeomType) const
 278 {
 279     if (mElements[pGeomType])
 280     {
 281         return mElements[pGeomType]->getNumberOfElements();
 282     }
 283     else
 284     {
 285         return 0;
 286     }
 287 }
 288
 289 int Mesh::getNumberOfElements() const
 290 {
 291     int accum = 0;
 292
 293     for (int i = 0; i < eMaxMedMesh; ++i)
 294     {
 295         if (mElements[i])
 296         {
 297             accum += mElements[i]->getNumberOfElements();
 298         }
 299     }
 300     return accum;
 301 }
 302
 303 Profil* Mesh::getProfil(const std::string pProfilName)
 304 {
 305     Profil* profile = NULL;
 306     std::vector<Profil*>::iterator it;
 307     if (pProfilName.size() > 0)
 308     {
 309         for (it = mProfils.begin(); it != mProfils.end(); ++it)
 310         {
 311             if ((*it)->getName() == pProfilName)
 312             {
 313                 profile = (*it);
 314                 break;
 315             }
 316         }
 317     }
 318     return profile;
 319 }
 320
 321 Mesh* Mesh::createFromSetOfElements(const std::set<med_int>* pSetOfElements, const char* pNewMeshName)
 322 {
 323     if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName should not be NULL", __FILE__, __LINE__);
 324
 325     //---------------------------------------------------------------------
 326     // Create a new mesh
 327     //---------------------------------------------------------------------
 328     Mesh* mesh = new Mesh();
 329
 330     //---------------------------------------------------------------------
 331     // Build name of the new mesh
 332     //---------------------------------------------------------------------
 333     strcpy(mesh->mMeshName, pNewMeshName);
 334
 335     MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
 336
 337     //---------------------------------------------------------------------
 338     // Fill general infos
 339     //---------------------------------------------------------------------
 340     strcpy(mesh->mMeshUName, mMeshUName);
 341     strcpy(mesh->mMeshDesc, mMeshDesc);
 342
 343     mesh->mMeshDim = mMeshDim;
 344     mesh->mMeshType = mMeshType;
 345
 346     MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
 347     MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
 348     MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
 349     MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
 350
 351     //---------------------------------------------------------------------
 352     // Build nodes and elements
 353     //---------------------------------------------------------------------
 354     // get all elements involved
 355         for (int i = 0; i < eMaxMedMesh; ++i)
 356         {
 357                 if (pSetOfElements[i].size() != 0)
 358                 {
 359                 mesh->mElements[i] = mElements[i]->extractSubSet(pSetOfElements[i]);
 360                 MULTIPR_LOG((*(mesh->mElements[i])) << endl);
 361                 }
 362         }
 363
 364     // get all nodes involved
 365         set<med_int> setOfNodes;
 366         for (int i = 0; i < eMaxMedMesh; ++i)
 367         {
 368                 if (mElements[i] != NULL && mesh->mElements[i] != NULL)
 369                 {
 370                 const set<med_int>& curSetOfNodes = mesh->mElements[i]->getSetOfNodes();
 371                         setOfNodes.insert(curSetOfNodes.begin(), curSetOfNodes.end());
 372                 }
 373         }
 374     mesh->mNodes = mNodes->extractSubSet(setOfNodes);
 375     MULTIPR_LOG((*(mesh->mNodes)) << endl);
 376
 377     //---------------------------------------------------------------------
 378     // Remap nodes
 379     //---------------------------------------------------------------------
 380         for (int i = 0; i < eMaxMedMesh; ++i)
 381         {
 382                 if (mElements[i] != NULL && mesh->mElements[i] != NULL)
 383                 {
 384                         mesh->mElements[i]->remap(setOfNodes);
 385                         MULTIPR_LOG((*(mesh->mElements[i])) << endl);
 386                 }
 387         }
 388
 389
 390     //---------------------------------------------------------------------
 391     // Build families
 392     //---------------------------------------------------------------------
 393     MULTIPR_LOG("Build fam.:" << endl);
 394     // get families of nodes
 395     {
 396         set<med_int> famOfNodes = mesh->mNodes->getSetOfFamilies();
 397         for (set<med_int>::iterator itFam = famOfNodes.begin() ; itFam != famOfNodes.end() ; itFam++)
 398         {
 399             Family* famSrc = mFamIdToFam[*itFam];
 400             if (mesh->mFamIdToFam.find(famSrc->getId()) == mesh->mFamIdToFam.end())
 401             {
 402                 Family* famDest = famSrc->extractGroup(NULL);
 403                 mesh->mFamilies.push_back(famDest);
 404                 mesh->mFamIdToFam[famDest->getId()] = famDest;
 405             }
 406         }
 407     }
 408
 409     // get families of elements
 410     {
 411         set<med_int> famOfElt;
 412                 for (int i = 0; i < eMaxMedMesh; ++i)
 413                 {
 414                         if (mElements[i] != NULL && mesh->mElements[i] != NULL)
 415                         {
 416                                 set<med_int> curSetOfFamilies = mesh->mElements[i]->getSetOfFamilies();
 417                                 famOfElt.insert(curSetOfFamilies.begin(), curSetOfFamilies.end());
 418                         }
 419                 }
 420         for (set<med_int>::iterator itFam = famOfElt.begin() ; itFam != famOfElt.end() ; itFam++)
 421         {
 422             Family* famSrc = mFamIdToFam[*itFam];
 423             if (mesh->mFamIdToFam.find(famSrc->getId()) == mesh->mFamIdToFam.end())
 424             {
 425                 Family* famDest = famSrc->extractGroup(NULL);
 426                 mesh->mFamilies.push_back(famDest);
 427                 mesh->mFamIdToFam[famDest->getId()] = famDest;
 428             }
 429         }
 430     }
 431
 432     MULTIPR_LOG("Finalize:");
 433
 434     // fill families with elements and build groups
 435     //mesh->finalizeFamiliesAndGroups();
 436
 437     MULTIPR_LOG("OK\n");
 438
 439         //---------------------------------------------------------------------
 440     // Build profils.
 441     //---------------------------------------------------------------------
 442     for (std::vector<Profil*>::iterator it = mProfils.begin(); it != mProfils.end(); ++it)
 443     {
 444         Profil* lProfil = new Profil();
 445         lProfil->create((*it)->getName());
 446         std::set<med_int> lSet;
 447         if ((*it)->getBinding() == OnNodes)
 448         {
 449             (*it)->extractSetOfElement(setOfNodes, lSet);
 450         }
 451         else
 452         {
 453             (*it)->extractSetOfElement(pSetOfElements[lProfil->getGeomIdx()], lSet);
 454         }
 455         if (lSet.size() == 0)
 456         {
 457             delete lProfil;
 458             continue;
 459         }
 460         lProfil->set(lSet);
 461         mesh->mProfils.push_back(lProfil);
 462     }
 463
 464     //---------------------------------------------------------------------
 465     // Create new fields and collect Gauss
 466     //---------------------------------------------------------------------
 467     // for each field
 468     set<string> newSetOfGauss;
 469     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 470     {
 471         Field* currentField = mFields[itField];
 472
 473         Field* newField = NULL;
 474         if (currentField->isFieldOnNodes())
 475         {
 476             newField = currentField->extractSubSet(setOfNodes);
 477         }
 478         else
 479         {
 480                         if (pSetOfElements[currentField->getGeomIdx()].size() != 0)
 481                         {
 482                 newField = currentField->extractSubSet(pSetOfElements[currentField->getGeomIdx()]);
 483                         }
 484         }
 485
 486         if (newField != NULL && !newField->isEmpty())
 487         {
 488             mesh->mFields.push_back(newField);
 489             newField->getSetOfGaussLoc(newSetOfGauss);
 490         }
 491     }
 492     MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
 493
 494     //---------------------------------------------------------------------
 495     // Build Gauss infos
 496     //---------------------------------------------------------------------
 497     for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
 498     {
 499         const string& keyName = (*itSet);
 500
 501         GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
 502         if (gaussLoc != NULL)
 503         {
 504             GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
 505             mesh->mGaussLoc.push_back(copyGaussLoc);
 506             mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
 507         }
 508     }
 509
 510     //---------------------------------------------------------------------
 511     // Compute bbox
 512     //---------------------------------------------------------------------
 513     mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
 514
 515     return mesh;
 516 }
 517
 518
 519 Mesh* Mesh::createFromGroup(const Group* pGroup, const char* pNewMeshName)
 520 {
 521     if (pGroup == NULL) throw NullArgumentException("pGroup should not be NULL", __FILE__, __LINE__);
 522     if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName should not be NULL", __FILE__, __LINE__);
 523     if (strlen(pNewMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("pNewMeshName length too long", __FILE__, __LINE__);
 524     //---------------------------------------------------------------------
 525     // Create a new mesh
 526     //---------------------------------------------------------------------
 527     Mesh* mesh = new Mesh();
 528
 529     //---------------------------------------------------------------------
 530     // Build name of the new mesh
 531     //---------------------------------------------------------------------
 532     strcpy(mesh->mMeshName, pNewMeshName);
 533
 534     MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
 535
 536     //---------------------------------------------------------------------
 537     // Fill general infos
 538     //---------------------------------------------------------------------
 539     strcpy(mesh->mMeshUName, mMeshUName);
 540     strcpy(mesh->mMeshDesc, mMeshDesc);
 541
 542     mesh->mMeshDim = mMeshDim;
 543     mesh->mMeshType = mMeshType;
 544
 545     MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
 546     MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
 547     MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
 548     MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
 549
 550     //---------------------------------------------------------------------
 551     // Build nodes and elements
 552     //---------------------------------------------------------------------
 553     // Get all elements involved
 554         std::set< med_int>* setOfEltList = new std::set< med_int>[eMaxMedMesh];
 555         for (int i = 0; i < eMaxMedMesh; ++i)
 556         {
 557                 if (mElements[i] != NULL)
 558                 {
 559                         const set<med_int> setOfElt = pGroup->getSetOfElt((eMeshType)i);
 560                         mesh->mElements[i] = mElements[i]->extractSubSet(setOfElt);
 561                         setOfEltList[i] = setOfElt;
 562                 }
 563         }
 564
 565     // Get all nodes involved
 566         // The nodes a common for all elements so we don't need to store them separately.
 567         set<med_int> setOfNodes;
 568         for (int i = 0; i < eMaxMedMesh; ++i)
 569         {
 570                 if (mesh->mElements[i] != NULL)
 571                 {
 572                         const set<med_int>& curSetOfNodes = mesh->mElements[i]->getSetOfNodes();
 573                         setOfNodes.insert(curSetOfNodes.begin(), curSetOfNodes.end());
 574                 }
 575         }
 576     mesh->mNodes = mNodes->extractSubSet(setOfNodes);
 577
 578         // We need this for the optimized memory management.
 579         this->mGaussIndex.push_back(IndexPair(setOfEltList, setOfNodes));
 580     //---------------------------------------------------------------------
 581     // Remap nodes
 582     //---------------------------------------------------------------------
 583         for (int i = 0; i < eMaxMedMesh; ++i)
 584         {
 585                 if (mesh->mElements[i] != NULL)
 586                 {
 587                         mesh->mElements[i]->remap(setOfNodes);
 588                 }
 589         }
 590
 591     //---------------------------------------------------------------------
 592     // Build families
 593     //---------------------------------------------------------------------
 594     MULTIPR_LOG("Build fam.:" << endl);
 595     // get families of nodes
 596     {
 597                 set<med_int> famOfNodes = mesh->mNodes->getSetOfFamilies();
 598                 for (set<med_int>::iterator itFam = famOfNodes.begin() ; itFam != famOfNodes.end() ; itFam++)
 599         {
 600             Family* famSrc = mFamIdToFam[*itFam];
 601             Family* famDest = famSrc->extractGroup(pGroup->getName().c_str());
 602             mesh->mFamilies.push_back(famDest);
 603         }
 604         }
 605
 606     // get families of elements
 607     {
 608                 set<med_int> famOfElt;
 609                 for (int i = 0; i < eMaxMedMesh; ++i)
 610                 {
 611                         if (mesh->mElements[i] != NULL)
 612                         {
 613                                 set<med_int> curSetOfFamilies = mesh->mElements[i]->getSetOfFamilies();
 614                                 famOfElt.insert(curSetOfFamilies.begin(), curSetOfFamilies.end());
 615                         }
 616                 }
 617         for (set<med_int>::iterator itFam = famOfElt.begin() ; itFam != famOfElt.end() ; itFam++)
 618         {
 619             Family* famSrc = mFamIdToFam[*itFam];
 620             Family* famDest = famSrc->extractGroup(pGroup->getName().c_str());
 621             mesh->mFamilies.push_back(famDest);
 622         }
 623     }
 624
 625     MULTIPR_LOG("Finalize:");
 626
 627     // fill families with elements and build groups
 628     mesh->finalizeFamiliesAndGroups();
 629
 630     MULTIPR_LOG("OK\n");
 631
 632     //---------------------------------------------------------------------
 633     // Build profils.
 634     //---------------------------------------------------------------------
 635     for (std::vector<Profil*>::iterator it = mProfils.begin(); it != mProfils.end(); ++it)
 636     {
 637         Profil* lProfil = new Profil();
 638         lProfil->create((*it)->getName());
 639         std::set<med_int> lSet;
 640         if ((*it)->getBinding() == OnNodes)
 641         {
 642             (*it)->extractSetOfElement(setOfNodes, lSet);
 643         }
 644         else
 645         {
 646             (*it)->extractSetOfElement(setOfEltList[lProfil->getGeomIdx()], lSet);
 647         }
 648         if (lSet.size() == 0)
 649         {
 650             delete lProfil;
 651             continue;
 652         }
 653         lProfil->set(lSet);
 654         mesh->mProfils.push_back(lProfil);
 655     }
 656
 657     //---------------------------------------------------------------------
 658     // Create new fields and collect Gauss
 659     //---------------------------------------------------------------------
 660     // for each field
 661     set<string> newSetOfGauss;
 662     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 663     {
 664         Field* currentField = mFields[itField];
 665
 666         Field* newField;
 667         if (currentField->isFieldOnNodes())
 668         {
 669             newField = currentField->extractSubSet(setOfNodes);
 670         }
 671         else
 672         {
 673                         if (setOfEltList[currentField->getGeomIdx()].size() != 0)
 674                         {
 675                 newField = currentField->extractSubSet(setOfEltList[currentField->getGeomIdx()]);
 676                         }
 677         }
 678
 679         if (!newField->isEmpty())
 680         {
 681             mesh->mFields.push_back(newField);
 682             newField->getSetOfGaussLoc(newSetOfGauss);
 683         }
 684     }
 685         // Get gauss locs for optimized fields reading.
 686         if (mFields.size() == 0)
 687         {
 688                 for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 689                 {
 690                         const string& gaussLocName = mGaussLoc[itGaussLoc]->getName();
 691
 692                         if (gaussLocName.length() != 0)
 693                         {
 694                                 newSetOfGauss.insert(gaussLocName);
 695                         }
 696                 }
 697         }
 698     MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
 699
 700     //---------------------------------------------------------------------
 701     // Build Gauss infos
 702     //---------------------------------------------------------------------
 703         for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
 704     {
 705         const string& keyName = (*itSet);
 706
 707         GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
 708         if (gaussLoc != NULL)
 709         {
 710             GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
 711             mesh->mGaussLoc.push_back(copyGaussLoc);
 712             mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
 713         }
 714     }
 715
 716     //---------------------------------------------------------------------
 717     // Compute bbox
 718     //---------------------------------------------------------------------
 719     mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
 720
 721     return mesh;
 722 }
 723
 724
 725 Mesh* Mesh::createFromFamily(const Family* pFamily, const char* pNewMeshName)
 726 {
 727     if (pFamily == NULL) throw NullArgumentException("pFamily should not be NULL", __FILE__, __LINE__);
 728     if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName should not be NULL", __FILE__, __LINE__);
 729     if (strlen(pNewMeshName) > MED_TAILLE_NOM)
 730     {
 731         char msg[256];
 732         sprintf(msg, "pNewMeshName length (=%d) too long: max size allowed is %d", strlen(pNewMeshName), MED_TAILLE_NOM);
 733         throw IllegalArgumentException(msg, __FILE__, __LINE__);
 734     }
 735     //---------------------------------------------------------------------
 736     // Create a new mesh
 737     //---------------------------------------------------------------------
 738     Mesh* mesh = new Mesh();
 739
 740     //---------------------------------------------------------------------
 741     // Build name of the new mesh
 742     //---------------------------------------------------------------------
 743     strcpy(mesh->mMeshName, pNewMeshName);
 744
 745     MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
 746
 747     //---------------------------------------------------------------------
 748     // Fill general infos
 749     //---------------------------------------------------------------------
 750     strcpy(mesh->mMeshUName, mMeshUName);
 751     strcpy(mesh->mMeshDesc, mMeshDesc);
 752
 753     mesh->mMeshDim = mMeshDim;
 754     mesh->mMeshType = mMeshType;
 755
 756     MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
 757     MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
 758     MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
 759     MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
 760
 761     //---------------------------------------------------------------------
 762     // Build nodes and elements
 763     //---------------------------------------------------------------------
 764     // Get all elements involved
 765         std::set< med_int>* setOfEltList = new std::set< med_int>[eMaxMedMesh];
 766         for (int i = 0; i < eMaxMedMesh; ++i)
 767         {
 768                 if (mElements[i] != NULL)
 769                 {
 770                         const set<med_int> setOfElt = pFamily->getSetOfElt((eMeshType)i);
 771                         mesh->mElements[i] = mElements[i]->extractSubSet(setOfElt);
 772                         setOfEltList[i] = setOfElt;
 773                 }
 774         }
 775
 776     // Get all nodes involved
 777         // The nodes a common for all elements so we don't need to store them separately.
 778         set<med_int> setOfNodes;
 779         for (int i = 0; i < eMaxMedMesh; ++i)
 780         {
 781                 if (mesh->mElements[i] != NULL)
 782                 {
 783                         const set<med_int>& curSetOfNodes = mesh->mElements[i]->getSetOfNodes();
 784                         setOfNodes.insert(curSetOfNodes.begin(), curSetOfNodes.end());
 785                 }
 786         }
 787     mesh->mNodes = mNodes->extractSubSet(setOfNodes);
 788
 789         // We need this for the optimized memory management.
 790         this->mGaussIndex.push_back(IndexPair(setOfEltList, setOfNodes));
 791     //---------------------------------------------------------------------
 792     // Remap nodes
 793     //---------------------------------------------------------------------
 794         for (int i = 0; i < eMaxMedMesh; ++i)
 795         {
 796                 if (mesh->mElements[i] != NULL)
 797                 {
 798                         mesh->mElements[i]->remap(setOfNodes);
 799                 }
 800         }
 801
 802     //---------------------------------------------------------------------
 803     // Build families
 804     //---------------------------------------------------------------------
 805     MULTIPR_LOG("Build fam.:" << endl);
 806     // get families of nodes
 807         Family*    lFam = new Family(*pFamily);
 808     mesh->mFamilies.push_back(lFam);
 809
 810     MULTIPR_LOG("Finalize:");
 811
 812     // fill families with elements and build groups
 813     mesh->finalizeFamiliesAndGroups();
 814
 815     MULTIPR_LOG("OK\n");
 816
 817     //---------------------------------------------------------------------
 818     // Build profils.
 819     //---------------------------------------------------------------------
 820     for (std::vector<Profil*>::iterator it = mProfils.begin(); it != mProfils.end(); ++it)
 821     {
 822         Profil* lProfil = new Profil();
 823         lProfil->create((*it)->getName());
 824         std::set<med_int> lSet;
 825         if ((*it)->getBinding() == OnNodes)
 826         {
 827             (*it)->extractSetOfElement(setOfNodes, lSet);
 828         }
 829         else
 830         {
 831             (*it)->extractSetOfElement(setOfEltList[lProfil->getGeomIdx()], lSet);
 832         }
 833         if (lSet.size() == 0)
 834         {
 835             delete lProfil;
 836             continue;
 837         }
 838         lProfil->set(lSet);
 839         mesh->mProfils.push_back(lProfil);
 840     }
 841
 842     //---------------------------------------------------------------------
 843     // Create new fields and collect Gauss
 844     //---------------------------------------------------------------------
 845     // for each field
 846     set<string> newSetOfGauss;
 847     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 848     {
 849         Field* currentField = mFields[itField];
 850
 851         Field* newField;
 852         if (currentField->isFieldOnNodes())
 853         {
 854             newField = currentField->extractSubSet(setOfNodes);
 855         }
 856         else
 857         {
 858                         if (setOfEltList[currentField->getGeomIdx()].size() != 0)
 859                         {
 860                 newField = currentField->extractSubSet(setOfEltList[currentField->getGeomIdx()]);
 861                         }
 862         }
 863
 864         if (!newField->isEmpty())
 865         {
 866             mesh->mFields.push_back(newField);
 867             newField->getSetOfGaussLoc(newSetOfGauss);
 868         }
 869     }
 870         // Get gauss locs for optimized fields reading.
 871         if (mFields.size() == 0)
 872         {
 873                 for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 874                 {
 875                         const string& gaussLocName = mGaussLoc[itGaussLoc]->getName();
 876
 877                         if (gaussLocName.length() != 0)
 878                         {
 879                                 newSetOfGauss.insert(gaussLocName);
 880                         }
 881                 }
 882         }
 883     MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
 884
 885     //---------------------------------------------------------------------
 886     // Build Gauss infos
 887     //---------------------------------------------------------------------
 888         for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
 889     {
 890         const string& keyName = (*itSet);
 891
 892         GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
 893         if (gaussLoc != NULL)
 894         {
 895             GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
 896             mesh->mGaussLoc.push_back(copyGaussLoc);
 897             mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
 898         }
 899     }
 900
 901     //---------------------------------------------------------------------
 902     // Compute bbox
 903     //---------------------------------------------------------------------
 904     mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
 905
 906     return mesh;
 907 }
 908
 909 Mesh* Mesh::mergePartial(vector<Mesh*> pMeshes, const char* pFieldName, int pFieldIt)
 910 {
 911     if (pMeshes.size() == 0) throw IllegalArgumentException("list must contain one mesh at least", __FILE__, __LINE__);
 912
 913     //---------------------------------------------------------------------
 914     // Create a new mesh
 915     //---------------------------------------------------------------------
 916     Mesh* mesh = new Mesh();
 917
 918     //---------------------------------------------------------------------
 919     // Build name of the new mesh
 920     //---------------------------------------------------------------------
 921     strcpy(mesh->mMeshName, mMeshName);
 922
 923     //---------------------------------------------------------------------
 924     // Merge general infos
 925     //---------------------------------------------------------------------
 926     strcpy(mesh->mMeshUName, mMeshUName);
 927     strcpy(mesh->mMeshDesc, mMeshDesc);
 928
 929     mesh->mMeshDim  = mMeshDim;
 930     mesh->mMeshType = mMeshType;
 931
 932     //---------------------------------------------------------------------
 933     // Merge nodes and elements
 934     //---------------------------------------------------------------------
 935     vector<Nodes*>    nodes;
 936     vector<Elements*> elements[eMaxMedMesh];
 937     vector<int>       offsets[eMaxMedMesh];
 938
 939     int offset[eMaxMedMesh];
 940     for (unsigned j = 0; j < eMaxMedMesh; ++j)
 941     {
 942         offset[j] = mNodes->getNumberOfNodes();
 943     }
 944
 945     for (unsigned i = 0 ; i < pMeshes.size() ; i++)
 946     {
 947         if (mMeshDim != pMeshes[i]->mMeshDim) throw IllegalStateException("all meshes should have same dimension", __FILE__, __LINE__);
 948         if (mMeshType != pMeshes[i]->mMeshType) throw IllegalStateException("all meshes should have same type", __FILE__, __LINE__);
 949
 950
 951         nodes.push_back(pMeshes[i]->mNodes);
 952         for (unsigned j = 0; j < eMaxMedMesh; ++j)
 953         {
 954
 955             if (pMeshes[i]->mElements[j] != NULL)
 956             {
 957                 elements[j].push_back(pMeshes[i]->mElements[j]);
 958                 offsets[j].push_back(offset[j]);
 959             }
 960             offset[j] += pMeshes[i]->mNodes->getNumberOfNodes();
 961         }
 962     }
 963
 964     mesh->mNodes = mNodes->mergePartial(nodes);
 965     for (unsigned i = 0; i < eMaxMedMesh; ++i)
 966     {
 967         if (elements[i].size() != 0)
 968         {
 969             if (mElements[i] == NULL)
 970             {
 971                 mElements[i] = new Elements(CELL_TYPES[i]);
 972                 mElements[i]->mergePartial(elements[i].front(), offsets[i].front());
 973                 elements[i].erase(elements[i].begin());
 974             }
 975             mesh->mElements[i] = mElements[i]->mergePartial(elements[i], offsets[i]);
 976         }
 977         else if (mElements[i] != NULL)
 978         {
 979             mesh->mElements[i] = mElements[i];
 980         }
 981     }
 982
 983
 984     //---------------------------------------------------------------------
 985     // Merge families
 986     //---------------------------------------------------------------------
 987     for (unsigned i = 0 ; i < mFamilies.size() ; i++)
 988     {
 989         Family* family = new Family(*(mFamilies[i]));
 990         mesh->mFamilies.push_back(family);
 991         mesh->mFamIdToFam.insert(make_pair(family->getId(), family));
 992     }
 993
 994     for (unsigned j = 0 ; j < pMeshes.size() ; j++)
 995     {
 996         for (unsigned i = 0 ; i < pMeshes[j]->mFamilies.size() ; i++)
 997         {
 998             // test if there is a fimaly with the same id
 999             map<med_int, Family*>::iterator itFam = mesh->mFamIdToFam.find(pMeshes[j]->mFamilies[i]->getId());
1000
1001             if (itFam == mesh->mFamIdToFam.end())
1002             {
1003                 // id not found: create a new family
1004                 Family* family = new Family(*(pMeshes[j]->mFamilies[i]));
1005                 mesh->mFamilies.push_back(family);
1006                 mesh->mFamIdToFam.insert(make_pair(family->getId(), family));
1007             }
1008         }
1009     }
1010
1011     //---------------------------------------------------------------------
1012     // Merge fields
1013     //---------------------------------------------------------------------
1014     // for each mesh
1015
1016     vector<Field*> fields;
1017     for (unsigned i = 0 ; i < pMeshes.size() ; i++)
1018     {
1019         for (std::vector<Field*>::iterator it = pMeshes[i]->mFields.begin();
1020             it != pMeshes[i]->mFields.end(); ++it)
1021         {
1022             if (strcmp((*it)->getName(), pFieldName) == 0)
1023             {
1024                 fields.push_back(*it);
1025                 break;
1026             }
1027         }
1028     }
1029     bool    hasField = false;
1030     for (std::vector<Field*>::iterator it = mFields.begin();
1031         it != mFields.end(); ++it)
1032     {
1033         if (strcmp((*it)->getName(), pFieldName) == 0)
1034         {
1035             Field* field = (*it)->merge(fields, pFieldIt);
1036             mesh->mFields.push_back(field);
1037             hasField = true;
1038             break;
1039         }
1040     }
1041     if (hasField == false)
1042     {
1043         Field*  lField = fields.back();
1044         fields.pop_back();
1045         Field* field = lField->merge(fields, pFieldIt);
1046         mesh->mFields.push_back(field);
1047     }
1048
1049     //---------------------------------------------------------------------
1050     // Merge Gauss infos
1051     //---------------------------------------------------------------------
1052     // WARNING: assume Gauss infos are the same for the two meshes.
1053     for (unsigned i = 0 ; i < mGaussLoc.size() ; i++)
1054     {
1055         GaussLoc* copyGaussLoc = new GaussLoc(*(mGaussLoc[i]));
1056         mesh->mGaussLoc.push_back(copyGaussLoc);
1057         mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
1058     }
1059
1060     return mesh;
1061 }
1062
1063
1064 MeshDis* Mesh::splitGroupsOfElements()
1065 {
1066     MeshDis* meshDis = new MeshDis();
1067     meshDis->setSequentialMEDFilename(mMEDfilename);
1068
1069     // get prefix from the original MED filename
1070     string strPrefix = removeExtension(mMEDfilename, ".med");
1071     int numGroup = 1;
1072         int numFam = 1;
1073
1074     // for each group
1075     for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
1076     {
1077         Group* currentGroup = mGroups[itGroup];
1078
1079         // skip this group if it is a group of nodes
1080         if (currentGroup->isGroupOfNodes())
1081         {
1082                         this->mGaussIndex.push_back(IndexPair());
1083                         continue;
1084         }
1085
1086         char strPartName[256];
1087         sprintf(strPartName, "%s_%d", mMeshName, numGroup);
1088
1089         char strMEDfilename[256];
1090                 char strMedGroup[256];
1091                 int i;
1092                 for (i = 0; currentGroup->getName()[i] && currentGroup->getName()[i] != ' '; ++i)
1093                 {
1094                         strMedGroup[i] = currentGroup->getName()[i];
1095                 }
1096                 strMedGroup[i] = '\0';
1097         sprintf(strMEDfilename, "%s_%s.med", strPrefix.c_str(), strMedGroup);
1098
1099         Mesh* mesh = createFromGroup(currentGroup, mMeshName);
1100
1101         // skip the group which contain all the others groups, even if it contains only 1 group
1102                 if (mesh->getNumberOfElements() == this->getNumberOfElements())
1103         {
1104                         for (int i = 0; i < eMaxMedMesh; ++i)
1105                         {
1106                                 this->mGaussIndex.back().first[i].clear();
1107                         }
1108                         this->mGaussIndex.back().second.clear();
1109             delete mesh;
1110                         mesh = NULL;
1111                         continue;
1112                 }
1113         meshDis->addMesh(
1114             MeshDisPart::MULTIPR_WRITE_MESH,
1115             mMeshName,
1116             numGroup,
1117             strPartName,
1118             "localhost",
1119             strMEDfilename,
1120             mesh);
1121
1122         numGroup++;
1123     }
1124         if (mGroups.size() == 0)
1125         {
1126                 for (unsigned itFam = 0; itFam < mFamilies.size(); ++itFam)
1127                 {
1128
1129                         Family* currentFam = mFamilies[itFam];
1130
1131                         // skip this family if it is a family of nodes
1132                         if (currentFam->isFamilyOfNodes())
1133                         {
1134                                 this->mGaussIndex.push_back(IndexPair());
1135                                 continue;
1136                         }
1137
1138                         char strPartName[256];
1139                         sprintf(strPartName, "%s_%d", mMeshName, numGroup);
1140
1141                         char strMEDfilename[256];
1142                         char strMedFam[256];
1143                         int i;
1144                         for (i = 0; currentFam->getName()[i] && currentFam->getName()[i] != ' '; ++i)
1145                         {
1146                                 strMedFam[i] = currentFam->getName()[i];
1147                         }
1148                         strMedFam[i] = '\0';
1149                         sprintf(strMEDfilename, "%s_%s.med", strPrefix.c_str(), strMedFam);
1150
1151                         Mesh* mesh = createFromFamily(currentFam, mMeshName);
1152
1153                         meshDis->addMesh(
1154                                 MeshDisPart::MULTIPR_WRITE_MESH,
1155                                 mMeshName,
1156                                 numFam,
1157                                 strPartName,
1158                                 "localhost",
1159                                 strMEDfilename,
1160                                 mesh);
1161
1162                         numFam++;
1163                 }
1164         }
1165
1166     return meshDis;
1167 }
1168
1169
1170 Mesh* Mesh::decimate(
1171     const char* pFilterName,
1172     const char* pArgv,
1173     const char* pNameNewMesh)
1174 {
1175     //---------------------------------------------------------------------
1176     // Check parameters
1177     //---------------------------------------------------------------------
1178     if (pFilterName == NULL) throw NullArgumentException("pFilterName should not be NULL", __FILE__, __LINE__);
1179     if (pArgv == NULL) throw NullArgumentException("pArgv should not be NULL", __FILE__, __LINE__);
1180     if (pNameNewMesh == NULL) throw NullArgumentException("pNameNewMesh should not be NULL", __FILE__, __LINE__);
1181
1182     //---------------------------------------------------------------------
1183     // Instanciate filter used for decimation
1184     //---------------------------------------------------------------------
1185     DecimationFilter* filter = DecimationFilter::create(pFilterName);
1186
1187     //---------------------------------------------------------------------
1188     // Create new mesh by decimating current one
1189     //---------------------------------------------------------------------
1190     Mesh* decimatedMesh = filter->apply(this, pArgv, pNameNewMesh);
1191
1192     //---------------------------------------------------------------------
1193     // Cleans
1194     //---------------------------------------------------------------------
1195     delete filter;
1196
1197     return decimatedMesh;
1198 }
1199
1200
1201
1202 void Mesh::getAllPointsOfField(Field* pField, int pTimeStepIt, std::vector<PointOfField>& pPoints, eMeshType pGeomType)
1203 {
1204     //---------------------------------------------------------------------
1205     // Check arguments
1206     //---------------------------------------------------------------------
1207     if (pField == NULL) throw NullArgumentException("field should not be NULL", __FILE__, __LINE__);
1208     if (pTimeStepIt < 1) throw IllegalArgumentException("invalid field iteration; should be >= 1", __FILE__, __LINE__);
1209
1210     if (mMeshDim != 3) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
1211     if (pField->getType() != MED_FLOAT64) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
1212     if (pField->getNumberOfComponents() != 1) throw UnsupportedOperationException("field have more than 1 component (vectorial field, expected scalar field)", __FILE__, __LINE__);
1213
1214     //---------------------------------------------------------------------
1215     // Get profile
1216     //---------------------------------------------------------------------
1217
1218     const std::string&  profilName = pField->getProfil(pTimeStepIt);
1219     std::vector<Profil*>::iterator it;
1220     Profil* profile = NULL;
1221     int count = 0;
1222     if (profilName.size() > 0)
1223     {
1224         for (it = mProfils.begin(); it != mProfils.end(); ++it)
1225         {
1226             if ((*it)->getName() == profilName)
1227             {
1228                 profile = (*it);
1229                 break;
1230             }
1231         }
1232         if (it == mProfils.end()) throw IllegalStateException("Can't find the profile in the mesh.", __FILE__, __LINE__);
1233
1234     }
1235
1236     //---------------------------------------------------------------------
1237     // Collect points
1238     //---------------------------------------------------------------------
1239
1240     if (pField->isFieldOnNodes())
1241     {
1242         //-------------------------------------------------------------
1243         // Case 1: field of nodes
1244         //-------------------------------------------------------------
1245         if (mNodes == NULL) throw IllegalStateException("no nodes in the current mesh", __FILE__, __LINE__);
1246
1247         // for each node
1248         for (int itNode = 0, size = mNodes->getNumberOfNodes() ; itNode < size ; itNode++)
1249         {
1250             if (profile && profile->find(itNode) == false)
1251             {
1252                 continue;
1253             }
1254             // collect coordinates and value of the point
1255             const med_float* coo = mNodes->getCoordinates(itNode);
1256
1257             const med_float* val =
1258                 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, count + 1));
1259             // add new point
1260             pPoints.push_back(PointOfField(coo[0], coo[1], coo[2], val[0]));
1261             ++count;
1262         }
1263     }
1264     else
1265     {
1266         //-------------------------------------------------------------
1267         // Case 2: field of elements
1268         //-------------------------------------------------------------
1269
1270         if (mElements[pGeomType] == NULL) throw IllegalStateException("no elements in the current mesh", __FILE__, __LINE__);
1271
1272         const string& nameGaussLoc = pField->getNameGaussLoc(pTimeStepIt);
1273         GaussLoc* gaussLoc = getGaussLocByName(nameGaussLoc.c_str());
1274         if (gaussLoc == NULL) throw IllegalStateException("no Gauss localization for these elements", __FILE__, __LINE__);
1275
1276         int numGauss = pField->getNumberOfGaussPointsByElement(pTimeStepIt);
1277         int size = gaussLoc->getDim() * gaussLoc->getNumGaussPoints();
1278         med_float* cooGaussPts = new med_float[size];
1279
1280         int dim = mElements[pGeomType]->getTypeOfPrimitives() / 100;
1281         //int numNodes = mElements[pGeomType]->getTypeOfPrimitives() % 100;
1282         size = dim * numGauss;
1283         med_float* cooNodes = new med_float[size];
1284
1285         // for each elements
1286         for (int itElt = 0, size = mElements[pGeomType]->getNumberOfElements() ; itElt < size ; itElt++)
1287         {
1288             if (profile && profile->find(itElt) == false)
1289             {
1290                 continue;
1291             }
1292             // get coordinates of nodes of the current elements
1293             mElements[pGeomType]->getCoordinates(itElt, mNodes, cooNodes, numGauss);
1294
1295             // compute coordinates of gauss points
1296             gaussLoc->getCoordGaussPoints(cooNodes, cooGaussPts, numGauss);
1297
1298             const med_float* val =
1299                 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, count + 1));
1300
1301             // for each point of Gauss of the element
1302             med_float* srcCoo = cooGaussPts;
1303             for (int itPtGauss = 0 ; itPtGauss < numGauss ; itPtGauss++)
1304             {
1305                 pPoints.push_back(PointOfField(srcCoo[0], srcCoo[1], srcCoo[2], val[itPtGauss]));
1306                 srcCoo += 3;
1307             }
1308             ++count;
1309         }
1310
1311         delete[] cooNodes;
1312         delete[] cooGaussPts;
1313     }
1314 }
1315
1316
1317 float Mesh::evalDefaultRadius(int pN) const
1318 {
1319     if (mFields.size() == 0) return 1.0f;
1320
1321     //---------------------------------------------------------------------
1322     // Compute default radius
1323     //---------------------------------------------------------------------
1324
1325     med_float volumeBBox =
1326         (mMeshBBoxMax[0] - mMeshBBoxMin[0]) *
1327         (mMeshBBoxMax[1] - mMeshBBoxMin[1]) *
1328         (mMeshBBoxMax[2] - mMeshBBoxMin[2]);
1329
1330     if (isnan(volumeBBox))
1331     {
1332         return 1.0f;
1333     }
1334
1335     const med_float k = 0.8; // considered 80% of the volume
1336
1337     // get number of gauss points in the field
1338     try
1339     {
1340         Field* anyField = mFields[mFields.size()-1];
1341         int numTimeSteps = anyField->getNumberOfTimeSteps();
1342
1343         int numGaussPoints = getNumberOfElements() * anyField->getNumberOfGaussPointsByElement(numTimeSteps-1);
1344
1345         med_float radius = med_float(pow( (3.0/4.0) * pN * k * volumeBBox / (3.1415 * numGaussPoints), 1.0/3.0));
1346
1347         return float(radius);
1348     }
1349     catch (...)
1350     {
1351         return 1.0f;
1352     }
1353 }
1354
1355 void Mesh::_openMEDFile(const char* pMEDfilename, med_mode_acces pMEDModeAccess)
1356 {
1357     reset();
1358
1359     //---------------------------------------------------------------------
1360     // Check arguments
1361     //---------------------------------------------------------------------
1362     MULTIPR_LOG("Check arguments: ");
1363     if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
1364     MULTIPR_LOG("OK\n");
1365
1366     strncpy(mMEDfilename, pMEDfilename, 256);
1367
1368     //---------------------------------------------------------------------
1369     // Open MED file (READ_ONLY)
1370     //---------------------------------------------------------------------
1371     MULTIPR_LOG("Open MED file: ");
1372     mMEDfile = MEDouvrir(mMEDfilename, pMEDModeAccess);
1373     if (mMEDfile <= 0) throw FileNotFoundException("MED file not found", __FILE__, __LINE__);
1374     MULTIPR_LOG("OK\n");
1375
1376     //---------------------------------------------------------------------
1377     // Check valid HDF format
1378     //---------------------------------------------------------------------
1379     MULTIPR_LOG("Format HDF: ");
1380     if (MEDformatConforme(mMEDfilename) != 0) throw IOException("invalid file", __FILE__, __LINE__);
1381     MULTIPR_LOG("OK\n");
1382
1383     //---------------------------------------------------------------------
1384     // Get MED version
1385     //---------------------------------------------------------------------
1386     MULTIPR_LOG("MED version: ");
1387     med_int verMajor, verMinor, verRelease;
1388     med_err ret = MEDversionLire(mMEDfile, &verMajor, &verMinor, &verRelease);
1389     if (ret != 0) throw IOException("error while reading MED version", __FILE__, __LINE__);
1390     MULTIPR_LOG(verMajor << "." << verMinor << "." << verRelease << ": OK\n");
1391   }
1392
1393
1394 void Mesh::_readSequentialMED(const char* pMeshName, bool pReadFields)
1395 {
1396     med_err ret = 1;
1397
1398     //---------------------------------------------------------------------
1399     // Check arguments
1400     //---------------------------------------------------------------------
1401     MULTIPR_LOG("Check arguments: ");
1402     if (pMeshName == NULL) throw NullArgumentException("pMeshName should not be NULL", __FILE__, __LINE__);
1403     MULTIPR_LOG("OK\n");
1404
1405     strncpy(mMeshName, pMeshName, MED_TAILLE_NOM);
1406
1407     //---------------------------------------------------------------------
1408     // Read profils
1409     //---------------------------------------------------------------------
1410     MULTIPR_LOG("Profils: ");
1411     med_int nbProfils = MEDnProfil(mMEDfile);
1412     for (med_int i = 1; i <= nbProfils; ++i)
1413     {
1414         Profil* profil = new Profil();
1415         profil->readMED(mMEDfile, i);
1416         profil->readProfilBinding(mMEDfile, const_cast<char*>(pMeshName));
1417         this->mProfils.push_back(profil);
1418     }
1419     MULTIPR_LOG("OK\n");
1420
1421     //---------------------------------------------------------------------
1422     // Read all Gauss localizations
1423     //---------------------------------------------------------------------
1424     MULTIPR_LOG("Gauss localizations: ");
1425     readGaussLoc();
1426
1427     //---------------------------------------------------------------------
1428     // Read scalars (should be 0)
1429     //---------------------------------------------------------------------
1430     MULTIPR_LOG("Scalars: ");
1431     med_int nbScalars = MEDnScalaire(mMEDfile);
1432     if (nbScalars != 0) throw UnsupportedOperationException("scalars information not yet implemented", __FILE__, __LINE__);
1433     MULTIPR_LOG(nbScalars << ": OK\n");
1434
1435     //---------------------------------------------------------------------
1436     // Find the mesh
1437     //---------------------------------------------------------------------
1438     // read number of meshes
1439     MULTIPR_LOG("Num meshes: ");
1440     med_int nbMeshes = MEDnMaa(mMEDfile);
1441     if (nbMeshes <= 0) throw IOException("i/o error while reading number of meshes in MED file", __FILE__, __LINE__);
1442     MULTIPR_LOG(nbMeshes << ": OK\n");
1443
1444     med_int meshIndex = -1;
1445     // iteration over mesh to find the mesh we want
1446     // for each mesh in the file (warning: first mesh is number 1)
1447     for (int itMesh = 1 ; itMesh <= nbMeshes ; itMesh++)
1448     {
1449         char meshName[MED_TAILLE_NOM + 1];
1450
1451         ret = MEDmaaInfo(
1452             mMEDfile,
1453             itMesh,
1454             meshName,
1455             &mMeshDim,
1456             &mMeshType,
1457             mMeshDesc);
1458
1459         if (ret != 0) throw IOException("i/o error while reading mesh information in MED file", __FILE__, __LINE__);
1460         MULTIPR_LOG("Mesh: |" << meshName << "|");
1461
1462         // test if the current mesh is the mesh we want
1463         if (strcmp(pMeshName, meshName) == 0)
1464         {
1465             // *** mesh found ***
1466             MULTIPR_LOG(" OK (found)\n");
1467             meshIndex = itMesh;
1468             break;
1469         }
1470         else
1471         {
1472             // not the mesh we want: skip this mesh
1473             MULTIPR_LOG(" skipped\n");
1474         }
1475     }
1476
1477     if (meshIndex == -1)
1478     {
1479         throw IllegalStateException("mesh not found in the given MED file", __FILE__, __LINE__);
1480     }
1481
1482     //---------------------------------------------------------------------
1483     // Check mesh validity
1484     //---------------------------------------------------------------------
1485     // dimension of the mesh must be 3 (= 3D mesh)
1486     MULTIPR_LOG("Mesh is 3D: ");
1487     if (mMeshDim != 3) throw UnsupportedOperationException("dimension of the mesh should be 3; other dimension not yet implemented", __FILE__, __LINE__);
1488     MULTIPR_LOG("OK\n");
1489
1490     // mesh must not be a grid
1491     MULTIPR_LOG("Mesh is not a grid: ");
1492     if (mMeshType != MED_NON_STRUCTURE)
1493         throw UnsupportedOperationException("grid not supported", __FILE__, __LINE__);
1494     MULTIPR_LOG("OK\n");
1495
1496     med_connectivite connectivite = MED_NOD; // NODAL CONNECTIVITY ONLY
1497         // Read all the supported geometry type.
1498         for (int itCell = 0; itCell < eMaxMedMesh ; ++itCell)
1499     {
1500         med_int meshNumCells = MEDnEntMaa(
1501             mMEDfile,
1502             mMeshName,
1503             MED_CONN,
1504             MED_MAILLE,
1505             CELL_TYPES[itCell],
1506             connectivite);
1507
1508
1509                 //---------------------------------------------------------------------
1510                 // Read elements
1511                 //---------------------------------------------------------------------
1512         if (meshNumCells > 0)
1513         {
1514                         mElements[itCell] = new Elements();
1515                         mElements[itCell]->readMED(mMEDfile, mMeshName, mMeshDim, MED_MAILLE, CELL_TYPES[itCell]);
1516         }
1517         }
1518     // everything is OK...
1519
1520     //---------------------------------------------------------------------
1521     // Check num joint = 0
1522     //---------------------------------------------------------------------
1523     MULTIPR_LOG("Num joints: ");
1524     med_int numJoints = MEDnJoint(mMEDfile, mMeshName);
1525     MULTIPR_LOG(numJoints << ": OK\n");
1526
1527     //---------------------------------------------------------------------
1528     // Check num equivalence = 0
1529     //---------------------------------------------------------------------
1530     MULTIPR_LOG("Num equivalences: ");
1531     med_int numEquiv = MEDnEquiv(mMEDfile, mMeshName);
1532     MULTIPR_LOG(numEquiv << ": OK\n");
1533
1534     //---------------------------------------------------------------------
1535     // Read nodes
1536     //---------------------------------------------------------------------
1537     mNodes = new Nodes();
1538     mNodes->readMED(mMEDfile, mMeshName, mMeshDim);
1539     mNodes->getBBox(mMeshBBoxMin, mMeshBBoxMax);
1540
1541     if (mNodes->getNumberOfNodes() != 0)
1542     {
1543
1544         //---------------------------------------------------------------------
1545         // Read families
1546         //---------------------------------------------------------------------
1547         readFamilies();
1548         finalizeFamiliesAndGroups();
1549
1550         //---------------------------------------------------------------------
1551         // Read fields
1552         //---------------------------------------------------------------------
1553         if (pReadFields)
1554                 {
1555                         readFields();
1556                 }
1557     }
1558
1559     //---------------------------------------------------------------------
1560     // Close the MED file
1561     //---------------------------------------------------------------------
1562     MULTIPR_LOG("Close MED file: ");
1563     ret = MEDfermer(mMEDfile);
1564     if (ret != 0) throw IOException("i/o error while closing MED file", __FILE__, __LINE__);
1565     MULTIPR_LOG("OK\n");
1566 }
1567
1568
1569 void Mesh::readSequentialMED(const char* pMEDfilename, const char* pMeshName, bool pReadFields)
1570 {
1571     _openMEDFile(pMEDfilename);
1572     _readSequentialMED(pMeshName, pReadFields);
1573 }
1574
1575
1576 void Mesh::readSequentialMED(const char* pMEDfilename, med_int pMeshNumber, bool pReadFields)
1577 {
1578     _openMEDFile(pMEDfilename);
1579
1580     //---------------------------------------------------------------------
1581     // Find the mesh
1582     //---------------------------------------------------------------------
1583     // read number of meshes
1584     MULTIPR_LOG("Num meshes: ");
1585     med_int nbMeshes = MEDnMaa(mMEDfile);
1586     if (nbMeshes <= 0) throw IOException("i/o error while reading number of meshes in MED file", __FILE__, __LINE__);
1587     MULTIPR_LOG(nbMeshes << ": OK\n");
1588
1589     med_int meshDim;
1590     med_maillage meshType;
1591     char meshDesc[MED_TAILLE_DESC + 1];
1592     char meshName[MED_TAILLE_NOM + 1];
1593
1594     med_err ret = MEDmaaInfo(mMEDfile, pMeshNumber, meshName, &meshDim, &meshType, meshDesc);
1595     if (ret != 0) throw IOException("i/o error while reading mesh information in MED file", __FILE__, __LINE__);
1596
1597     _readSequentialMED(meshName, pReadFields);
1598 }
1599
1600
1601 void Mesh::writeMED(const char* pMEDfilename)
1602 {
1603     writeMED(pMEDfilename, getName());
1604 }
1605
1606
1607 void Mesh::writeMED(const char* pMEDfilename, const char* pMeshName)
1608 {
1609     bool noMesh = true;
1610     MULTIPR_LOG("Write MED: " << pMEDfilename << endl);
1611     if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
1612     if (strlen(pMEDfilename) == 0) throw IllegalArgumentException("pMEDfilename size is 0", __FILE__, __LINE__);
1613
1614     if (pMeshName == NULL) throw NullArgumentException("pMeshName should not be NULL", __FILE__, __LINE__);
1615     if (strlen(pMeshName) == 0) throw IllegalArgumentException("pMeshName size is 0", __FILE__, __LINE__);
1616
1617     remove(pMEDfilename);
1618
1619     char aMeshName[MED_TAILLE_NOM + 1];
1620     strncpy(aMeshName, pMeshName, MED_TAILLE_NOM);
1621
1622     //---------------------------------------------------------------------
1623     // Create the new MED file (WRITE_ONLY)
1624     //---------------------------------------------------------------------
1625     med_idt newMEDfile = MEDouvrir(const_cast<char*>(pMEDfilename), MED_CREATION);
1626     if (newMEDfile == -1) throw IOException("", __FILE__, __LINE__);
1627
1628     //---------------------------------------------------------------------
1629     // Write scalars
1630     //---------------------------------------------------------------------
1631     // no scalars to write
1632
1633     //---------------------------------------------------------------------
1634     // Create mesh: must be created first
1635     //---------------------------------------------------------------------
1636     med_err ret = MEDmaaCr(
1637         newMEDfile,
1638         aMeshName,
1639         mMeshDim,
1640         MED_NON_STRUCTURE,
1641         mMeshDesc);
1642
1643     if (ret != 0) throw IOException("", __FILE__, __LINE__);
1644     MULTIPR_LOG("    Create mesh: |" << aMeshName << "|: OK" << endl);
1645
1646     //---------------------------------------------------------------------
1647     // Write nodes and elements (mesh must exist)
1648     //---------------------------------------------------------------------
1649     if (mNodes == NULL) throw IllegalStateException("mNodes should not be NULL", __FILE__, __LINE__);
1650     mNodes->writeMED(newMEDfile, aMeshName);
1651     MULTIPR_LOG("    Write nodes: ok" << endl);
1652
1653         for (int i = 0; i < eMaxMedMesh; ++i)
1654         {
1655                 if (mElements[i] != NULL)
1656                 {
1657                         noMesh = false;
1658                         mElements[i]->writeMED(newMEDfile, aMeshName, mMeshDim);
1659                 }
1660         }
1661         if (noMesh == true)
1662         {
1663                 //throw IllegalStateException("No mesh writen.", __FILE__, __LINE__);
1664                 return ;
1665         }
1666
1667         MULTIPR_LOG("    write elt: ok" << endl);
1668
1669     //---------------------------------------------------------------------
1670     // Write families (mesh must exist)
1671     //---------------------------------------------------------------------
1672     // MED assume there is a family 0 in the file.
1673     bool    createFamZero = true;
1674     for(std::vector<Family*>::iterator it =  mFamilies.begin();
1675         it != mFamilies.end(); ++it)
1676     {
1677         if ((*it)->getId() == 0)
1678         {
1679             createFamZero = false;
1680             break;
1681         }
1682     }
1683     if (createFamZero)
1684     {
1685         Family  famZero;
1686         famZero.setId(0);
1687         strcpy(const_cast<char*>(famZero.getName()), "FAMILLE_ZERO");
1688         famZero.writeMED(newMEDfile, aMeshName);
1689     }
1690
1691     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; ++itFam)
1692     {
1693         Family* fam = mFamilies[itFam];
1694                 fam->writeMED(newMEDfile, aMeshName);
1695     }
1696     MULTIPR_LOG("    Write families: ok" << endl);
1697
1698     //---------------------------------------------------------------------
1699     // Write profil
1700     //---------------------------------------------------------------------
1701     for (unsigned itProf = 0; itProf < mProfils.size(); ++itProf)
1702     {
1703         Profil* prof = mProfils[itProf];
1704         prof->writeMED(newMEDfile);
1705     }
1706
1707     //---------------------------------------------------------------------
1708     // Write Gauss localization (must be written before fields)
1709     //---------------------------------------------------------------------
1710     for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ;
1711         itGaussLoc++)
1712     {
1713         GaussLoc* gaussLoc = mGaussLoc[itGaussLoc];
1714         gaussLoc->writeMED(newMEDfile);
1715     }
1716
1717     MULTIPR_LOG("    Write Gauss: ok" << endl);
1718
1719     //---------------------------------------------------------------------
1720     // Write fields
1721     //---------------------------------------------------------------------
1722     std::set<std::string> fieldNameList;
1723     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
1724     {
1725         Field* field = mFields[itField];
1726         if (fieldNameList.find(std::string(field->getName())) != fieldNameList.end())
1727         {
1728             field->writeMED(newMEDfile, aMeshName, false);
1729         }
1730         else
1731         {
1732             fieldNameList.insert(std::string(field->getName()));
1733             field->writeMED(newMEDfile, aMeshName);
1734         }
1735     }
1736
1737     MULTIPR_LOG("    Write fields: ok" << endl);
1738
1739     //---------------------------------------------------------------------
1740     // Close the new MED file
1741     //---------------------------------------------------------------------
1742     ret = MEDfermer(newMEDfile);
1743     if (ret != 0) throw IOException("", __FILE__, __LINE__);
1744 }
1745
1746 void Mesh::readGaussLoc()
1747 {
1748     MULTIPR_LOG("Gauss ref: ");
1749     med_int numGauss = MEDnGauss(mMEDfile);
1750     if (numGauss < 0) throw IOException("", __FILE__, __LINE__);
1751     MULTIPR_LOG(numGauss << ": OK\n");
1752
1753     for (int itGauss = 1 ; itGauss <= numGauss ; itGauss++)
1754     {
1755         GaussLoc* gaussLoc = new GaussLoc();
1756         gaussLoc->readMED(mMEDfile, itGauss);
1757         MULTIPR_LOG((*gaussLoc) << endl);
1758         mGaussLoc.push_back(gaussLoc);
1759         mGaussLocNameToGaussLoc.insert(make_pair(gaussLoc->getName(), gaussLoc));
1760     }
1761 }
1762
1763 void Mesh::readFamilies()
1764 {
1765     med_int numFamilies = MEDnFam(mMEDfile, mMeshName);
1766     if (numFamilies <= 0) throw IOException("", __FILE__, __LINE__);
1767
1768     for (int itFam = 1 ; itFam <= numFamilies ; ++itFam)
1769     {
1770         Family* fam = new Family();
1771         fam->readMED(mMEDfile, mMeshName, itFam);
1772         mFamilies.push_back(fam);
1773     }
1774 }
1775
1776
1777 void Mesh::finalizeFamiliesAndGroups()
1778 {
1779     if (mFamilies.size() == 0)
1780     {
1781         return ;
1782     }
1783     //---------------------------------------------------------------------
1784     // Build mapping between family id and pointers towards families
1785     //---------------------------------------------------------------------
1786     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; ++itFam)
1787     {
1788         Family* fam  = mFamilies[itFam];
1789         mFamIdToFam.insert(make_pair(fam->getId(), fam));
1790     }
1791     //---------------------------------------------------------------------
1792     // Fill families of nodes
1793     //---------------------------------------------------------------------
1794     for (int itNode = 1 ; itNode <= mNodes->getNumberOfNodes() ; ++itNode)
1795     {
1796         // get family of the ith nodes
1797         int famIdent = mNodes->getFamIdent(itNode - 1); // MED nodes start at 1
1798
1799                 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1800
1801         if (itFam == mFamIdToFam.end())
1802         {
1803             char msg[256];
1804             sprintf(msg, "wrong family of nodes for node #%d: family %d not found", itNode, famIdent);
1805             throw IllegalStateException(msg, __FILE__, __LINE__);
1806         }
1807
1808         Family* fam = (*itFam).second;
1809
1810         // add the current node to its family
1811         fam->insertElt(itNode);
1812         fam->setIsFamilyOfNodes(true);
1813     }
1814     //---------------------------------------------------------------------
1815     // Fill families of elements
1816     //---------------------------------------------------------------------
1817         for (int itMesh = 0; itMesh < eMaxMedMesh; ++itMesh)
1818         {
1819                 if (mElements[itMesh] != NULL)
1820                 {
1821                 for (int itElt = 1 ; itElt <= mElements[itMesh]->getNumberOfElements() ; itElt++)
1822                 {
1823                         // get family of the ith element (MED index start at 1)
1824                         int famIdent = mElements[itMesh]->getFamilyIdentifier(itElt - 1);
1825
1826                                 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1827
1828                         if (itFam == mFamIdToFam.end())
1829                         {
1830                         char msg[256];
1831                         sprintf(msg, "wrong family of elements for element #%d: family %d not found", itElt, famIdent);
1832                         throw IllegalStateException(msg, __FILE__, __LINE__);
1833                         }
1834
1835                         Family* fam = (*itFam).second;
1836
1837                         // add the current element to its family
1838                         fam->insertElt( itElt, (eMeshType)itMesh);
1839                         fam->setIsFamilyOfNodes(false);
1840                 }
1841                 }
1842         }
1843     //---------------------------------------------------------------------
1844     // Build groups
1845     //---------------------------------------------------------------------
1846     // for each family
1847     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1848     {
1849         mFamilies[itFam]->buildGroups(mGroups, mGroupNameToGroup);
1850     }
1851 }
1852
1853
1854 void Mesh::readFields()
1855 {
1856     //---------------------------------------------------------------------
1857     // Read number of fields
1858     //---------------------------------------------------------------------
1859     MULTIPR_LOG("Read fields: ");
1860     med_int numFields = MEDnChamp(mMEDfile, 0);
1861     if (numFields <= 0) throw IOException("", __FILE__, __LINE__);
1862     MULTIPR_LOG(numFields << ": OK\n");
1863
1864     //---------------------------------------------------------------------
1865     // Iterate over fields
1866     //---------------------------------------------------------------------
1867     // for each field, read number of components and others infos
1868         for (int itField = 1 ; itField <= numFields ; itField++)
1869     {
1870                 for (int i = 0; i < eMaxMedMesh; ++i)
1871                 {
1872                 Field* field = new Field();
1873                 field->readMED(mMEDfile, itField, mMeshName, CELL_TYPES[i]);
1874
1875             // if the nth field does not apply on our mesh => slip it
1876                 if (field->isEmpty())
1877                 {
1878                     delete field;
1879                 }
1880                 else
1881                 {
1882                                 mFields.push_back(field);
1883                                 // ReadMed will always work with fields on node so we need to stop the first time.
1884                                 // ie : CELL_TYPES[i] is not used in this case.
1885                                 if (field->isFieldOnNodes())
1886                                 {
1887                                         break;
1888                                 }
1889                 }
1890                 }
1891     }
1892 }
1893
1894
1895 ostream& operator<<(ostream& pOs, Mesh& pM)
1896 {
1897     pOs << "Mesh: " << endl;
1898     pOs << "    MED file =|" << pM.mMEDfilename << "|" << endl;
1899     pOs << "    Name     =|" << pM.mMeshName << "|" << endl;
1900     pOs << "    Unv name =|" << pM.mMeshUName << "|" << endl;
1901     pOs << "    Desc     =|" << pM.mMeshDesc << "|" << endl;
1902     pOs << "    Dim      =" << pM.mMeshDim << endl;
1903     pOs << "    Type     =" << ((pM.mMeshType == MED_STRUCTURE)?"STRUCTURE":"NON_STRUCTURE") << endl;
1904     pOs << "    BBox     =[" << pM.mMeshBBoxMin[0] << " ; " << pM.mMeshBBoxMax[0] << "] x [" << pM.mMeshBBoxMin[1] << " ; " << pM.mMeshBBoxMax[1] << "] x [" << pM.mMeshBBoxMin[2] << " ; " << pM.mMeshBBoxMax[2] << "]" << endl;
1905
1906     int numFamOfNodes = 0;
1907     for (unsigned i = 0 ; i < pM.mFamilies.size() ; i++)
1908     {
1909         if (pM.mFamilies[i]->isFamilyOfNodes())
1910         {
1911             numFamOfNodes++;
1912         }
1913     }
1914
1915     int numGroupsOfNodes = 0;
1916     for (unsigned i = 0 ; i < pM.mGroups.size() ; i++)
1917     {
1918         if (pM.mGroups[i]->isGroupOfNodes())
1919         {
1920             numGroupsOfNodes++;
1921         }
1922     }
1923
1924     if (pM.mFlagPrintAll)
1925     {
1926         cout << (*(pM.mNodes)) << endl;
1927         for (int i = 0; i < eMaxMedMesh; ++i)
1928                 {
1929                         if (pM.mElements[i] != NULL)
1930                         {
1931                         cout << (*(pM.mElements[i])) << endl;
1932                         }
1933                 }
1934
1935         pOs << "    Families : #=" << pM.mFamilies.size() << " (nodes=" << numFamOfNodes << " ; elements=" << (pM.mFamilies.size() - numFamOfNodes) << ")" << endl;
1936         for (unsigned i = 0 ; i < pM.mFamilies.size() ; i++)
1937         {
1938             cout << (*(pM.mFamilies[i])) << endl;
1939         }
1940
1941         pOs << "    Groups   : #=" << pM.mGroups.size() << " (nodes=" << numGroupsOfNodes << " ; elements=" << (pM.mGroups.size() - numGroupsOfNodes) << ")" << endl;
1942         for (unsigned i = 0 ; i < pM.mGroups.size() ; i++)
1943         {
1944             cout << (*(pM.mGroups[i])) << endl;
1945         }
1946
1947         pOs << "    Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1948         for (unsigned i = 0 ; i < pM.mGaussLoc.size() ; i++)
1949         {
1950             cout << (*(pM.mGaussLoc[i])) << endl;
1951         }
1952
1953         pOs << "    Fields   : #=" << pM.mFields.size() << endl;
1954         for (unsigned i = 0 ; i < pM.mFields.size() ; i++)
1955         {
1956             cout << (*(pM.mFields[i])) << endl;
1957         }
1958     }
1959     else
1960     {
1961         if (pM.mNodes != NULL)
1962         {
1963             pOs << "    Nodes    : #=" << pM.mNodes->getNumberOfNodes() << endl;
1964         }
1965         for (int i = 0; i < eMaxMedMesh; ++i)
1966                 {
1967                         if (pM.mElements[i] != NULL)
1968                         {
1969                                 const set<med_int>& setOfNodes = pM.mElements[i]->getSetOfNodes();
1970                                 if (setOfNodes.size() == 0)
1971                                 {
1972                                         pOs << "    Elt      : #=" << pM.mElements[i]->getNumberOfElements() << endl;
1973                                 }
1974                                 else
1975                                 {
1976                                         set<med_int>::iterator itNode = setOfNodes.end();
1977                                         itNode--;
1978                                         pOs << "    Elt      : #=" << pM.mElements[i]->getNumberOfElements() << " node_id_min=" << (*(setOfNodes.begin())) << " node_id_max=" << (*itNode) << endl;
1979                                 }
1980                         }
1981         }
1982         pOs << "    Families : #=" << pM.mFamilies.size() << " (nodes=" << numFamOfNodes << " ; elements=" << (pM.mFamilies.size() - numFamOfNodes) << ")" << endl;
1983         pOs << "    Groups   : #=" << pM.mGroups.size() << " (nodes=" << numGroupsOfNodes << " ; elements=" << (pM.mGroups.size() - numGroupsOfNodes) << ")" << endl;
1984         pOs << "    Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1985         pOs << "    Fields   : #=" << pM.mFields.size() << endl;
1986     }
1987
1988     return pOs;
1989 }
1990
1991
1992 } // namespace  multipr
1993
1994 // EOF