src/MULTIPR/MULTIPR_Field.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_Field.cxx
   6  *
   7  * \brief   see MULTIPR_Field.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_Field.hxx"
  19 #include "MULTIPR_Exceptions.hxx"
  20
  21 #include <iostream>
  22
  23 using namespace std;
  24
  25
  26 namespace multipr
  27 {
  28
  29
  30 //*****************************************************************************
  31 // Class Field implementation
  32 //*****************************************************************************
  33
  34 Field::Field()
  35 {
  36     reset();
  37 }
  38
  39
  40 Field::~Field()
  41 {
  42     reset();
  43 }
  44
  45
  46 void Field::reset()
  47 {
  48     mName[0]       = '\0';
  49     mEntity        = MED_NOEUD;
  50     mGeom          = MED_NONE;
  51     mType          = MED_FLOAT64;
  52     mSizeOfType    = 8;
  53     mNumComponents = 0;
  54     mStrComponent  = "";
  55     mStrUnit       = "";
  56
  57     mNGauss.clear();
  58     mDT.clear();
  59     mNumDT.clear();
  60     mDTUnit.clear();
  61     mNumO.clear();
  62     mGaussLoc.clear();
  63     mProfil.clear();
  64     mSizeOfData.clear();
  65     mNVal.clear();
  66
  67     for (unsigned it = 0 ; it < mVal.size() ; it++)
  68     {
  69         delete[] mVal[it];
  70     }
  71     mVal.clear();
  72
  73     mFlagPrintAll = false;
  74 }
  75
  76
  77 bool Field::isEmpty() const
  78 {
  79     return (mNGauss.size() == 0);
  80 }
  81
  82
  83 int Field::getNumberOfGaussPointsByElement(int pTimeStepIt) const
  84 {
  85     if ((pTimeStepIt < 1) || (pTimeStepIt > int(mNGauss.size()))) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
  86
  87     return mNGauss[pTimeStepIt - 1];
  88 }
  89
  90
  91 const string& Field::getNameGaussLoc(int pTimeStepIt) const
  92 {
  93     if ((pTimeStepIt < 1) || (pTimeStepIt > int(mNGauss.size()))) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
  94
  95     return mGaussLoc[pTimeStepIt - 1];
  96 }
  97
  98
  99 const unsigned char* Field::getValue(int pTimeStepIt, int pIndex) const
 100 {
 101     if ((pTimeStepIt < 1) || (pTimeStepIt > int(mNGauss.size()))) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 102     if ((pIndex < 1) || (pIndex > mNVal[pTimeStepIt - 1] / mNGauss[pTimeStepIt - 1])) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 103
 104     // size of one data = size of type * number of components * number of Gauss points
 105     int sizeOfOneData = mSizeOfType * mNumComponents * mNGauss[pTimeStepIt - 1];
 106
 107     unsigned char* ret = mVal[pTimeStepIt - 1] + (pIndex - 1) * sizeOfOneData;
 108
 109     return ret;
 110 }
 111
 112
 113 Field* Field::extractSubSet(const set<med_int>& pSetIndices) const
 114 {
 115     Field* subset = new Field();
 116
 117     memcpy(subset->mName, mName, MED_TAILLE_NOM + 1);
 118     subset->mEntity        = mEntity;
 119     subset->mGeom          = mGeom;
 120     subset->mType          = mType;
 121     subset->mSizeOfType    = mSizeOfType;
 122     subset->mNumComponents = mNumComponents;
 123     subset->mStrComponent  = mStrComponent;
 124     subset->mStrUnit       = mStrUnit;
 125
 126     subset->mNGauss        = mNGauss;
 127     subset->mDT            = mDT;
 128     subset->mNumDT         = mNumDT;
 129     subset->mDTUnit        = mDTUnit;
 130     subset->mNumO          = mNumO;
 131     subset->mGaussLoc      = mGaussLoc;
 132     subset->mProfil        = mProfil;
 133
 134     // for each time step
 135     for (unsigned itTimeStep = 0 ; itTimeStep < mNGauss.size() ; itTimeStep++)
 136     {
 137         if (mProfil[itTimeStep].size() != 0) throw UnsupportedOperationException("", __FILE__, __LINE__);
 138         // WARNING : do not manage profil for the moment
 139         // if there is a profil,
 140         // 1. we should set mProfil to NO_PROFIL for this time_step
 141         // 2. we should extract data according to the profil
 142
 143         int nval = pSetIndices.size() * subset->mNGauss[itTimeStep];
 144         subset->mNVal.push_back(nval);
 145
 146         int sizeOfData = nval * mSizeOfType * mNumComponents;
 147         subset->mSizeOfData.push_back(sizeOfData);
 148
 149         unsigned char* data = new unsigned char[sizeOfData];
 150         unsigned char* dest = data;
 151         int sizeOfOneData = mSizeOfType * mNumComponents * subset->mNGauss[itTimeStep];
 152
 153         // for each element to extract
 154         for (set<med_int>::iterator itSet = pSetIndices.begin() ; itSet != pSetIndices.end() ; itSet++)
 155         {
 156             int indexElt = (*itSet);
 157
 158             // MED index start at 1.
 159             if (indexElt < 1) throw new IllegalArgumentException("", __FILE__, __LINE__);
 160
 161             unsigned char* src = mVal[itTimeStep] + (indexElt - 1) * sizeOfOneData;
 162             memcpy(dest, src, sizeOfOneData);
 163
 164             dest += sizeOfOneData;
 165         }
 166         subset->mVal.push_back(data);
 167     }
 168
 169     return subset;
 170 }
 171
 172
 173 void Field::getSetOfGaussLoc(set<string>& pSetOfGaussLoc) const
 174 {
 175     for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 176     {
 177         const string& gaussLocName = mGaussLoc[itGaussLoc];
 178
 179         if (gaussLocName.length() != 0)
 180         {
 181             pSetOfGaussLoc.insert(gaussLocName);
 182         }
 183     }
 184 }
 185
 186
 187 void Field::readMED(med_idt pMEDfile, med_int pIndex, char* pMeshName)
 188 {
 189     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 190     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 191     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 192     if (pIndex < 1) throw IllegalArgumentException("", __FILE__, __LINE__);
 193
 194     reset();
 195
 196     mNumComponents = MEDnChamp(pMEDfile, pIndex);
 197
 198     if (mNumComponents < 0) throw IOException("", __FILE__, __LINE__);
 199
 200     char* strComponent = new char[mNumComponents * MED_TAILLE_PNOM + 1];
 201     char* strUnit      = new char[mNumComponents * MED_TAILLE_PNOM + 1];
 202
 203     strComponent[0] = '\0';
 204     strUnit[0] = '\0';
 205
 206     med_err ret = MEDchampInfo(
 207         pMEDfile,
 208         pIndex,
 209         mName,
 210         &(mType),
 211         strComponent,
 212         strUnit,
 213         mNumComponents);
 214
 215     if (ret != 0) throw IOException("", __FILE__, __LINE__);
 216
 217     mStrComponent = strComponent;
 218     mStrUnit = strUnit;
 219
 220     delete[] strUnit;
 221     delete[] strComponent;
 222
 223     switch (mType)
 224     {
 225         case MED_FLOAT64: mSizeOfType = 8; break;
 226         case MED_INT64: mSizeOfType = 8; break;
 227         case MED_INT32: mSizeOfType = 4; break;
 228         case MED_INT: mSizeOfType = 4; break;
 229         default: throw IllegalStateException("should not be there", __FILE__, __LINE__);
 230     }
 231
 232     //---------------------------------------------------------------------
 233     // Read fields over nodes
 234     //---------------------------------------------------------------------
 235     bool fieldOnNodes = false;
 236     {
 237         med_int numTimeStepNodes = MEDnPasdetemps(
 238             pMEDfile,
 239             mName,
 240             MED_NOEUD,
 241             (med_geometrie_element) 0);
 242
 243         if (numTimeStepNodes < 0) throw IOException("", __FILE__, __LINE__);
 244
 245         if (numTimeStepNodes != 0)
 246         {
 247             fieldOnNodes = true;
 248             mEntity = MED_NOEUD;
 249             mGeom = (med_geometrie_element) 0;
 250             readMEDtimeSteps(pMEDfile, numTimeStepNodes, pMeshName);
 251         }
 252     }
 253
 254     //---------------------------------------------------------------------
 255     // Read fields over elements
 256     //---------------------------------------------------------------------
 257     {
 258         med_int numTimeStepElt = MEDnPasdetemps(
 259             pMEDfile,
 260             mName,
 261             MED_MAILLE,
 262             MED_TETRA10);
 263
 264         if (numTimeStepElt  < 0) throw IOException("", __FILE__, __LINE__);
 265
 266         if (numTimeStepElt != 0)
 267         {
 268             if (fieldOnNodes) throw IllegalStateException("", __FILE__, __LINE__);
 269
 270             mEntity = MED_MAILLE;
 271             mGeom = MED_TETRA10;
 272             readMEDtimeSteps(pMEDfile, numTimeStepElt, pMeshName);
 273         }
 274     }
 275 }
 276
 277
 278 void Field::readMEDtimeSteps(med_idt pMEDfile, med_int pNumberOfTimeSteps, char* pMeshName)
 279 {
 280     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 281     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 282     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 283     if (pNumberOfTimeSteps < 0) throw IllegalArgumentException("", __FILE__, __LINE__);
 284
 285     char strEntity[8];
 286     switch (mEntity)
 287     {
 288         case MED_ARETE:
 289         case MED_FACE:
 290         case MED_MAILLE: strcpy(strEntity, "CELLS"); break;
 291         case MED_NOEUD: strcpy(strEntity, "NODES"); break;
 292         default: strcpy(strEntity, "UNKNOWN"); break;
 293     }
 294
 295     // iterates over time step
 296     for (int itTimeStep = 1 ; itTimeStep <= pNumberOfTimeSteps ; itTimeStep++)
 297     {
 298         med_int ngauss;
 299         med_int numdt;
 300         med_int numo;
 301         char dtunit[MED_TAILLE_PNOM + 1];
 302         med_float dt;
 303         char maa[MED_TAILLE_NOM + 1];
 304         med_booleen local;
 305         med_int nmaa;
 306
 307         med_err ret = MEDpasdetempsInfo(
 308             pMEDfile,
 309             mName,
 310             mEntity,
 311             mGeom,
 312             itTimeStep,
 313                 &ngauss,
 314             &numdt,
 315             &numo,
 316             dtunit,
 317             &dt,
 318             maa,
 319             &local,
 320             &nmaa);
 321
 322         if (ret != 0) throw IOException("i/o error while reading #timesteps in MED file", __FILE__, __LINE__);
 323
 324         // mesh must be local
 325         if (local != MED_VRAI) throw IllegalStateException("only local fields are currently supported", __FILE__, __LINE__);
 326
 327         // #mesh must be 1
 328         if (nmaa != 1) throw IllegalStateException("field shoud be associated with 1 mesh only", __FILE__, __LINE__);
 329
 330         // mesh must be pMeshName
 331         // if field does not apply on the current mesh, skip
 332         if (strcmp(maa, pMeshName) != 0)
 333         {
 334             continue;
 335         }
 336
 337         mNGauss.push_back(ngauss);
 338         mDT.push_back(dt);
 339         mNumDT.push_back(numdt);
 340         mDTUnit.push_back(dtunit);
 341         mNumO.push_back(numo);
 342
 343         med_int nval = MEDnVal(
 344             pMEDfile,
 345             mName,
 346             mEntity,
 347             mGeom,
 348             numdt,
 349             numo,
 350             pMeshName,
 351             MED_GLOBAL);
 352
 353         if (nval < 0) throw IOException("i/o error while reading field in MED file", __FILE__, __LINE__);
 354
 355         mNVal.push_back(nval);
 356
 357         char gaussLocName[MED_TAILLE_NOM + 1];
 358         char profilName[MED_TAILLE_NOM + 1];
 359         int sizeOfData = mSizeOfType * mNumComponents * nval;
 360         mSizeOfData.push_back(sizeOfData);
 361         unsigned char* fieldData = new unsigned char[sizeOfData];
 362
 363         ret = MEDchampLire(
 364             pMEDfile,
 365             pMeshName,
 366             mName,
 367             fieldData,
 368             MED_FULL_INTERLACE,
 369             MED_ALL,
 370             gaussLocName,
 371             profilName,
 372             MED_COMPACT, // should be: MED_GLOBAL, MED_NO_PFLMOD, MED_COMPACT
 373             mEntity,
 374             mGeom,
 375             numdt,
 376             numo);
 377
 378         if (ret != 0) throw IOException("i/o error while reading field in MED file", __FILE__, __LINE__);
 379
 380         mGaussLoc.push_back(gaussLocName);
 381         mProfil.push_back(profilName);
 382         mVal.push_back(fieldData);
 383     }
 384 }
 385
 386
 387 void Field::writeMED(med_idt pMEDfile, char* pMeshName)
 388 {
 389     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 390     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 391     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 392     if (mNumComponents < 1) throw IllegalStateException("", __FILE__, __LINE__);
 393
 394     med_err ret = MEDchampCr(
 395         pMEDfile,
 396         mName,                                        // name of the field
 397         mType,                                        // type of data (MED_FLOAT64, MED_INT32, etc.)
 398         const_cast<char*>(mStrComponent.c_str()),     // name of components
 399         const_cast<char*>(mStrUnit.c_str()),          // name of units
 400         mNumComponents);                              // number of components
 401
 402     if (ret != 0) throw IOException("i/o error while creating field in MED file", __FILE__, __LINE__);
 403
 404     // for each time step
 405     for (unsigned i = 0 ; i < mNGauss.size() ; i++)
 406     {
 407         // skip if no values
 408         if (mNVal[i] == 0) continue;
 409
 410         ret = MEDchampEcr(
 411             pMEDfile,
 412             pMeshName,          // name of the mesh (first call to MEDchampEcr => name of reference)
 413             mName,              // name of the field
 414             mVal[i],            // data (= values)
 415             MED_FULL_INTERLACE, // data organization
 416             mNVal[i],           // number of values
 417             const_cast<char*>(mGaussLoc[i].c_str()), // name of Gauss reference
 418             MED_ALL,            // components to be selected
 419             const_cast<char*>(mProfil[i].c_str()), // name of profil
 420             MED_GLOBAL,         // how to read data: MED_NO_PFLMOD,MED_COMPACT,MED_GLOBAL
 421             mEntity,            // type of entity (MED_NOEUD, MED_MAILLE, etc.)
 422             mGeom,              // type of geometry (TETRA10, etc.)
 423             mNumDT[i],          // time step iteration
 424             const_cast<char*>(mDTUnit[i].c_str()), // unit of time step
 425             mDT[i],             // time key
 426             mNumO[i]);          // order number
 427
 428         if (ret != 0) throw IOException("i/o error while writing field in MED file", __FILE__, __LINE__);
 429     }
 430 }
 431
 432
 433 ostream& operator<<(ostream& pOs, Field& pF)
 434 {
 435     char strEntity[16];
 436     switch (pF.mEntity)
 437     {
 438         case MED_MAILLE: strcpy(strEntity, "MED_MAILLE"); break;
 439         case MED_FACE:   strcpy(strEntity, "MED_FACE"); break;
 440         case MED_ARETE:  strcpy(strEntity, "MED_ARETE"); break;
 441         case MED_NOEUD:  strcpy(strEntity, "MED_NOEUD"); break;
 442         default:         strcpy(strEntity, "UNKNOWN"); break;
 443     }
 444
 445     char strType[16];
 446     switch (pF.mType)
 447     {
 448         case MED_FLOAT64: strcpy(strType, "MED_FLOAT64"); break;
 449         case MED_INT32:   strcpy(strType, "MED_INT32"); break;
 450         case MED_INT64:   strcpy(strType, "MED_INT64"); break;
 451         case MED_INT:     strcpy(strType, "MED_INT"); break;
 452         default:          strcpy(strType, "UNKNOWN"); break;
 453     }
 454
 455     pOs << "Field: " << endl;
 456     pOs << "    Name       =|" << pF.mName << "|" << endl;
 457     pOs << "    Entity     =" << strEntity << endl;
 458     pOs << "    Geom       =" << pF.mGeom << endl;
 459     pOs << "    Type       =" << strType << " (size=" << pF.mSizeOfType << ")" << endl;
 460     pOs << "    #Components=" << pF.mNumComponents << endl;
 461     pOs << "        Name component=|" << pF.mStrComponent << "|" << endl;
 462     pOs << "        Unit component=|" << pF.mStrUnit << "|" << endl;
 463     pOs << "    #Time steps=" << pF.mNGauss.size() << endl;
 464
 465     for (unsigned itTimeStep = 0 ; itTimeStep < pF.mNGauss.size() ; itTimeStep++)
 466     {
 467         pOs << "        Time=" << pF.mDT[itTimeStep];
 468         pOs << "  it=" << pF.mNumDT[itTimeStep];
 469         pOs << "  order=" << pF.mNumO[itTimeStep];
 470         pOs << "  #gauss=" << pF.mNGauss[itTimeStep];
 471         pOs << "  #val=" << pF.mNVal[itTimeStep];
 472         pOs << "  sizeof_val=" << pF.mSizeOfData[itTimeStep];
 473         pOs << "  gauss_loc=|" << ((pF.mGaussLoc[itTimeStep].size() == 0)?"NONE":pF.mGaussLoc[itTimeStep]) << "|  size=" << pF.mGaussLoc[itTimeStep].size();
 474         pOs << "  profil=|" << ((pF.mProfil[itTimeStep].size() == 0)?"NONE":pF.mProfil[itTimeStep]) << "|  size=" << pF.mProfil[itTimeStep].size() << endl;
 475
 476         if (pF.mFlagPrintAll)
 477         {
 478             cout << "    Values: ";
 479             switch (pF.mType)
 480             {
 481             case MED_FLOAT64:
 482                 {
 483                     med_float* src = reinterpret_cast<med_float*>(pF.mVal[itTimeStep]);
 484                     for (int itVal = 0 ; itVal < pF.mNVal[itTimeStep] * pF.mNumComponents ; itVal++)
 485                     {
 486                         cout << src[itVal] << " ";
 487                     }
 488                 }
 489                 break;
 490             case MED_INT:
 491             case MED_INT32:
 492                 {
 493                     med_int* src = reinterpret_cast<med_int*>(pF.mVal[itTimeStep]);
 494                     for (int itVal = 0 ; itVal < pF.mNVal[itTimeStep] * pF.mNumComponents ; itVal++)
 495                     {
 496                         cout << src[itVal] << " ";
 497                     }
 498                 }
 499                 break;
 500             case MED_INT64:
 501                 // not yet implemented
 502                 throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 503             default:
 504                 // should not be there
 505                 throw IllegalStateException("should not be there", __FILE__, __LINE__);
 506             }
 507             cout << endl;
 508         }
 509
 510     }
 511
 512     return pOs;
 513 }
 514
 515
 516 } // namespace  multipr
 517
 518 // EOF