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 (int 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 (int 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
 471         // get prefix from the original MED filename
 472         string strPrefix = removeExtension(mMEDfilename, ".med");
 473
 474         int numGroup = 1;
 475
 476         // for each group
 477         for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
 478         {
 479                 Group* currentGroup = mGroups[itGroup];
 480
 481                 // skip this group if it is a group of nodes
 482                 if (currentGroup->isGroupOfNodes())
 483                 {
 484                         continue;
 485                 }
 486
 487                 char strPartName[256];
 488                 sprintf(strPartName, "%s_%d", mMeshName, numGroup);
 489
 490                 char strMEDfilename[256];
 491                 sprintf(strMEDfilename, "%s_grain%d.med", strPrefix.c_str(), numGroup);
 492
 493                 Mesh* mesh = createFromGroup(currentGroup, mMeshName);
 494
 495                 // skip the group which contain all the others groups, even it contains only 1 group
 496                 if ((mesh->mElements->getNumberOfElements() == mElements->getNumberOfElements()) && (mGroups.size() > 1))
 497                 {
 498                         delete mesh;
 499                         continue;
 500                 }
 501
 502                 meshDis->addMesh(
 503                         MeshDisPart::MULTIPR_WRITE_MESH,
 504                         mMeshName,
 505                         numGroup,
 506                         strPartName,
 507                         "localhost",
 508                         strMEDfilename,
 509                         mesh);
 510
 511                 numGroup++;
 512         }
 513
 514         return meshDis;
 515 }
 516
 517
 518 Mesh* Mesh::decimate(
 519         const char* pFilterName,
 520         const char* pArgv,
 521         const char* pNameNewMesh)
 522 {
 523         //---------------------------------------------------------------------
 524         // Check parameters
 525         //---------------------------------------------------------------------
 526         if (pFilterName == NULL) throw NullArgumentException("pFilterName should not be NULL", __FILE__, __LINE__);
 527         if (pArgv == NULL) throw NullArgumentException("pArgv should not be NULL", __FILE__, __LINE__);
 528         if (pNameNewMesh == NULL) throw NullArgumentException("pNameNewMesh should not be NULL", __FILE__, __LINE__);
 529
 530         //---------------------------------------------------------------------
 531         // Instanciate filter used for decimation
 532         //---------------------------------------------------------------------
 533         DecimationFilter* filter = DecimationFilter::create(pFilterName);
 534
 535         //---------------------------------------------------------------------
 536         // Create new mesh by decimating current one
 537         //---------------------------------------------------------------------
 538         Mesh* decimatedMesh = filter->apply(this, pArgv, pNameNewMesh);
 539
 540         //---------------------------------------------------------------------
 541         // Cleans
 542         //---------------------------------------------------------------------
 543         delete filter;
 544
 545         return decimatedMesh;
 546 }
 547
 548
 549
 550 void Mesh::getAllPointsOfField(Field* pField, int pTimeStepIt, std::vector<PointOfField>& pPoints)
 551 {
 552         //---------------------------------------------------------------------
 553         // Check arguments
 554         //---------------------------------------------------------------------
 555
 556         if (pField == NULL) throw NullArgumentException("field should not be NULL", __FILE__, __LINE__);
 557         if (pTimeStepIt < 1) throw IllegalArgumentException("invalid field iteration; should be >= 1", __FILE__, __LINE__);
 558
 559         if (mMeshDim != 3) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 560         if (pField->getType() != MED_FLOAT64) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 561         if (pField->getNumberOfComponents() != 1) throw UnsupportedOperationException("field have more than 1 component (vectorial field, expected scalar field)", __FILE__, __LINE__);
 562
 563         //---------------------------------------------------------------------
 564         // Collect points
 565         //---------------------------------------------------------------------
 566
 567         if (pField->isFieldOnNodes())
 568         {
 569                 //-------------------------------------------------------------
 570                 // Case 1: field of nodes
 571                 //-------------------------------------------------------------
 572                 if (mNodes == NULL) throw IllegalStateException("no nodes in the current mesh", __FILE__, __LINE__);
 573
 574                 // for each node
 575                 for (int itNode = 0, size = mNodes->getNumberOfNodes() ; itNode < size ; itNode++)
 576                 {
 577                         // collect coordinates and value of the point
 578                         const med_float* coo = mNodes->getCoordinates(itNode);
 579
 580                         const med_float* val =
 581                                 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itNode + 1));
 582
 583                         // add new point
 584                         pPoints.push_back(PointOfField(coo[0], coo[1], coo[2], val[0]));
 585                 }
 586         }
 587         else
 588         {
 589                 //-------------------------------------------------------------
 590                 // Case 2: field of elements
 591                 //-------------------------------------------------------------
 592
 593                 if (mElements == NULL) throw IllegalStateException("no elements in the current mesh", __FILE__, __LINE__);
 594                 if (mElements->getTypeOfPrimitives() != MED_TETRA10) throw UnsupportedOperationException("only support TETRA10 mesh", __FILE__, __LINE__);
 595
 596                 const string& nameGaussLoc = pField->getNameGaussLoc(pTimeStepIt);
 597                 GaussLoc* gaussLoc = getGaussLocByName(nameGaussLoc.c_str());
 598                 if (gaussLoc == NULL) throw IllegalStateException("no Gauss localization for these elements", __FILE__, __LINE__);
 599
 600                 int numGauss = pField->getNumberOfGaussPointsByElement(pTimeStepIt);
 601
 602                 int size = gaussLoc->getDim() * gaussLoc->getNumGaussPoints();
 603                 med_float* cooGaussPts = new med_float[size];
 604
 605                 int dim = mElements->getTypeOfPrimitives() / 100;
 606                 int numNodes = mElements->getTypeOfPrimitives() % 100;
 607                 size = dim * numNodes;
 608                 med_float* cooNodes = new med_float[size];
 609
 610                 // for each elements
 611                 for (int itElt = 0, size = mElements->getNumberOfElements() ; itElt < size ; itElt++)
 612                 {
 613                         // get coordinates of nodes of the current elements
 614                         // OPTIMIZATION: ASSUME TETRA10: ONLY GETS THE 4 FIRST NODES OF EACH ELEMENT
 615                         mElements->getCoordinates(itElt, mNodes, cooNodes, 4);
 616
 617                         // compute coordinates of gauss points
 618                         gaussLoc->getCoordGaussPoints(cooNodes, cooGaussPts);
 619
 620                         //printArray2D(cooGaussPts, 5, 3, "Gauss pt");
 621
 622                         const med_float* val =
 623                                 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itElt + 1));
 624
 625                         // for each point of Gauss of the element
 626                         med_float* srcCoo = cooGaussPts;
 627                         for (int itPtGauss = 0 ; itPtGauss < numGauss ; itPtGauss++)
 628                         {
 629                                 pPoints.push_back(PointOfField(srcCoo[0], srcCoo[1], srcCoo[2], val[itPtGauss]));
 630                                 srcCoo += 3;
 631                         }
 632                 }
 633
 634                 delete[] cooNodes;
 635                 delete[] cooGaussPts;
 636         }
 637 }
 638
 639
 640 float Mesh::evalDefaultRadius(int pN) const
 641 {
 642         if (mFields.size() == 0) return 1.0f;
 643
 644         //---------------------------------------------------------------------
 645         // Compute default radius
 646         //---------------------------------------------------------------------
 647
 648         med_float volumeBBox =
 649                 (mMeshBBoxMax[0] - mMeshBBoxMin[0]) *
 650                 (mMeshBBoxMax[1] - mMeshBBoxMin[1]) *
 651                 (mMeshBBoxMax[2] - mMeshBBoxMin[2]);
 652
 653         if (isnan(volumeBBox))
 654         {
 655                 return 1.0f;
 656         }
 657
 658         const med_float k = 0.8; // considered 80% of the volume
 659
 660         // get nunmber of gauss points in the field
 661         try
 662         {
 663                 Field* anyField = mFields[mFields.size()-1];
 664                 int numTimeSteps = anyField->getNumberOfTimeSteps();
 665
 666                 int numGaussPoints = getNumberOfElements() * anyField->getNumberOfGaussPointsByElement(numTimeSteps-1);
 667
 668                 med_float radius = med_float(pow( (3.0/4.0) * pN * k * volumeBBox / (3.1415 * numGaussPoints), 1.0/3.0));
 669
 670                 return float(radius);
 671         }
 672         catch (...)
 673         {
 674                 return 1.0f;
 675         }
 676 }
 677
 678
 679 med_geometrie_element CELL_TYPES[MED_NBR_GEOMETRIE_MAILLE] =
 680 {
 681         MED_POINT1,
 682         MED_SEG2,
 683         MED_SEG3,
 684         MED_TRIA3,
 685         MED_TRIA6,
 686         MED_QUAD4,
 687         MED_QUAD8,
 688         MED_TETRA4,
 689         MED_TETRA10,
 690         MED_HEXA8,
 691         MED_HEXA20,
 692         MED_PENTA6,
 693         MED_PENTA15
 694 };
 695
 696
 697 char CELL_NAMES[MED_NBR_GEOMETRIE_MAILLE][MED_TAILLE_NOM + 1] =
 698 {
 699         "MED_POINT1",
 700         "MED_SEG2",
 701         "MED_SEG3",
 702         "MED_TRIA3",
 703         "MED_TRIA6",
 704         "MED_QUAD4",
 705         "MED_QUAD8",
 706         "MED_TETRA4",
 707         "MED_TETRA10",
 708         "MED_HEXA8",
 709         "MED_HEXA20",
 710         "MED_PENTA6",
 711         "MED_PENTA15",
 712         "MED_PYRA5",
 713         "MED_PYRA13"
 714 };
 715
 716
 717 void Mesh::readSequentialMED(const char* pMEDfilename, const char* pMeshName)
 718 {
 719         //cout << "File: |" << multipr::getFilenameWithoutPath(pMEDfilename) << "|" << endl;
 720         //cout << "Path: |" << multipr::getPath(pMEDfilename) << "|" << endl;
 721
 722         reset();
 723
 724         //---------------------------------------------------------------------
 725         // Check arguments
 726         //---------------------------------------------------------------------
 727         MULTIPR_LOG("Check arguments: ");
 728         if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
 729         if (pMeshName == NULL) throw NullArgumentException("pMeshName should not be NULL", __FILE__, __LINE__);
 730         MULTIPR_LOG("OK\n");
 731
 732         strncpy(mMEDfilename, pMEDfilename, 256);
 733         strncpy(mMeshName, pMeshName, MED_TAILLE_NOM);
 734
 735         //---------------------------------------------------------------------
 736         // Open MED file (READ_ONLY)
 737         //---------------------------------------------------------------------
 738         MULTIPR_LOG("Open MED file: ");
 739         mMEDfile = MEDouvrir(mMEDfilename, MED_LECTURE); // open MED file for reading
 740         if (mMEDfile <= 0) throw FileNotFoundException("MED file not found", __FILE__, __LINE__);
 741         MULTIPR_LOG("OK\n");
 742
 743         //---------------------------------------------------------------------
 744         // Check valid HDF format
 745         //---------------------------------------------------------------------
 746         MULTIPR_LOG("Format HDF: ");
 747         if (MEDformatConforme(mMEDfilename) != 0) throw IOException("invalid file", __FILE__, __LINE__);
 748         MULTIPR_LOG("OK\n");
 749
 750         //---------------------------------------------------------------------
 751         // Get MED version
 752         //---------------------------------------------------------------------
 753         MULTIPR_LOG("MED version: ");
 754         med_int verMajor, verMinor, verRelease;
 755         med_err ret = MEDversionLire(mMEDfile, &verMajor, &verMinor, &verRelease);
 756         if (ret != 0) throw IOException("error while reading MED version", __FILE__, __LINE__);
 757         MULTIPR_LOG(verMajor << "." << verMinor << "." << verRelease << ": OK\n");
 758
 759         //---------------------------------------------------------------------
 760         // Check that there is no profil
 761         //---------------------------------------------------------------------
 762         MULTIPR_LOG("#profils must be 0: ");
 763         med_int nbProfils = MEDnProfil(mMEDfile);
 764         if (nbProfils != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 765         MULTIPR_LOG("OK\n");
 766
 767         //---------------------------------------------------------------------
 768         // Read all Gauss localizations
 769         //---------------------------------------------------------------------
 770         readGaussLoc();
 771
 772         //---------------------------------------------------------------------
 773         // Read scalars (should be 0)
 774         //---------------------------------------------------------------------
 775         MULTIPR_LOG("Scalars: ");
 776         med_int nbScalars = MEDnScalaire(mMEDfile);
 777         if (nbScalars != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 778         MULTIPR_LOG(nbScalars << ": OK\n");
 779
 780         //---------------------------------------------------------------------
 781         // Find the mesh
 782         //---------------------------------------------------------------------
 783         // read number of meshes
 784         MULTIPR_LOG("Num meshes: ");
 785         med_int nbMeshes = MEDnMaa(mMEDfile);
 786         if (nbMeshes <= 0) throw IOException("i/o error while reading number of meshes in MED file", __FILE__, __LINE__);
 787         MULTIPR_LOG(nbMeshes << ": OK\n");
 788
 789         med_int meshIndex = -1;
 790         // iteration over mesh to find the mesh we want
 791         // for each mesh in the file (warning: first mesh is number 1)
 792         for (int itMesh = 1 ; itMesh <= nbMeshes ; itMesh++)
 793         {
 794                 char meshName[MED_TAILLE_NOM + 1];
 795
 796                 ret = MEDmaaInfo(
 797                         mMEDfile,
 798                         itMesh,
 799                         meshName,
 800                         &mMeshDim,
 801                         &mMeshType,
 802                         mMeshDesc);
 803
 804                 if (ret != 0) throw IOException("i/o error while reading mesh information in MED file", __FILE__, __LINE__);
 805                 MULTIPR_LOG("Mesh: |" << meshName << "|");
 806
 807                 // test if the current mesh is the mesh we want
 808                 if (strcmp(pMeshName, meshName) == 0)
 809                 {
 810                         // *** mesh found ***
 811                         MULTIPR_LOG(" OK (found)\n");
 812                         meshIndex = itMesh;
 813                         break;
 814                 }
 815                 else
 816                 {
 817                         // not the mesh we want: skip this mesh
 818                         MULTIPR_LOG(" skipped\n");
 819                 }
 820         }
 821
 822         if (meshIndex == -1)
 823         {
 824                 throw IllegalStateException("mesh not found in the given MED file", __FILE__, __LINE__);
 825         }
 826
 827         //---------------------------------------------------------------------
 828         // Check mesh validity
 829         //---------------------------------------------------------------------
 830         // dimension of the mesh must be 3 (= 3D mesh)
 831         MULTIPR_LOG("Mesh is 3D: ");
 832         if (mMeshDim != 3) throw UnsupportedOperationException("dimension of the mesh should be 3; other dimension not yet implemented", __FILE__, __LINE__);
 833         MULTIPR_LOG("OK\n");
 834
 835         // mesh must not be a grid
 836         MULTIPR_LOG("Mesh is not a grid: ");
 837         if (mMeshType != MED_NON_STRUCTURE)
 838                 throw UnsupportedOperationException("grid not supported", __FILE__, __LINE__);
 839         MULTIPR_LOG("OK\n");
 840
 841         // mesh must only contain TETRA10 elements
 842         MULTIPR_LOG("Only TETRA10: ");
 843         med_connectivite connectivite = MED_NOD; // NODAL CONNECTIVITY ONLY
 844         bool onlyTETRA10 = true;
 845         int numTetra10 = -1;
 846         for (int itCell = 0 ; itCell < MED_NBR_GEOMETRIE_MAILLE ; itCell++)
 847         {
 848                 med_int meshNumCells = MEDnEntMaa(
 849                         mMEDfile,
 850                         mMeshName,
 851                         MED_CONN,
 852                         MED_MAILLE,
 853                         CELL_TYPES[itCell],
 854                         connectivite);
 855
 856                 if ((meshNumCells > 0) && (strcmp(CELL_NAMES[itCell], "MED_TETRA10") != 0))
 857                 {
 858                         onlyTETRA10 = false;
 859                         break;
 860                 }
 861                 if (strcmp(CELL_NAMES[itCell], "MED_TETRA10") == 0)
 862                 {
 863                         numTetra10 = meshNumCells;
 864                 }
 865         }
 866
 867         if (!onlyTETRA10) throw UnsupportedOperationException("mesh should only contain TETRA10 elements", __FILE__, __LINE__);
 868         MULTIPR_LOG(numTetra10 << ": OK\n");
 869
 870         // everything is OK...
 871
 872         //---------------------------------------------------------------------
 873         // Check num joint = 0
 874         //---------------------------------------------------------------------
 875         MULTIPR_LOG("Num joints: ");
 876         med_int numJoints = MEDnJoint(mMEDfile, mMeshName);
 877         MULTIPR_LOG(numJoints << ": OK\n");
 878
 879         //---------------------------------------------------------------------
 880         // Check num equivalence = 0
 881         //---------------------------------------------------------------------
 882         MULTIPR_LOG("Num equivalences: ");
 883         med_int numEquiv = MEDnEquiv(mMEDfile, mMeshName);
 884         MULTIPR_LOG(numEquiv << ": OK\n");
 885
 886         //---------------------------------------------------------------------
 887         // Read nodes
 888         //---------------------------------------------------------------------
 889         mNodes = new Nodes();
 890         mNodes->readMED(mMEDfile, mMeshName, mMeshDim);
 891         mNodes->getBBox(mMeshBBoxMin, mMeshBBoxMax);
 892
 893         //---------------------------------------------------------------------
 894         // Read elements
 895         //---------------------------------------------------------------------
 896         mElements = new Elements();
 897         mElements->readMED(mMEDfile, mMeshName, mMeshDim, MED_MAILLE, MED_TETRA10);
 898
 899         //---------------------------------------------------------------------
 900         // Read families
 901         //---------------------------------------------------------------------
 902         readFamilies();
 903         finalizeFamiliesAndGroups();
 904
 905         //---------------------------------------------------------------------
 906         // Read fields
 907         //---------------------------------------------------------------------
 908         readFields();
 909
 910         //---------------------------------------------------------------------
 911         // Close the MED file
 912         //---------------------------------------------------------------------
 913         MULTIPR_LOG("Close MED file: ");
 914         ret = MEDfermer(mMEDfile);
 915         if (ret != 0) throw IOException("i/o error while closing MED file", __FILE__, __LINE__);
 916         MULTIPR_LOG("OK\n");
 917 }
 918
 919
 920 void Mesh::writeMED(const char* pMEDfilename)
 921 {
 922         MULTIPR_LOG("Write MED: " << pMEDfilename << endl);
 923
 924         if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
 925         if (strlen(pMEDfilename) == 0) throw IllegalArgumentException("pMEDfilename size is 0", __FILE__, __LINE__);
 926
 927         remove(pMEDfilename);
 928
 929         //---------------------------------------------------------------------
 930         // Create the new MED file (WRITE_ONLY)
 931         //---------------------------------------------------------------------
 932         med_idt newMEDfile = MEDouvrir(const_cast<char*>(pMEDfilename), MED_CREATION);
 933         if (newMEDfile == -1) throw IOException("", __FILE__, __LINE__);
 934
 935         //---------------------------------------------------------------------
 936         // Write scalars
 937         //---------------------------------------------------------------------
 938         // no scalars to write
 939
 940         //---------------------------------------------------------------------
 941         // Create mesh: must be created first
 942         //---------------------------------------------------------------------
 943         med_err ret = MEDmaaCr(
 944                 newMEDfile,
 945                 mMeshName,
 946                 mMeshDim,
 947                 MED_NON_STRUCTURE,
 948                 mMeshDesc);
 949
 950         if (ret != 0) throw IOException("", __FILE__, __LINE__);
 951         MULTIPR_LOG("    Create mesh: |" << mMeshName << "|: OK" << endl);
 952
 953         //---------------------------------------------------------------------
 954         // Write nodes and elements (mesh must exist)
 955         //---------------------------------------------------------------------
 956         mNodes->writeMED(newMEDfile, mMeshName);
 957         MULTIPR_LOG("    Write nodes: ok" << endl);
 958         mElements->writeMED(newMEDfile, mMeshName, mMeshDim);
 959         MULTIPR_LOG("    write elt: ok" << endl);
 960
 961         //---------------------------------------------------------------------
 962         // Write families (mesh must exist)
 963         //---------------------------------------------------------------------
 964         for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
 965         {
 966                 Family* fam = mFamilies[itFam];
 967                 fam->writeMED(newMEDfile, mMeshName);
 968         }
 969         MULTIPR_LOG("    Write families: ok" << endl);
 970
 971         //---------------------------------------------------------------------
 972         // Write profil
 973         //---------------------------------------------------------------------
 974         // no profil
 975
 976         //---------------------------------------------------------------------
 977         // Write Gauss localization (must be written before fields)
 978         //---------------------------------------------------------------------
 979         for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 980         {
 981
 982                 GaussLoc* gaussLoc = mGaussLoc[itGaussLoc];
 983                 gaussLoc->writeMED(newMEDfile);
 984         }
 985         MULTIPR_LOG("    Write Gauss: ok" << endl);
 986
 987         //---------------------------------------------------------------------
 988         // Write fields
 989         //---------------------------------------------------------------------
 990         for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
 991         {
 992                 Field* field = mFields[itField];
 993                 field->writeMED(newMEDfile, mMeshName);
 994         }
 995         MULTIPR_LOG("    Write fields: ok" << endl);
 996
 997         //---------------------------------------------------------------------
 998         // Close the new MED file
 999         //---------------------------------------------------------------------
1000         ret = MEDfermer(newMEDfile);
1001         if (ret != 0) throw IOException("", __FILE__, __LINE__);
1002 }
1003
1004
1005 void Mesh::readGaussLoc()
1006 {
1007         MULTIPR_LOG("Gauss ref: ");
1008         med_int numGauss = MEDnGauss(mMEDfile);
1009         if (numGauss < 0) throw IOException("", __FILE__, __LINE__);
1010         MULTIPR_LOG(numGauss << ": OK\n");
1011
1012         for (int itGauss = 1 ; itGauss <= numGauss ; itGauss++)
1013         {
1014                 GaussLoc* gaussLoc = new GaussLoc();
1015                 gaussLoc->readMED(mMEDfile, itGauss);
1016
1017                 MULTIPR_LOG((*gaussLoc) << endl);
1018
1019                 mGaussLoc.push_back(gaussLoc);
1020                 mGaussLocNameToGaussLoc.insert(make_pair(gaussLoc->getName(), gaussLoc));
1021         }
1022 }
1023
1024
1025 void Mesh::readFamilies()
1026 {
1027         med_int numFamilies = MEDnFam(mMEDfile, mMeshName);
1028         if (numFamilies <= 0) throw IOException("", __FILE__, __LINE__);
1029
1030         for (int itFam = 1 ; itFam <= numFamilies ; itFam++)
1031         {
1032                 Family* fam = new Family();
1033                 fam->readMED(mMEDfile, mMeshName, itFam);
1034                 mFamilies.push_back(fam);
1035         }
1036 }
1037
1038
1039 void Mesh::finalizeFamiliesAndGroups()
1040 {
1041         //---------------------------------------------------------------------
1042         // Build mapping between family id and pointers towards families
1043         //---------------------------------------------------------------------
1044         for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1045         {
1046                 Family* fam  = mFamilies[itFam];
1047                 mFamIdToFam.insert(make_pair(fam->getId(), fam));
1048         }
1049
1050         //---------------------------------------------------------------------
1051         // Fill families of nodes
1052         //---------------------------------------------------------------------
1053         for (int itNode = 1 ; itNode <= mNodes->getNumberOfNodes() ; itNode++)
1054         {
1055                 // get family of the ith nodes
1056                 int famIdent = mNodes->getFamIdent(itNode - 1); // MED nodes start at 1
1057                 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1058
1059                 if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1060
1061                 Family* fam = (*itFam).second;
1062
1063                 // insert the current node to its family
1064                 fam->insertElt(itNode);
1065                 fam->setIsFamilyOfNodes(true);
1066         }
1067
1068         //---------------------------------------------------------------------
1069         // Fill families of elements
1070         //---------------------------------------------------------------------
1071         for (int itElt = 1 ; itElt <= mElements->getNumberOfElements() ; itElt++)
1072         {
1073                 // get family of the ith element (MED index start at 1)
1074                 int famIdent = mElements->getFamilyIdentifier(itElt - 1);
1075                 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1076
1077                 if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1078
1079                 Family* fam = (*itFam).second;
1080
1081                 // insert the current node its family
1082                 fam->insertElt(itElt);
1083                 fam->setIsFamilyOfNodes(false);
1084         }
1085
1086         //---------------------------------------------------------------------
1087         // Build groups
1088         //---------------------------------------------------------------------
1089         // for each family
1090         for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1091         {
1092                 mFamilies[itFam]->buildGroups(mGroups, mGroupNameToGroup);
1093         }
1094 }
1095
1096
1097 void Mesh::readFields()
1098 {
1099         //---------------------------------------------------------------------
1100         // Read number of fields
1101         //---------------------------------------------------------------------
1102         MULTIPR_LOG("Read fields: ");
1103         med_int numFields = MEDnChamp(mMEDfile, 0);
1104         if (numFields <= 0) throw IOException("", __FILE__, __LINE__);
1105         MULTIPR_LOG(numFields << ": OK\n");
1106
1107         //---------------------------------------------------------------------
1108         // Iterate over fields
1109         //---------------------------------------------------------------------
1110         // for each field, read number of components and others infos
1111         for (int itField = 1 ; itField <= numFields ; itField++)
1112         {
1113                 Field* field = new Field();
1114                 field->readMED(mMEDfile, itField, mMeshName);
1115
1116                 // if the nth field does not apply on our mesh => slip it
1117                 if (field->isEmpty())
1118                 {
1119                         delete field;
1120                 }
1121                 else
1122                 {
1123                         mFields.push_back(field);
1124                 }
1125         }
1126 }
1127
1128
1129 ostream& operator<<(ostream& pOs, Mesh& pM)
1130 {
1131         pOs << "Mesh: " << endl;
1132         pOs << "    MED file =|" << pM.mMEDfilename << "|" << endl;
1133         pOs << "    Name     =|" << pM.mMeshName << "|" << endl;
1134         pOs << "    Unv name =|" << pM.mMeshUName << "|" << endl;
1135         pOs << "    Desc     =|" << pM.mMeshDesc << "|" << endl;
1136         pOs << "    Dim      =" << pM.mMeshDim << endl;
1137         pOs << "    Type     =" << ((pM.mMeshType == MED_STRUCTURE)?"STRUCTURE":"NON_STRUCTURE") << endl;
1138         pOs << "    BBox     =[" << pM.mMeshBBoxMin[0] << " ; " << pM.mMeshBBoxMax[0] << "] x [" << pM.mMeshBBoxMin[1] << " ; " << pM.mMeshBBoxMax[1] << "] x [" << pM.mMeshBBoxMin[2] << " ; " << pM.mMeshBBoxMax[2] << "]" << endl;
1139
1140         if (pM.mFlagPrintAll)
1141         {
1142                 cout << (*(pM.mNodes)) << endl;
1143                 cout << (*(pM.mElements)) << endl;
1144
1145                 pOs << "    Families : #=" << pM.mFamilies.size() << endl;
1146                 for (unsigned i = 0 ; i < pM.mFamilies.size() ; i++)
1147                 {
1148                         cout << (*(pM.mFamilies[i])) << endl;
1149                 }
1150
1151                 pOs << "    Groups   : #=" << pM.mGroups.size() << endl;
1152                 for (unsigned i = 0 ; i < pM.mGroups.size() ; i++)
1153                 {
1154                         cout << (*(pM.mGroups[i])) << endl;
1155                 }
1156
1157                 pOs << "    Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1158                 for (unsigned i = 0 ; i < pM.mGaussLoc.size() ; i++)
1159                 {
1160                         cout << (*(pM.mGaussLoc[i])) << endl;
1161                 }
1162
1163                 pOs << "    Fields   : #=" << pM.mFields.size() << endl;
1164                 for (unsigned i = 0 ; i < pM.mFields.size() ; i++)
1165                 {
1166                         cout << (*(pM.mFields[i])) << endl;
1167                 }
1168         }
1169         else
1170         {
1171                 pOs << "    Nodes    : #=" << pM.mNodes->getNumberOfNodes() << endl;
1172
1173                 const set<med_int>& setOfNodes = pM.mElements->getSetOfNodes();
1174                 if (setOfNodes.size() == 0)
1175                 {
1176                         pOs << "    Elt      : #=" << pM.mElements->getNumberOfElements() << endl;
1177                 }
1178                 else
1179                 {
1180                         set<med_int>::iterator itNode = setOfNodes.end();
1181                         itNode--;
1182                         pOs << "    Elt      : #=" << pM.mElements->getNumberOfElements() << " node_id_min=" << (*(setOfNodes.begin())) << " node_id_max=" << (*itNode) << endl;
1183                 }
1184
1185                 pOs << "    Families : #=" << pM.mFamilies.size() << endl;
1186                 pOs << "    Groups   : #=" << pM.mGroups.size() << endl;
1187                 pOs << "    Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1188                 pOs << "    Fields   : #=" << pM.mFields.size() << endl;
1189         }
1190
1191         return pOs;
1192 }
1193
1194
1195 } // namespace  multipr
1196
1197 // EOF