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
  33 using namespace std;
  34
  35
  36 namespace multipr
  37 {
  38
  39
  40 //*****************************************************************************
  41 // Class Mesh implementation
  42 //*****************************************************************************
  43
  44 Mesh::Mesh()
  45 {
  46     mNodes    = NULL;
  47     mElements = NULL;
  48
  49     reset();
  50 }
  51
  52
  53 Mesh::~Mesh()
  54 {
  55     reset();
  56 }
  57
  58
  59 void Mesh::reset()
  60 {
  61     mMEDfilename[0] = '\0';
  62     mMEDfile        = 0;
  63
  64     mMeshName[0]    = '\0';
  65     mMeshUName[0]   = '\0';
  66     mMeshDesc[0]    = '\0';
  67     mMeshDim        = -1;
  68     mMeshType       = MED_NON_STRUCTURE;
  69
  70     for (int itDim = 0 ; itDim < 3 ; itDim++)
  71     {
  72         mMeshBBoxMin[itDim] = numeric_limits<med_float>::quiet_NaN();
  73         mMeshBBoxMax[itDim] = numeric_limits<med_float>::quiet_NaN();
  74     }
  75
  76     if (mNodes != NULL)    { delete mNodes;    mNodes = NULL; }
  77     if (mElements != NULL) { delete mElements; mElements = NULL; }
  78
  79     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
  80     {
  81         delete mFamilies[itFam];
  82     }
  83     mFamilies.clear();
  84     mFamIdToFam.clear();
  85
  86     for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
  87     {
  88         delete mGroups[itGroup];
  89     }
  90     mGroups.clear();
  91     mGroupNameToGroup.clear();
  92
  93     for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
  94     {
  95         delete mGaussLoc[itGaussLoc];
  96     }
  97     mGaussLoc.clear();
  98     mGaussLocNameToGaussLoc.clear();
  99
 100     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 101     {
 102         delete mFields[itField];
 103     }
 104     mFields.clear();
 105
 106     for (unsigned itProfil = 0 ; itProfil < mProfils.size() ; itProfil++)
 107     {
 108         delete mProfils[itProfil];
 109     }
 110     mProfils.clear();
 111
 112     mFlagPrintAll = false;
 113 }
 114
 115
 116 vector<string> Mesh::getNameScalarFields() const
 117 {
 118     vector<string> res;
 119
 120     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 121     {
 122         Field* currentField = mFields[itField];
 123
 124         // only get scalar fields, not vectorial fields
 125         // (because, currently, decimation can only be performed on scalar fields)
 126         if (currentField->getNumberOfComponents() == 1)
 127         {
 128             res.push_back(currentField->getName());
 129         }
 130     }
 131
 132     return res;
 133 }
 134
 135
 136 int Mesh::getTimeStamps(const char* pFieldName) const
 137 {
 138     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 139     {
 140         Field* currentField = mFields[itField];
 141         if (strcmp(currentField->getName(), pFieldName) == 0)
 142         {
 143             return currentField->getNumberOfTimeSteps();
 144         }
 145     }
 146
 147     return 0;
 148 }
 149
 150
 151 Field* Mesh::getFieldByName(const char* pFieldName) const
 152 {
 153     if (pFieldName == NULL) throw NullArgumentException("pFieldName should not be NULL", __FILE__, __LINE__);
 154
 155     Field* retField = NULL;
 156
 157     // for each field
 158     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 159     {
 160         Field* currentField = mFields[itField];
 161         if (strcmp(pFieldName, currentField->getName()) == 0)
 162         {
 163             // field found!
 164             retField = currentField;
 165             break;
 166         }
 167     }
 168
 169     return retField;
 170 }
 171
 172
 173 GaussLoc* Mesh::getGaussLocByName(const char* pGaussLocName) const
 174 {
 175     if (pGaussLocName == NULL) throw NullArgumentException("pGaussLocName should not be NULL", __FILE__, __LINE__);
 176
 177     map<string, GaussLoc*>::const_iterator itGaussLoc = mGaussLocNameToGaussLoc.find(pGaussLocName);
 178     GaussLoc* retGaussLoc = NULL;
 179
 180     if (itGaussLoc != mGaussLocNameToGaussLoc.end())
 181     {
 182         retGaussLoc = (*itGaussLoc).second;
 183     }
 184
 185     return retGaussLoc;
 186 }
 187
 188
 189 int Mesh::getNumberOfElements() const
 190 {
 191     if (mElements == NULL) throw IllegalStateException("", __FILE__, __LINE__);
 192
 193     return mElements->getNumberOfElements();
 194 }
 195
 196
 197 Mesh* Mesh::createFromSetOfElements(const std::set<med_int>& pSetOfElements, const char* pNewMeshName)
 198 {
 199     if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName", __FILE__, __LINE__);
 200
 201     //---------------------------------------------------------------------
 202     // Create a new mesh
 203     //---------------------------------------------------------------------
 204     Mesh* mesh = new Mesh();
 205
 206     //---------------------------------------------------------------------
 207     // Build name of the new mesh
 208     //---------------------------------------------------------------------
 209     strcpy(mesh->mMeshName, pNewMeshName);
 210
 211     MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
 212
 213     //---------------------------------------------------------------------
 214     // Fill general infos
 215     //---------------------------------------------------------------------
 216     strcpy(mesh->mMeshUName, mMeshUName);
 217     strcpy(mesh->mMeshDesc, mMeshDesc);
 218
 219     mesh->mMeshDim = mMeshDim;
 220     mesh->mMeshType = mMeshType;
 221
 222     MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
 223     MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
 224     MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
 225     MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
 226
 227     //---------------------------------------------------------------------
 228     // Build nodes and elements
 229     //---------------------------------------------------------------------
 230     // get all elements involved
 231     mesh->mElements = mElements->extractSubSet(pSetOfElements);
 232     MULTIPR_LOG((*(mesh->mElements)) << endl);
 233
 234     // get all nodes involved
 235     const set<med_int> setOfNodes = mesh->mElements->getSetOfNodes();
 236     mesh->mNodes = mNodes->extractSubSet(setOfNodes);
 237     MULTIPR_LOG((*(mesh->mNodes)) << endl);
 238
 239     //---------------------------------------------------------------------
 240     // Remap nodes
 241     //---------------------------------------------------------------------
 242     mesh->mElements->remap();
 243     MULTIPR_LOG((*(mesh->mElements)) << endl);
 244
 245     //---------------------------------------------------------------------
 246     // Build families
 247     //---------------------------------------------------------------------
 248     MULTIPR_LOG("Build fam.:" << endl);
 249     // get families of nodes
 250     {
 251         set<med_int> famOfNodes = mesh->mNodes->getSetOfFamilies();
 252         for (set<med_int>::iterator itFam = famOfNodes.begin() ; itFam != famOfNodes.end() ; itFam++)
 253         {
 254             Family* famSrc = mFamIdToFam[*itFam];
 255             cout << (*famSrc) << endl;
 256             Family* famDest = famSrc->extractGroup(NULL);
 257             mesh->mFamilies.push_back(famDest);
 258         }
 259     }
 260
 261     // get families of elements
 262     {
 263         set<med_int> famOfElt = mesh->mElements->getSetOfFamilies();
 264         for (set<med_int>::iterator itFam = famOfElt.begin() ; itFam != famOfElt.end() ; itFam++)
 265         {
 266             Family* famSrc = mFamIdToFam[*itFam];
 267             Family* famDest = famSrc->extractGroup(NULL);
 268             mesh->mFamilies.push_back(famDest);
 269         }
 270     }
 271
 272     MULTIPR_LOG("Finalize:");
 273
 274     // fill families with elements and build groups
 275     mesh->finalizeFamiliesAndGroups();
 276
 277     MULTIPR_LOG("OK\n");
 278
 279     //---------------------------------------------------------------------
 280     // Create new fields and collect Gauss
 281     //---------------------------------------------------------------------
 282     // for each field
 283     set<string> newSetOfGauss;
 284     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 285     {
 286         Field* currentField = mFields[itField];
 287
 288         Field* newField;
 289         if (currentField->isFieldOnNodes())
 290         {
 291             newField = currentField->extractSubSet(setOfNodes);
 292         }
 293         else
 294         {
 295             newField = currentField->extractSubSet(pSetOfElements);
 296         }
 297
 298         if (!newField->isEmpty())
 299         {
 300             mesh->mFields.push_back(newField);
 301             newField->getSetOfGaussLoc(newSetOfGauss);
 302         }
 303     }
 304     MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
 305
 306     //---------------------------------------------------------------------
 307     // Build Gauss infos
 308     //---------------------------------------------------------------------
 309     for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
 310     {
 311         const string& keyName = (*itSet);
 312
 313         GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
 314         if (gaussLoc != NULL)
 315         {
 316             GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
 317             mesh->mGaussLoc.push_back(copyGaussLoc);
 318             mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
 319         }
 320     }
 321
 322     //---------------------------------------------------------------------
 323     // Compute bbox
 324     //---------------------------------------------------------------------
 325     mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
 326
 327     return mesh;
 328 }
 329
 330
 331 Mesh* Mesh::createFromGroup(const Group* pGroup, const char* pNewMeshName)
 332 {
 333     if (pGroup == NULL) throw NullArgumentException("pGroup should not be NULL", __FILE__, __LINE__);
 334     if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName should not be NULL", __FILE__, __LINE__);
 335     if (strlen(pNewMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("pNewMeshName length too long", __FILE__, __LINE__);
 336
 337     //---------------------------------------------------------------------
 338     // Create a new mesh
 339     //---------------------------------------------------------------------
 340     Mesh* mesh = new Mesh();
 341
 342     //---------------------------------------------------------------------
 343     // Build name of the new mesh
 344     //---------------------------------------------------------------------
 345     strcpy(mesh->mMeshName, pNewMeshName);
 346
 347     MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
 348
 349     //---------------------------------------------------------------------
 350     // Fill general infos
 351     //---------------------------------------------------------------------
 352     strcpy(mesh->mMeshUName, mMeshUName);
 353     strcpy(mesh->mMeshDesc, mMeshDesc);
 354
 355     mesh->mMeshDim = mMeshDim;
 356     mesh->mMeshType = mMeshType;
 357
 358     MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
 359     MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
 360     MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
 361     MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
 362
 363     //---------------------------------------------------------------------
 364     // Build nodes and elements
 365     //---------------------------------------------------------------------
 366     // get all elements involved
 367     const set<med_int> setOfElt = pGroup->getSetOfElt();
 368     mesh->mElements = mElements->extractSubSet(setOfElt);
 369     MULTIPR_LOG((*(mesh->mElements)) << endl);
 370
 371     // get all nodes involved
 372     const set<med_int> setOfNodes = mesh->mElements->getSetOfNodes();
 373     mesh->mNodes = mNodes->extractSubSet(setOfNodes);
 374     MULTIPR_LOG((*(mesh->mNodes)) << endl);
 375
 376     //---------------------------------------------------------------------
 377     // Remap nodes
 378     //---------------------------------------------------------------------
 379     mesh->mElements->remap();
 380     MULTIPR_LOG((*(mesh->mElements)) << endl);
 381
 382     //---------------------------------------------------------------------
 383     // Build families
 384     //---------------------------------------------------------------------
 385     MULTIPR_LOG("Build fam.:" << endl);
 386     // get families of nodes
 387     {
 388         set<med_int> famOfNodes = mesh->mNodes->getSetOfFamilies();
 389         for (set<med_int>::iterator itFam = famOfNodes.begin() ; itFam != famOfNodes.end() ; itFam++)
 390         {
 391             Family* famSrc = mFamIdToFam[*itFam];
 392             Family* famDest = famSrc->extractGroup(pGroup->getName().c_str());
 393             mesh->mFamilies.push_back(famDest);
 394         }
 395     }
 396
 397     // get families of elements
 398     {
 399         set<med_int> famOfElt = mesh->mElements->getSetOfFamilies();
 400         for (set<med_int>::iterator itFam = famOfElt.begin() ; itFam != famOfElt.end() ; itFam++)
 401         {
 402             Family* famSrc = mFamIdToFam[*itFam];
 403             Family* famDest = famSrc->extractGroup(pGroup->getName().c_str());
 404             mesh->mFamilies.push_back(famDest);
 405         }
 406     }
 407
 408     MULTIPR_LOG("Finalize:");
 409
 410     // fill families with elements and build groups
 411     mesh->finalizeFamiliesAndGroups();
 412
 413     MULTIPR_LOG("OK\n");
 414
 415     //---------------------------------------------------------------------
 416     // Create new fields and collect Gauss
 417     //---------------------------------------------------------------------
 418     // for each field
 419     set<string> newSetOfGauss;
 420     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 421     {
 422         Field* currentField = mFields[itField];
 423
 424         Field* newField;
 425         if (currentField->isFieldOnNodes())
 426         {
 427             newField = currentField->extractSubSet(setOfNodes);
 428         }
 429         else
 430         {
 431             newField = currentField->extractSubSet(setOfElt);
 432         }
 433
 434         if (!newField->isEmpty())
 435         {
 436             mesh->mFields.push_back(newField);
 437             newField->getSetOfGaussLoc(newSetOfGauss);
 438         }
 439     }
 440     MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
 441
 442     //---------------------------------------------------------------------
 443     // Build Gauss infos
 444     //---------------------------------------------------------------------
 445     for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
 446     {
 447         const string& keyName = (*itSet);
 448
 449         GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
 450         if (gaussLoc != NULL)
 451         {
 452             GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
 453             mesh->mGaussLoc.push_back(copyGaussLoc);
 454             mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
 455         }
 456     }
 457
 458     //---------------------------------------------------------------------
 459     // Compute bbox
 460     //---------------------------------------------------------------------
 461     mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
 462
 463     return mesh;
 464 }
 465
 466
 467 MeshDis* Mesh::splitGroupsOfElements()
 468 {
 469     MeshDis* meshDis = new MeshDis();
 470     meshDis->setSequentialMEDFilename(mMEDfilename);
 471
 472     // get prefix from the original MED filename
 473     string strPrefix = removeExtension(mMEDfilename, ".med");
 474
 475     int numGroup = 1;
 476
 477     // for each group
 478     for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
 479     {
 480         Group* currentGroup = mGroups[itGroup];
 481
 482         // skip this group if it is a group of nodes
 483         if (currentGroup->isGroupOfNodes())
 484         {
 485             continue;
 486         }
 487
 488         char strPartName[256];
 489         sprintf(strPartName, "%s_%d", mMeshName, numGroup);
 490
 491         char strMEDfilename[256];
 492         sprintf(strMEDfilename, "%s_grain%d.med", strPrefix.c_str(), numGroup);
 493
 494         Mesh* mesh = createFromGroup(currentGroup, mMeshName);
 495
 496         // skip the group which contain all the others groups, even it contains only 1 group
 497         if ((mesh->mElements->getNumberOfElements() == mElements->getNumberOfElements()) && (mGroups.size() > 1))
 498         {
 499             delete mesh;
 500             continue;
 501         }
 502
 503         meshDis->addMesh(
 504             MeshDisPart::MULTIPR_WRITE_MESH,
 505             mMeshName,
 506             numGroup,
 507             strPartName,
 508             "localhost",
 509             strMEDfilename,
 510             mesh);
 511
 512         numGroup++;
 513     }
 514
 515     return meshDis;
 516 }
 517
 518
 519 Mesh* Mesh::decimate(
 520     const char* pFilterName,
 521     const char* pArgv,
 522     const char* pNameNewMesh)
 523 {
 524     //---------------------------------------------------------------------
 525     // Check parameters
 526     //---------------------------------------------------------------------
 527     if (pFilterName == NULL) throw NullArgumentException("pFilterName should not be NULL", __FILE__, __LINE__);
 528     if (pArgv == NULL) throw NullArgumentException("pArgv should not be NULL", __FILE__, __LINE__);
 529     if (pNameNewMesh == NULL) throw NullArgumentException("pNameNewMesh should not be NULL", __FILE__, __LINE__);
 530
 531     //---------------------------------------------------------------------
 532     // Instanciate filter used for decimation
 533     //---------------------------------------------------------------------
 534     DecimationFilter* filter = DecimationFilter::create(pFilterName);
 535
 536     //---------------------------------------------------------------------
 537     // Create new mesh by decimating current one
 538     //---------------------------------------------------------------------
 539     Mesh* decimatedMesh = filter->apply(this, pArgv, pNameNewMesh);
 540
 541     //---------------------------------------------------------------------
 542     // Cleans
 543     //---------------------------------------------------------------------
 544     delete filter;
 545
 546     return decimatedMesh;
 547 }
 548
 549
 550
 551 void Mesh::getAllPointsOfField(Field* pField, int pTimeStepIt, std::vector<PointOfField>& pPoints)
 552 {
 553     //---------------------------------------------------------------------
 554     // Check arguments
 555     //---------------------------------------------------------------------
 556
 557     if (pField == NULL) throw NullArgumentException("field should not be NULL", __FILE__, __LINE__);
 558     if (pTimeStepIt < 1) throw IllegalArgumentException("invalid field iteration; should be >= 1", __FILE__, __LINE__);
 559
 560     if (mMeshDim != 3) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 561     if (pField->getType() != MED_FLOAT64) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 562     if (pField->getNumberOfComponents() != 1) throw UnsupportedOperationException("field have more than 1 component (vectorial field, expected scalar field)", __FILE__, __LINE__);
 563
 564     //---------------------------------------------------------------------
 565     // Collect points
 566     //---------------------------------------------------------------------
 567
 568     if (pField->isFieldOnNodes())
 569     {
 570         //-------------------------------------------------------------
 571         // Case 1: field of nodes
 572         //-------------------------------------------------------------
 573         if (mNodes == NULL) throw IllegalStateException("no nodes in the current mesh", __FILE__, __LINE__);
 574
 575         // for each node
 576         for (int itNode = 0, size = mNodes->getNumberOfNodes() ; itNode < size ; itNode++)
 577         {
 578             // collect coordinates and value of the point
 579             const med_float* coo = mNodes->getCoordinates(itNode);
 580
 581             const med_float* val =
 582                 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itNode + 1));
 583
 584             // add new point
 585             pPoints.push_back(PointOfField(coo[0], coo[1], coo[2], val[0]));
 586         }
 587     }
 588     else
 589     {
 590         //-------------------------------------------------------------
 591         // Case 2: field of elements
 592         //-------------------------------------------------------------
 593
 594         if (mElements == NULL) throw IllegalStateException("no elements in the current mesh", __FILE__, __LINE__);
 595         if (mElements->getTypeOfPrimitives() != MED_TETRA10) throw UnsupportedOperationException("only support TETRA10 mesh", __FILE__, __LINE__);
 596
 597         const string& nameGaussLoc = pField->getNameGaussLoc(pTimeStepIt);
 598         GaussLoc* gaussLoc = getGaussLocByName(nameGaussLoc.c_str());
 599         if (gaussLoc == NULL) throw IllegalStateException("no Gauss localization for these elements", __FILE__, __LINE__);
 600
 601         int numGauss = pField->getNumberOfGaussPointsByElement(pTimeStepIt);
 602
 603         int size = gaussLoc->getDim() * gaussLoc->getNumGaussPoints();
 604         med_float* cooGaussPts = new med_float[size];
 605
 606         int dim = mElements->getTypeOfPrimitives() / 100;
 607         int numNodes = mElements->getTypeOfPrimitives() % 100;
 608         size = dim * numNodes;
 609         med_float* cooNodes = new med_float[size];
 610
 611         // for each elements
 612         for (int itElt = 0, size = mElements->getNumberOfElements() ; itElt < size ; itElt++)
 613         {
 614             // get coordinates of nodes of the current elements
 615             // OPTIMIZATION: ASSUME TETRA10: ONLY GETS THE 4 FIRST NODES OF EACH ELEMENT
 616             mElements->getCoordinates(itElt, mNodes, cooNodes, 4);
 617
 618             // compute coordinates of gauss points
 619             gaussLoc->getCoordGaussPoints(cooNodes, cooGaussPts);
 620
 621             //printArray2D(cooGaussPts, 5, 3, "Gauss pt"); // debug
 622
 623             const med_float* val =
 624                 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itElt + 1));
 625
 626             // for each point of Gauss of the element
 627             med_float* srcCoo = cooGaussPts;
 628             for (int itPtGauss = 0 ; itPtGauss < numGauss ; itPtGauss++)
 629             {
 630                 pPoints.push_back(PointOfField(srcCoo[0], srcCoo[1], srcCoo[2], val[itPtGauss]));
 631                 srcCoo += 3;
 632             }
 633         }
 634
 635         delete[] cooNodes;
 636         delete[] cooGaussPts;
 637     }
 638 }
 639
 640
 641 float Mesh::evalDefaultRadius(int pN) const
 642 {
 643     if (mFields.size() == 0) return 1.0f;
 644
 645     //---------------------------------------------------------------------
 646     // Compute default radius
 647     //---------------------------------------------------------------------
 648
 649     med_float volumeBBox =
 650         (mMeshBBoxMax[0] - mMeshBBoxMin[0]) *
 651         (mMeshBBoxMax[1] - mMeshBBoxMin[1]) *
 652         (mMeshBBoxMax[2] - mMeshBBoxMin[2]);
 653
 654     if (isnan(volumeBBox))
 655     {
 656         return 1.0f;
 657     }
 658
 659     const med_float k = 0.8; // considered 80% of the volume
 660
 661     // get nunmber of gauss points in the field
 662     try
 663     {
 664         Field* anyField = mFields[mFields.size()-1];
 665         int numTimeSteps = anyField->getNumberOfTimeSteps();
 666
 667         int numGaussPoints = getNumberOfElements() * anyField->getNumberOfGaussPointsByElement(numTimeSteps-1);
 668
 669         med_float radius = med_float(pow( (3.0/4.0) * pN * k * volumeBBox / (3.1415 * numGaussPoints), 1.0/3.0));
 670
 671         return float(radius);
 672     }
 673     catch (...)
 674     {
 675         return 1.0f;
 676     }
 677 }
 678
 679
 680 med_geometrie_element CELL_TYPES[MED_NBR_GEOMETRIE_MAILLE] =
 681 {
 682     MED_POINT1,
 683     MED_SEG2,
 684     MED_SEG3,
 685     MED_TRIA3,
 686     MED_TRIA6,
 687     MED_QUAD4,
 688     MED_QUAD8,
 689     MED_TETRA4,
 690     MED_TETRA10,
 691     MED_HEXA8,
 692     MED_HEXA20,
 693     MED_PENTA6,
 694     MED_PENTA15
 695 };
 696
 697
 698 char CELL_NAMES[MED_NBR_GEOMETRIE_MAILLE][MED_TAILLE_NOM + 1] =
 699 {
 700     "MED_POINT1",
 701     "MED_SEG2",
 702     "MED_SEG3",
 703     "MED_TRIA3",
 704     "MED_TRIA6",
 705     "MED_QUAD4",
 706     "MED_QUAD8",
 707     "MED_TETRA4",
 708     "MED_TETRA10",
 709     "MED_HEXA8",
 710     "MED_HEXA20",
 711     "MED_PENTA6",
 712     "MED_PENTA15",
 713     "MED_PYRA5",
 714     "MED_PYRA13"
 715 };
 716
 717
 718 void Mesh::readSequentialMED(const char* pMEDfilename, const char* pMeshName)
 719 {
 720     reset();
 721
 722     //---------------------------------------------------------------------
 723     // Check arguments
 724     //---------------------------------------------------------------------
 725     MULTIPR_LOG("Check arguments: ");
 726     if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
 727     if (pMeshName == NULL) throw NullArgumentException("pMeshName should not be NULL", __FILE__, __LINE__);
 728     MULTIPR_LOG("OK\n");
 729
 730     strncpy(mMEDfilename, pMEDfilename, 256);
 731     strncpy(mMeshName, pMeshName, MED_TAILLE_NOM);
 732
 733     //---------------------------------------------------------------------
 734     // Open MED file (READ_ONLY)
 735     //---------------------------------------------------------------------
 736     MULTIPR_LOG("Open MED file: ");
 737     mMEDfile = MEDouvrir(mMEDfilename, MED_LECTURE); // open MED file for reading
 738     if (mMEDfile <= 0) throw FileNotFoundException("MED file not found", __FILE__, __LINE__);
 739     MULTIPR_LOG("OK\n");
 740
 741     //---------------------------------------------------------------------
 742     // Check valid HDF format
 743     //---------------------------------------------------------------------
 744     MULTIPR_LOG("Format HDF: ");
 745     if (MEDformatConforme(mMEDfilename) != 0) throw IOException("invalid file", __FILE__, __LINE__);
 746     MULTIPR_LOG("OK\n");
 747
 748     //---------------------------------------------------------------------
 749     // Get MED version
 750     //---------------------------------------------------------------------
 751     MULTIPR_LOG("MED version: ");
 752     med_int verMajor, verMinor, verRelease;
 753     med_err ret = MEDversionLire(mMEDfile, &verMajor, &verMinor, &verRelease);
 754     if (ret != 0) throw IOException("error while reading MED version", __FILE__, __LINE__);
 755     MULTIPR_LOG(verMajor << "." << verMinor << "." << verRelease << ": OK\n");
 756
 757     //---------------------------------------------------------------------
 758     // Check that there is no profil
 759     //---------------------------------------------------------------------
 760     MULTIPR_LOG("#profils must be 0: ");
 761     med_int nbProfils = MEDnProfil(mMEDfile);
 762     if (nbProfils != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 763     MULTIPR_LOG("OK\n");
 764
 765     //---------------------------------------------------------------------
 766     // Read all Gauss localizations
 767     //---------------------------------------------------------------------
 768     readGaussLoc();
 769
 770     //---------------------------------------------------------------------
 771     // Read scalars (should be 0)
 772     //---------------------------------------------------------------------
 773     MULTIPR_LOG("Scalars: ");
 774     med_int nbScalars = MEDnScalaire(mMEDfile);
 775     if (nbScalars != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 776     MULTIPR_LOG(nbScalars << ": OK\n");
 777
 778     //---------------------------------------------------------------------
 779     // Find the mesh
 780     //---------------------------------------------------------------------
 781     // read number of meshes
 782     MULTIPR_LOG("Num meshes: ");
 783     med_int nbMeshes = MEDnMaa(mMEDfile);
 784     if (nbMeshes <= 0) throw IOException("i/o error while reading number of meshes in MED file", __FILE__, __LINE__);
 785     MULTIPR_LOG(nbMeshes << ": OK\n");
 786
 787     med_int meshIndex = -1;
 788     // iteration over mesh to find the mesh we want
 789     // for each mesh in the file (warning: first mesh is number 1)
 790     for (int itMesh = 1 ; itMesh <= nbMeshes ; itMesh++)
 791     {
 792         char meshName[MED_TAILLE_NOM + 1];
 793
 794         ret = MEDmaaInfo(
 795             mMEDfile,
 796             itMesh,
 797             meshName,
 798             &mMeshDim,
 799             &mMeshType,
 800             mMeshDesc);
 801
 802         if (ret != 0) throw IOException("i/o error while reading mesh information in MED file", __FILE__, __LINE__);
 803         MULTIPR_LOG("Mesh: |" << meshName << "|");
 804
 805         // test if the current mesh is the mesh we want
 806         if (strcmp(pMeshName, meshName) == 0)
 807         {
 808             // *** mesh found ***
 809             MULTIPR_LOG(" OK (found)\n");
 810             meshIndex = itMesh;
 811             break;
 812         }
 813         else
 814         {
 815             // not the mesh we want: skip this mesh
 816             MULTIPR_LOG(" skipped\n");
 817         }
 818     }
 819
 820     if (meshIndex == -1)
 821     {
 822         throw IllegalStateException("mesh not found in the given MED file", __FILE__, __LINE__);
 823     }
 824
 825     //---------------------------------------------------------------------
 826     // Check mesh validity
 827     //---------------------------------------------------------------------
 828     // dimension of the mesh must be 3 (= 3D mesh)
 829     MULTIPR_LOG("Mesh is 3D: ");
 830     if (mMeshDim != 3) throw UnsupportedOperationException("dimension of the mesh should be 3; other dimension not yet implemented", __FILE__, __LINE__);
 831     MULTIPR_LOG("OK\n");
 832
 833     // mesh must not be a grid
 834     MULTIPR_LOG("Mesh is not a grid: ");
 835     if (mMeshType != MED_NON_STRUCTURE)
 836         throw UnsupportedOperationException("grid not supported", __FILE__, __LINE__);
 837     MULTIPR_LOG("OK\n");
 838
 839     // mesh must only contain TETRA10 elements
 840     MULTIPR_LOG("Only TETRA10: ");
 841     med_connectivite connectivite = MED_NOD; // NODAL CONNECTIVITY ONLY
 842     bool onlyTETRA10 = true;
 843     int numTetra10 = -1;
 844     for (int itCell = 0 ; itCell < MED_NBR_GEOMETRIE_MAILLE ; itCell++)
 845     {
 846         med_int meshNumCells = MEDnEntMaa(
 847             mMEDfile,
 848             mMeshName,
 849             MED_CONN,
 850             MED_MAILLE,
 851             CELL_TYPES[itCell],
 852             connectivite);
 853
 854         if ((meshNumCells > 0) && (strcmp(CELL_NAMES[itCell], "MED_TETRA10") != 0))
 855         {
 856             onlyTETRA10 = false;
 857             break;
 858         }
 859         if (strcmp(CELL_NAMES[itCell], "MED_TETRA10") == 0)
 860         {
 861             numTetra10 = meshNumCells;
 862         }
 863     }
 864
 865     if (!onlyTETRA10) throw UnsupportedOperationException("mesh should only contain TETRA10 elements", __FILE__, __LINE__);
 866     MULTIPR_LOG(numTetra10 << ": OK\n");
 867
 868     // everything is OK...
 869
 870     //---------------------------------------------------------------------
 871     // Check num joint = 0
 872     //---------------------------------------------------------------------
 873     MULTIPR_LOG("Num joints: ");
 874     med_int numJoints = MEDnJoint(mMEDfile, mMeshName);
 875     MULTIPR_LOG(numJoints << ": OK\n");
 876
 877     //---------------------------------------------------------------------
 878     // Check num equivalence = 0
 879     //---------------------------------------------------------------------
 880     MULTIPR_LOG("Num equivalences: ");
 881     med_int numEquiv = MEDnEquiv(mMEDfile, mMeshName);
 882     MULTIPR_LOG(numEquiv << ": OK\n");
 883
 884     //---------------------------------------------------------------------
 885     // Read nodes
 886     //---------------------------------------------------------------------
 887     mNodes = new Nodes();
 888     mNodes->readMED(mMEDfile, mMeshName, mMeshDim);
 889     mNodes->getBBox(mMeshBBoxMin, mMeshBBoxMax);
 890
 891     //---------------------------------------------------------------------
 892     // Read elements
 893     //---------------------------------------------------------------------
 894     mElements = new Elements();
 895     mElements->readMED(mMEDfile, mMeshName, mMeshDim, MED_MAILLE, MED_TETRA10);
 896
 897     //---------------------------------------------------------------------
 898     // Read families
 899     //---------------------------------------------------------------------
 900     readFamilies();
 901     finalizeFamiliesAndGroups();
 902
 903     //---------------------------------------------------------------------
 904     // Read fields
 905     //---------------------------------------------------------------------
 906     readFields();
 907
 908     //---------------------------------------------------------------------
 909     // Close the MED file
 910     //---------------------------------------------------------------------
 911     MULTIPR_LOG("Close MED file: ");
 912     ret = MEDfermer(mMEDfile);
 913     if (ret != 0) throw IOException("i/o error while closing MED file", __FILE__, __LINE__);
 914     MULTIPR_LOG("OK\n");
 915 }
 916
 917
 918 void Mesh::writeMED(const char* pMEDfilename)
 919 {
 920     MULTIPR_LOG("Write MED: " << pMEDfilename << endl);
 921
 922     if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
 923     if (strlen(pMEDfilename) == 0) throw IllegalArgumentException("pMEDfilename size is 0", __FILE__, __LINE__);
 924
 925     remove(pMEDfilename);
 926
 927     //---------------------------------------------------------------------
 928     // Create the new MED file (WRITE_ONLY)
 929     //---------------------------------------------------------------------
 930     med_idt newMEDfile = MEDouvrir(const_cast<char*>(pMEDfilename), MED_CREATION);
 931     if (newMEDfile == -1) throw IOException("", __FILE__, __LINE__);
 932
 933     //---------------------------------------------------------------------
 934     // Write scalars
 935     //---------------------------------------------------------------------
 936     // no scalars to write
 937
 938     //---------------------------------------------------------------------
 939     // Create mesh: must be created first
 940     //---------------------------------------------------------------------
 941     med_err ret = MEDmaaCr(
 942         newMEDfile,
 943         mMeshName,
 944         mMeshDim,
 945         MED_NON_STRUCTURE,
 946         mMeshDesc);
 947
 948     if (ret != 0) throw IOException("", __FILE__, __LINE__);
 949     MULTIPR_LOG("    Create mesh: |" << mMeshName << "|: OK" << endl);
 950
 951     //---------------------------------------------------------------------
 952     // Write nodes and elements (mesh must exist)
 953     //---------------------------------------------------------------------
 954     mNodes->writeMED(newMEDfile, mMeshName);
 955     MULTIPR_LOG("    Write nodes: ok" << endl);
 956     mElements->writeMED(newMEDfile, mMeshName, mMeshDim);
 957     MULTIPR_LOG("    write elt: ok" << endl);
 958
 959     //---------------------------------------------------------------------
 960     // Write families (mesh must exist)
 961     //---------------------------------------------------------------------
 962     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
 963     {
 964         Family* fam = mFamilies[itFam];
 965         fam->writeMED(newMEDfile, mMeshName);
 966     }
 967     MULTIPR_LOG("    Write families: ok" << endl);
 968
 969     //---------------------------------------------------------------------
 970     // Write profil
 971     //---------------------------------------------------------------------
 972     // no profil
 973
 974     //---------------------------------------------------------------------
 975     // Write Gauss localization (must be written before fields)
 976     //---------------------------------------------------------------------
 977     for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 978     {
 979
 980         GaussLoc* gaussLoc = mGaussLoc[itGaussLoc];
 981         gaussLoc->writeMED(newMEDfile);
 982     }
 983     MULTIPR_LOG("    Write Gauss: ok" << endl);
 984
 985     //---------------------------------------------------------------------
 986     // Write fields
 987     //---------------------------------------------------------------------
 988     for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 989     {
 990         Field* field = mFields[itField];
 991         field->writeMED(newMEDfile, mMeshName);
 992     }
 993     MULTIPR_LOG("    Write fields: ok" << endl);
 994
 995     //---------------------------------------------------------------------
 996     // Close the new MED file
 997     //---------------------------------------------------------------------
 998     ret = MEDfermer(newMEDfile);
 999     if (ret != 0) throw IOException("", __FILE__, __LINE__);
1000 }
1001
1002
1003 void Mesh::readGaussLoc()
1004 {
1005     MULTIPR_LOG("Gauss ref: ");
1006     med_int numGauss = MEDnGauss(mMEDfile);
1007     if (numGauss < 0) throw IOException("", __FILE__, __LINE__);
1008     MULTIPR_LOG(numGauss << ": OK\n");
1009
1010     for (int itGauss = 1 ; itGauss <= numGauss ; itGauss++)
1011     {
1012         GaussLoc* gaussLoc = new GaussLoc();
1013         gaussLoc->readMED(mMEDfile, itGauss);
1014
1015         MULTIPR_LOG((*gaussLoc) << endl);
1016
1017         mGaussLoc.push_back(gaussLoc);
1018         mGaussLocNameToGaussLoc.insert(make_pair(gaussLoc->getName(), gaussLoc));
1019     }
1020 }
1021
1022
1023 void Mesh::readFamilies()
1024 {
1025     med_int numFamilies = MEDnFam(mMEDfile, mMeshName);
1026     if (numFamilies <= 0) throw IOException("", __FILE__, __LINE__);
1027
1028     for (int itFam = 1 ; itFam <= numFamilies ; itFam++)
1029     {
1030         Family* fam = new Family();
1031         fam->readMED(mMEDfile, mMeshName, itFam);
1032         mFamilies.push_back(fam);
1033     }
1034 }
1035
1036
1037 void Mesh::finalizeFamiliesAndGroups()
1038 {
1039     //---------------------------------------------------------------------
1040     // Build mapping between family id and pointers towards families
1041     //---------------------------------------------------------------------
1042     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1043     {
1044         Family* fam  = mFamilies[itFam];
1045         mFamIdToFam.insert(make_pair(fam->getId(), fam));
1046     }
1047
1048     //---------------------------------------------------------------------
1049     // Fill families of nodes
1050     //---------------------------------------------------------------------
1051     for (int itNode = 1 ; itNode <= mNodes->getNumberOfNodes() ; itNode++)
1052     {
1053         // get family of the ith nodes
1054         int famIdent = mNodes->getFamIdent(itNode - 1); // MED nodes start at 1
1055         map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1056
1057         if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1058
1059         Family* fam = (*itFam).second;
1060
1061         // insert the current node to its family
1062         fam->insertElt(itNode);
1063         fam->setIsFamilyOfNodes(true);
1064     }
1065
1066     //---------------------------------------------------------------------
1067     // Fill families of elements
1068     //---------------------------------------------------------------------
1069     for (int itElt = 1 ; itElt <= mElements->getNumberOfElements() ; itElt++)
1070     {
1071         // get family of the ith element (MED index start at 1)
1072         int famIdent = mElements->getFamilyIdentifier(itElt - 1);
1073         map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1074
1075         if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1076
1077         Family* fam = (*itFam).second;
1078
1079         // insert the current node its family
1080         fam->insertElt(itElt);
1081         fam->setIsFamilyOfNodes(false);
1082     }
1083
1084     //---------------------------------------------------------------------
1085     // Build groups
1086     //---------------------------------------------------------------------
1087     // for each family
1088     for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1089     {
1090         mFamilies[itFam]->buildGroups(mGroups, mGroupNameToGroup);
1091     }
1092 }
1093
1094
1095 void Mesh::readFields()
1096 {
1097     //---------------------------------------------------------------------
1098     // Read number of fields
1099     //---------------------------------------------------------------------
1100     MULTIPR_LOG("Read fields: ");
1101     med_int numFields = MEDnChamp(mMEDfile, 0);
1102     if (numFields <= 0) throw IOException("", __FILE__, __LINE__);
1103     MULTIPR_LOG(numFields << ": OK\n");
1104
1105     //---------------------------------------------------------------------
1106     // Iterate over fields
1107     //---------------------------------------------------------------------
1108     // for each field, read number of components and others infos
1109     for (int itField = 1 ; itField <= numFields ; itField++)
1110     {
1111         Field* field = new Field();
1112         field->readMED(mMEDfile, itField, mMeshName);
1113
1114         // if the nth field does not apply on our mesh => slip it
1115         if (field->isEmpty())
1116         {
1117             delete field;
1118         }
1119         else
1120         {
1121             mFields.push_back(field);
1122         }
1123     }
1124 }
1125
1126
1127 ostream& operator<<(ostream& pOs, Mesh& pM)
1128 {
1129     pOs << "Mesh: " << endl;
1130     pOs << "    MED file =|" << pM.mMEDfilename << "|" << endl;
1131     pOs << "    Name     =|" << pM.mMeshName << "|" << endl;
1132     pOs << "    Unv name =|" << pM.mMeshUName << "|" << endl;
1133     pOs << "    Desc     =|" << pM.mMeshDesc << "|" << endl;
1134     pOs << "    Dim      =" << pM.mMeshDim << endl;
1135     pOs << "    Type     =" << ((pM.mMeshType == MED_STRUCTURE)?"STRUCTURE":"NON_STRUCTURE") << endl;
1136     pOs << "    BBox     =[" << pM.mMeshBBoxMin[0] << " ; " << pM.mMeshBBoxMax[0] << "] x [" << pM.mMeshBBoxMin[1] << " ; " << pM.mMeshBBoxMax[1] << "] x [" << pM.mMeshBBoxMin[2] << " ; " << pM.mMeshBBoxMax[2] << "]" << endl;
1137
1138     if (pM.mFlagPrintAll)
1139     {
1140         cout << (*(pM.mNodes)) << endl;
1141         cout << (*(pM.mElements)) << endl;
1142
1143         pOs << "    Families : #=" << pM.mFamilies.size() << endl;
1144         for (unsigned i = 0 ; i < pM.mFamilies.size() ; i++)
1145         {
1146             cout << (*(pM.mFamilies[i])) << endl;
1147         }
1148
1149         pOs << "    Groups   : #=" << pM.mGroups.size() << endl;
1150         for (unsigned i = 0 ; i < pM.mGroups.size() ; i++)
1151         {
1152             cout << (*(pM.mGroups[i])) << endl;
1153         }
1154
1155         pOs << "    Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1156         for (unsigned i = 0 ; i < pM.mGaussLoc.size() ; i++)
1157         {
1158             cout << (*(pM.mGaussLoc[i])) << endl;
1159         }
1160
1161         pOs << "    Fields   : #=" << pM.mFields.size() << endl;
1162         for (unsigned i = 0 ; i < pM.mFields.size() ; i++)
1163         {
1164             cout << (*(pM.mFields[i])) << endl;
1165         }
1166     }
1167     else
1168     {
1169         pOs << "    Nodes    : #=" << pM.mNodes->getNumberOfNodes() << endl;
1170
1171         const set<med_int>& setOfNodes = pM.mElements->getSetOfNodes();
1172         if (setOfNodes.size() == 0)
1173         {
1174             pOs << "    Elt      : #=" << pM.mElements->getNumberOfElements() << endl;
1175         }
1176         else
1177         {
1178             set<med_int>::iterator itNode = setOfNodes.end();
1179             itNode--;
1180             pOs << "    Elt      : #=" << pM.mElements->getNumberOfElements() << " node_id_min=" << (*(setOfNodes.begin())) << " node_id_max=" << (*itNode) << endl;
1181         }
1182
1183         pOs << "    Families : #=" << pM.mFamilies.size() << endl;
1184         pOs << "    Groups   : #=" << pM.mGroups.size() << endl;
1185         pOs << "    Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1186         pOs << "    Fields   : #=" << pM.mFields.size() << endl;
1187     }
1188
1189     return pOs;
1190 }
1191
1192
1193 } // namespace  multipr
1194
1195 // EOF