src/MULTIPR/MULTIPR_Elements.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_Elements.cxx
   6  *
   7  * \brief   see MULTIPR_Elements.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_Elements.hxx"
  19 #include "MULTIPR_Nodes.hxx"
  20 #include "MULTIPR_Exceptions.hxx"
  21
  22 #include <iostream>
  23 #include <set>
  24 #include <map>
  25
  26 using namespace std;
  27
  28
  29 namespace multipr
  30 {
  31
  32
  33 //*****************************************************************************
  34 // Class Elements implementation
  35 //*****************************************************************************
  36
  37 Elements::Elements()
  38 {
  39     mId       = NULL;
  40     mFamIdent = NULL;
  41     mNames    = NULL;
  42     mCon      = NULL;
  43     mFlagPrintAll = false;
  44     reset();
  45 }
  46
  47 Elements::Elements(med_geometrie_element pGeomType)
  48 {
  49     mId       = NULL;
  50     mFamIdent = NULL;
  51     mNames    = NULL;
  52     mCon      = NULL;
  53     mNum      = 0;
  54     mGeom          = pGeomType;
  55     mNumNodesByElt = mGeom % 100;
  56     mDim           = mGeom / 100;
  57     mEntity        = MED_MAILLE;
  58     mFlagPrintAll  = false;
  59 }
  60
  61 Elements::~Elements()
  62 {
  63     reset();
  64 }
  65
  66
  67 void Elements::reset()
  68 {
  69     mNum           = 0;
  70     mEntity        = MED_MAILLE;
  71     mGeom          = MED_NONE;
  72     mNumNodesByElt = 0;
  73     mDim           = 0;
  74
  75     if (mId != NULL)       { delete[] mId;       mId       = NULL; }
  76     if (mFamIdent != NULL) { delete[] mFamIdent; mFamIdent = NULL; }
  77     if (mNames != NULL)    { delete[] mNames;    mNames    = NULL; }
  78     if (mCon != NULL)      { delete[] mCon;      mCon      = NULL; }
  79
  80     mSetOfNodes.clear();
  81
  82     mFlagPrintAll = false;
  83 }
  84
  85
  86 med_int Elements::getFamilyIdentifier(med_int pIndex) const
  87 {
  88     if ((pIndex < 0) || (pIndex >= mNum)) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
  89
  90     return mFamIdent[pIndex];
  91 }
  92
  93
  94 const med_int* Elements::getConnectivity(int pIndex) const
  95 {
  96     if ((pIndex < 0) || (pIndex >= mNum)) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
  97
  98     return mCon + mNumNodesByElt * pIndex;
  99 }
 100
 101
 102 void Elements::getCoordinates(med_int pIndexElt, const Nodes* pNodes, med_float* pCoo, int pFirst) const
 103 {
 104     if ((pIndexElt < 0) || (pIndexElt >= mNum)) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 105     if (pNodes == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 106     if (pCoo == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 107
 108     // get the list of nodes of the element
 109     const med_int* con = getConnectivity(pIndexElt);
 110
 111     med_float* destCoo = pCoo;
 112     int size = sizeof(med_float) * mDim;
 113
 114     // for each node of the element
 115     int n = (mNumNodesByElt < pFirst) ? mNumNodesByElt : pFirst;
 116     for (int i = 0 ; i < n ; i++)
 117     {
 118         // get index of node (MED index start at 1)
 119         med_int indexNode = con[i] - 1;
 120
 121         // get coordinates of this node
 122         const med_float* srcCoo = pNodes->getCoordinates(indexNode);
 123
 124         // copy coordinates to destCoo
 125         memcpy(destCoo, srcCoo, size);
 126
 127         // prepare next point
 128         destCoo += mDim;
 129     }
 130 }
 131
 132
 133 Elements* Elements::extractSubSet(const set<med_int>& pSetIndices) const
 134 {
 135     Elements* subset = new Elements();
 136
 137     //---------------------------------------------------------------------
 138     // Copy parameters
 139     //---------------------------------------------------------------------
 140     subset->mNum           = pSetIndices.size();
 141     subset->mEntity        = mEntity;
 142     subset->mGeom          = mGeom;          // e.g. 310 for a TETRA10
 143     subset->mNumNodesByElt = mNumNodesByElt; // e.g. 10 for a TETRA10
 144     subset->mDim           = mDim;           // e.g. 3 for a TETRA10
 145
 146     //---------------------------------------------------------------------
 147     // Allocate arrays
 148     //---------------------------------------------------------------------
 149     subset->mFamIdent     = new med_int[subset->mNum];
 150     subset->mCon          = new med_int[mNumNodesByElt * subset->mNum];
 151
 152     //---------------------------------------------------------------------
 153     // Copy subset of familys id and connectivity.
 154     //---------------------------------------------------------------------
 155     med_int* destCon = subset->mCon;
 156     set<med_int>::iterator itSet = pSetIndices.begin();
 157     for (int itElt = 0 ; itElt < subset->mNum && itSet != pSetIndices.end(); itElt++)
 158     {
 159         med_int srcIndex = (*itSet) - 1; // MED index start at 1
 160                 subset->mFamIdent[itElt] = mFamIdent[srcIndex];
 161
 162         med_int* srcCon = mCon + srcIndex * mNumNodesByElt;
 163         memcpy(destCon, srcCon, sizeof(med_int) * mNumNodesByElt);
 164
 165         destCon += mNumNodesByElt;
 166         itSet++;
 167     }
 168     //---------------------------------------------------------------------
 169     // Copy subset of identifiers if necessary
 170     //---------------------------------------------------------------------
 171     if (isIdentifiers())
 172     {
 173         itSet = pSetIndices.begin();
 174         subset->mId = new med_int[subset->mNum];
 175         for (int itElt = 0 ; itElt < subset->mNum && itSet != pSetIndices.end(); itElt++)
 176         {
 177             med_int srcIndex = (*itSet) - 1; // MED index start at 1
 178             subset->mId[itElt] = mId[srcIndex];
 179
 180             itSet++;
 181         }
 182     }
 183
 184     //---------------------------------------------------------------------
 185     // Copy subset of names if necessary
 186     //---------------------------------------------------------------------
 187     if (isNames())
 188     {
 189         itSet = pSetIndices.begin();
 190         subset->mNames = new char[MED_TAILLE_PNOM * subset->mNum + 1];
 191         char* destPtr = subset->mNames;
 192         for (int itElt = 0 ; itElt < subset->mNum && itSet != pSetIndices.end(); itElt++)
 193         {
 194             med_int srcIndex = (*itSet) - 1; // MED index start at 1
 195             char* srcPtr = mNames + srcIndex * MED_TAILLE_PNOM;
 196             memcpy(destPtr, srcPtr, MED_TAILLE_PNOM);
 197
 198             destPtr += MED_TAILLE_PNOM;
 199             itSet++;
 200         }
 201         subset->mNames[MED_TAILLE_PNOM * subset->mNum] = '\0';
 202     }
 203
 204     return subset;
 205 }
 206
 207 void        Elements::extractSubSetFromNodes(const std::set<med_int>& pSetOfNodes,
 208                                    std::set<med_int>& pSubSetOfElements) const
 209 {
 210     if (&pSetOfNodes == &pSubSetOfElements) throw IllegalStateException("pSetOfNodes and pSubSetOfElements must be different !", __FILE__, __LINE__);
 211
 212     int     numOfNodes = pSetOfNodes.size();
 213     bool    keepElt = false;
 214     for (int itElt = 0; itElt < mNum; ++itElt)
 215     {
 216         keepElt = false;
 217         for (std::set<med_int>::iterator itNode = pSetOfNodes.begin();
 218              itNode != pSetOfNodes.end() && keepElt == false; ++itNode)
 219         {
 220             for (int itCon = 0; itCon < mNumNodesByElt && keepElt == false; ++itCon)
 221             {
 222                 if ((*itNode) == mCon[itElt * mNumNodesByElt + itCon])
 223                 {
 224                     keepElt = true;
 225                 }
 226             }
 227         }
 228         if (keepElt)
 229         {
 230             pSubSetOfElements.insert(itElt + 1);
 231         }
 232     }
 233 }
 234
 235 const set<med_int>& Elements::getSetOfNodes()
 236 {
 237     // lazy get: test if mSetOfNodes has already been built
 238     if (mSetOfNodes.size() == 0)
 239     {
 240         buildSetOfNodes();
 241     }
 242
 243     return mSetOfNodes;
 244 }
 245
 246
 247 set<med_int> Elements::getSetOfFamilies() const
 248 {
 249     // set of families is empty at the beginning
 250     set<med_int> setOfFamilies;
 251
 252     // for each element, add its family to the set
 253     for (int itElt = 0 ; itElt < mNum ; itElt++)
 254     {
 255         setOfFamilies.insert(mFamIdent[itElt]);
 256     }
 257
 258     return setOfFamilies;
 259 }
 260
 261
 262 void Elements::remap()
 263 {
 264         int itNode, size;
 265     // build set of nodes if necessary
 266     if (mSetOfNodes.size() == 0)
 267     {
 268         buildSetOfNodes();
 269     }
 270
 271     // build the map for indices convertion
 272     map<med_int, med_int> mapOldIndexToNewIndex;
 273     med_int newIndex = 1; // MED index start at 1
 274     for (set<med_int>::iterator itSet = mSetOfNodes.begin(); itSet != mSetOfNodes.end() ; itSet++)
 275     {
 276         med_int oldIndex = (*itSet);
 277         mapOldIndexToNewIndex.insert(make_pair(oldIndex, newIndex));
 278         newIndex++;
 279     }
 280
 281     // for each node referenced by this set of elements
 282     for (itNode = 0, size = mNum * mNumNodesByElt ; itNode < size ; itNode++)
 283     {
 284         // convert old index to new index (remap)
 285         mCon[itNode] = mapOldIndexToNewIndex[mCon[itNode]];
 286     }
 287
 288     buildSetOfNodes();
 289 }
 290
 291 void Elements::remap(std::set<med_int>& pSetOfNodes)
 292 {
 293         int itNode, size;
 294
 295         // build the map for indices convertion
 296     map<med_int, med_int> mapOldIndexToNewIndex;
 297     med_int newIndex = 1; // MED index start at 1
 298         itNode = 1;
 299     for (std::set<med_int>::iterator it = pSetOfNodes.begin(); it != pSetOfNodes.end(); ++it)
 300     {
 301         med_int oldIndex = *it;
 302         mapOldIndexToNewIndex.insert(make_pair(oldIndex, itNode));
 303                 ++itNode;
 304     }
 305         // for each node referenced by this set of elements
 306     for (itNode = 0, size = mNum * mNumNodesByElt ; itNode < size ; itNode++)
 307     {
 308         // convert old index to new index (remap).
 309         mCon[itNode] = mapOldIndexToNewIndex[mCon[itNode]];;
 310     }
 311     buildSetOfNodes();
 312
 313 }
 314
 315 Elements* Elements::mergePartial(Elements* pElements, int pOffset)
 316 {
 317     Elements* elements = new Elements();
 318
 319     //---------------------------------------------------------------------
 320     // Copy parameters
 321     //---------------------------------------------------------------------
 322     elements->mNum           = mNum + pElements->mNum;
 323
 324     if (mEntity != pElements->mEntity) throw IllegalStateException("entity should be the same", __FILE__, __LINE__);
 325     elements->mEntity        = mEntity;
 326
 327     elements->mGeom          = mGeom;          // e.g. 310 for a TETRA10
 328     elements->mNumNodesByElt = mNumNodesByElt; // e.g. 10 for a TETRA10
 329     if (mDim != pElements->mDim) throw IllegalStateException("dimension should be the same", __FILE__, __LINE__);
 330     elements->mDim           = mDim;           // e.g. 3 for a TETRA10
 331
 332     elements->mId            = NULL;
 333     elements->mNames         = NULL;
 334
 335     //---------------------------------------------------------------------
 336     // Allocate arrays
 337     //---------------------------------------------------------------------
 338     elements->mFamIdent     = new med_int[elements->mNum];
 339     elements->mCon          = new med_int[mNumNodesByElt * elements->mNum];
 340
 341     //---------------------------------------------------------------------
 342     // Copy familys id and connectivity.
 343     //---------------------------------------------------------------------
 344     memcpy(elements->mFamIdent, mFamIdent, mNum * sizeof(med_int));
 345     memcpy(elements->mFamIdent + mNum, pElements->mFamIdent, pElements->mNum * sizeof(med_int));
 346
 347     memcpy(elements->mCon, mCon, mNum * mNumNodesByElt * sizeof(med_int));
 348
 349     med_int* dst = elements->mCon + mNum * mNumNodesByElt;
 350     for (int i = 0, n = pElements->mNum * mNumNodesByElt ; i < n ; i++)
 351     {
 352         dst[i] = pElements->mCon[i] + pOffset;
 353     }
 354
 355     //cout << "WARNING: do no build set of nodes" << endl;
 356     //elements->buildSetOfNodes();
 357
 358     return elements;
 359 }
 360
 361
 362 Elements* Elements::mergePartial(vector<Elements*> pElements, vector<int>& pOffsets)
 363 {
 364     if (pElements.size() == 0) throw IllegalArgumentException("pElements should contain at least 1 element", __FILE__, __LINE__);
 365
 366     Elements* elements = new Elements();
 367
 368     //---------------------------------------------------------------------
 369     // Copy parameters
 370     //---------------------------------------------------------------------
 371     elements->mNum = mNum;
 372     for (unsigned i = 0 ; i < pElements.size() ; i++)
 373     {
 374         elements->mNum += pElements[i]->mNum;
 375
 376         if (mEntity != pElements[i]->mEntity) throw IllegalStateException("entity should be the same", __FILE__, __LINE__);
 377         if (mDim != pElements[i]->mDim) throw IllegalStateException("dimension should be the same", __FILE__, __LINE__);
 378     }
 379
 380     elements->mEntity        = mEntity;
 381     elements->mGeom          = mGeom;          // e.g. 310 for a TETRA10
 382     elements->mNumNodesByElt = mNumNodesByElt; // e.g. 10 for a TETRA10
 383     elements->mDim           = mDim;           // e.g. 3 for a TETRA10
 384
 385     elements->mId            = NULL;
 386     elements->mNames         = NULL;
 387
 388     //---------------------------------------------------------------------
 389     // Allocate arrays
 390     //---------------------------------------------------------------------
 391     elements->mFamIdent     = new med_int[elements->mNum];
 392     elements->mCon          = new med_int[mNumNodesByElt * elements->mNum];
 393
 394     //---------------------------------------------------------------------
 395     // Copy familys id and connectivity.
 396     //---------------------------------------------------------------------
 397     memcpy(elements->mFamIdent, mFamIdent, mNum * sizeof(med_int));
 398     memcpy(elements->mCon, mCon, mNum * mNumNodesByElt * sizeof(med_int));
 399
 400     int offsetFamIdent = mNum;
 401     int offsetCon      = mNum;
 402     for (unsigned j = 0 ; j < pElements.size() ; j++)
 403     {
 404         memcpy(elements->mFamIdent + offsetFamIdent, pElements[j]->mFamIdent, pElements[j]->mNum * sizeof(med_int));  offsetFamIdent += pElements[j]->mNum;
 405
 406         med_int* dst = elements->mCon + offsetCon * mNumNodesByElt;
 407         for (int i = 0, n = pElements[j]->mNum * mNumNodesByElt ; i < n ; i++)
 408         {
 409             dst[i] = pElements[j]->mCon[i] + pOffsets[j];
 410         }
 411         offsetCon += pElements[j]->mNum;
 412     }
 413
 414     //cout << "WARNING: do no build set of nodes" << endl;
 415     //elements->buildSetOfNodes();
 416
 417     return elements;
 418 }
 419
 420
 421 void Elements::buildSetOfNodes()
 422 {
 423     if (mSetOfNodes.size() != 0)
 424     {
 425         mSetOfNodes.clear();
 426     }
 427
 428     // for each node referenced by this set of elements
 429     for (int itNode = 0, size = mNum * mNumNodesByElt ; itNode < size ; itNode++)
 430     {
 431         // add the node to the set
 432         mSetOfNodes.insert(mCon[itNode]);
 433     }
 434 }
 435
 436
 437 void Elements::readMED(
 438     med_idt               pMEDfile,
 439     char*                 pMeshName,
 440     med_int               pMeshDim,
 441     med_entite_maillage   pEntity,
 442     med_geometrie_element pGeom)
 443 {
 444     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 445     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 446     if ((pMeshDim < 1) || (pMeshDim > 3)) throw IllegalArgumentException("", __FILE__, __LINE__);
 447
 448     reset();
 449
 450     mEntity        = pEntity;
 451     mGeom          = pGeom;
 452     mNumNodesByElt = mGeom % 100;
 453     mDim           = mGeom / 100;
 454
 455     mNum = MEDnEntMaa(
 456         pMEDfile,
 457         pMeshName,
 458         MED_CONN,  // type of information requested = CONNECTIVITY
 459         pEntity,
 460         pGeom,
 461         MED_NOD);  // nodal connectivity
 462
 463     if (mNum < 0) throw IOException("i/o error while reading information about elements in MED file", __FILE__, __LINE__);
 464
 465     if (mNum == 0)
 466     {
 467         // empty mesh
 468         return;
 469     }
 470
 471     mCon      = new med_int[mNumNodesByElt * mNum];
 472     mNames    = new char[MED_TAILLE_PNOM * mNum + 1];
 473     mId       = new med_int[mNum];
 474     mFamIdent = new med_int[mNum];
 475     med_booleen isNames;
 476     med_booleen isIdentifiers;
 477
 478     med_err ret = MEDelementsLire(
 479         pMEDfile,
 480         pMeshName,
 481         pMeshDim,
 482         mCon,
 483         MED_FULL_INTERLACE,
 484         mNames,
 485         &isNames,
 486         mId,
 487         &isIdentifiers,
 488         mFamIdent,
 489         mNum,
 490         mEntity,
 491         mGeom,
 492         MED_NOD); // NODAL CONNECTIVITY
 493
 494     if (ret != 0) throw IOException("i/o error while reading elements in MED file", __FILE__, __LINE__);
 495
 496     if (!isNames)
 497     {
 498         delete[] mNames;
 499         mNames = NULL;
 500     }
 501
 502     if (!isIdentifiers)
 503     {
 504         delete[] mId;
 505         mId = NULL;
 506     }
 507 }
 508
 509
 510 void Elements::writeMED(med_idt pMEDfile, char* pMeshName, med_int pMeshDim)
 511 {
 512     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 513     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 514     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 515     if ((pMeshDim < 1) || (pMeshDim > 3)) throw IllegalArgumentException("", __FILE__, __LINE__);
 516
 517     // special case: if no elements => do nothing
 518     if (mNum == 0) return;
 519
 520     med_err ret = MEDelementsEcr(
 521         pMEDfile,
 522         pMeshName,
 523         pMeshDim,
 524         mCon,
 525         MED_FULL_INTERLACE,
 526         mNames,
 527         isNames()?MED_VRAI:MED_FAUX,
 528         mId,
 529         isIdentifiers()?MED_VRAI:MED_FAUX,
 530         mFamIdent,
 531         mNum,
 532         mEntity,
 533         mGeom,
 534         MED_NOD); // NODAL CONNECTIVITY
 535
 536     if (ret != 0) throw IOException("i/o error while writing elements in MED file", __FILE__, __LINE__);
 537 }
 538
 539
 540 ostream& operator<<(ostream& pOs, Elements& pE)
 541 {
 542     char strEntity[16];
 543     switch (pE.mEntity)
 544     {
 545         case MED_MAILLE: strcpy(strEntity, "MED_MAILLE"); break;
 546         case MED_FACE:   strcpy(strEntity, "MED_FACE"); break;
 547         case MED_ARETE:  strcpy(strEntity, "MED_ARETE"); break;
 548         case MED_NOEUD:  strcpy(strEntity, "MED_NOEUD"); break;
 549         default:         strcpy(strEntity, "UNKNOWN"); break;
 550     }
 551
 552     pOs << "Elements: " << endl;
 553     pOs << "    #number  =" << pE.mNum << endl;
 554     pOs << "    Entity   =" << strEntity << endl;
 555     pOs << "    Geom     =" << pE.mGeom << endl;
 556     pOs << "    Has names?" << (pE.isNames()?"yes":"no") << endl;
 557     pOs << "    Has id   ?" << (pE.isIdentifiers()?"yes":"no") << endl;
 558
 559     {
 560         set<med_int> setOfFam = pE.getSetOfFamilies();
 561         if (setOfFam.size() == 0)
 562         {
 563             pOs << "    Families: #fam=0" << endl;
 564         }
 565         else
 566         {
 567             set<med_int>::iterator itFam = setOfFam.end();
 568             itFam--;
 569             pOs << "    Families: #fam=" << setOfFam.size() << " id_min=" << (*(setOfFam.begin())) << " id_max=" << (*itFam) << endl;
 570
 571             if (pE.mFlagPrintAll)
 572             {
 573                 for (int i = 0 ; i < pE.mNum; i++)
 574                 {
 575                     pOs << "        Elt " << (i + 1) << ": " << pE.mFamIdent[i] << endl;
 576                 }
 577             }
 578         }
 579     }
 580
 581     {
 582         set<med_int> setOfNodes = pE.getSetOfNodes();
 583         if (setOfNodes.size() == 0)
 584         {
 585             pOs << "    Connectivity: #nodes=0" << endl;
 586         }
 587         else
 588         {
 589             set<med_int>::iterator itNode = setOfNodes.end();
 590             itNode--;
 591             pOs << "    Connectivity: #nodes=" << setOfNodes.size() << " id_min=" << (*(setOfNodes.begin())) << " id_max=" << (*itNode) << endl;
 592
 593             if (pE.mFlagPrintAll)
 594             {
 595                 for (int i = 0 ; i < pE.mNum ; i++)
 596                 {
 597                     pOs << "        Elt " << (i + 1) << ": ";
 598                     for (int j = 0 ; j < pE.mNumNodesByElt ; j++)
 599                     {
 600                         pOs << pE.mCon[i * pE.mNumNodesByElt + j] << " ";
 601                     }
 602                     pOs << endl;
 603                 }
 604             }
 605         }
 606     }
 607
 608     if (pE.mFlagPrintAll)
 609     {
 610
 611         if (pE.isIdentifiers())
 612         {
 613             pOs << "    Num (identifier): " << endl;
 614             for (int i = 0 ; i < pE.mNum; i++)
 615             {
 616                 pOs << "        Elt " << (i + 1) << ": " << pE.mId[i] << " " << endl;
 617             }
 618         }
 619
 620         if (pE.isNames())
 621         {
 622             pE.mNames[MED_TAILLE_PNOM * pE.mNum] = '\0';
 623             pOs << "    Names: |" << pE.mNames << "|" << endl;
 624         }
 625
 626     }
 627
 628     return pOs;
 629 }
 630
 631
 632 } // namespace  multipr
 633
 634 // EOF