src/MULTIPR/MULTIPR_Mesh.cxx

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