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