src/MULTIPR/MULTIPR_Field.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Partitioning/decimation module for the SALOME v3.2 platform
  20 //
  21 /**
  22  * \file    MULTIPR_Field.cxx
  23  *
  24  * \brief   see MULTIPR_Field.hxx
  25  *
  26  * \author  Olivier LE ROUX - CS, Virtual Reality Dpt
  27  *
  28  * \date    01/2007
  29  */
  30
  31 //*****************************************************************************
  32 // Includes section
  33 //*****************************************************************************
  34
  35 #include "MULTIPR_Field.hxx"
  36 #include "MULTIPR_MeshDis.hxx"
  37 #include "MULTIPR_Mesh.hxx"
  38 #include "MULTIPR_Family.hxx"
  39 #include "MULTIPR_Exceptions.hxx"
  40 #include "MULTIPR_Elements.hxx"
  41 #include "MULTIPR_Profil.hxx"
  42
  43 #include <iostream>
  44 #include <limits>
  45
  46 using namespace std;
  47
  48
  49 namespace multipr
  50 {
  51
  52 //*****************************************************************************
  53 // Class Field implementation
  54 //*****************************************************************************
  55
  56 Field::Field()
  57 {
  58     reset();
  59 }
  60
  61
  62 Field::~Field()
  63 {
  64     reset();
  65 }
  66
  67
  68 void Field::reset()
  69 {
  70     mName[0]       = '\0';
  71     mEntity        = MED_NOEUD;
  72     mGeom          = MED_NONE;
  73         mGeomIdx           = eMaxMedMesh;
  74     mType          = MED_FLOAT64;
  75     mSizeOfType    = 8;
  76     mNumComponents = 0;
  77     mStrComponent  = "";
  78     mStrUnit       = "";
  79
  80     mNGauss.clear();
  81     mDT.clear();
  82     mNumDT.clear();
  83     mDTUnit.clear();
  84     mNumO.clear();
  85     mGaussLoc.clear();
  86     mProfil.clear();
  87     mSizeOfData.clear();
  88     mNVal.clear();
  89
  90     for (unsigned it = 0 ; it < mVal.size() ; it++)
  91     {
  92         delete[] mVal[it];
  93     }
  94     mVal.clear();
  95
  96     mFlagPrintAll = false;
  97 }
  98
  99
 100 bool Field::isEmpty() const
 101 {
 102     return (mNGauss.size() == 0);
 103 }
 104
 105
 106 int Field::getNumberOfGaussPointsByElement(int pTimeStepIt) const
 107 {
 108     if ((pTimeStepIt < 1) || (pTimeStepIt > int(mNGauss.size()))) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 109
 110     return mNGauss[pTimeStepIt - 1];
 111 }
 112
 113
 114 const string& Field::getNameGaussLoc(int pTimeStepIt) const
 115 {
 116     if ((pTimeStepIt < 1) || (pTimeStepIt > int(mNGauss.size()))) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 117
 118     return mGaussLoc[pTimeStepIt - 1];
 119 }
 120
 121
 122 const unsigned char* Field::getValue(int pTimeStepIt, int pIndex) const
 123 {
 124     if ((pTimeStepIt < 1) || (pTimeStepIt > int(mNGauss.size()))) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 125     if ((pIndex < 1) || (pIndex > mNVal[pTimeStepIt - 1] / mNGauss[pTimeStepIt - 1])) throw IndexOutOfBoundsException("", __FILE__, __LINE__);
 126
 127     // size of one data = size of type * number of components * number of Gauss points
 128     int sizeOfOneData = mSizeOfType * mNumComponents * mNGauss[pTimeStepIt - 1];
 129
 130     unsigned char* ret = mVal[pTimeStepIt - 1] + (pIndex - 1) * sizeOfOneData;
 131
 132     return ret;
 133 }
 134
 135 const std::string& Field::getProfil(int pTimeStep) const
 136 {
 137     if ((pTimeStep < 1) || (pTimeStep > int(mProfil.size()))) throw IndexOutOfBoundsException("Unknown time step", __FILE__, __LINE__);
 138
 139     return mProfil[pTimeStep - 1];
 140 }
 141
 142 Field* Field::extractSubSet(const set<med_int>& pSetIndices) const
 143 {
 144     Field* subset = new Field();
 145
 146     memcpy(subset->mName, mName, MED_TAILLE_NOM + 1);
 147     subset->mEntity        = mEntity;
 148     subset->mGeom          = mGeom;
 149     subset->mType          = mType;
 150     subset->mSizeOfType    = mSizeOfType;
 151     subset->mNumComponents = mNumComponents;
 152     subset->mStrComponent  = mStrComponent;
 153     subset->mStrUnit       = mStrUnit;
 154
 155     subset->mNGauss        = mNGauss;
 156     subset->mDT            = mDT;
 157     subset->mNumDT         = mNumDT;
 158     subset->mDTUnit        = mDTUnit;
 159     subset->mNumO          = mNumO;
 160     subset->mGaussLoc      = mGaussLoc;
 161     subset->mProfil        = mProfil;
 162
 163     // for each time step
 164     for (unsigned itTimeStep = 0 ; itTimeStep < mNGauss.size() ; itTimeStep++)
 165     {
 166         //if (mProfil[itTimeStep].size() != 0) throw UnsupportedOperationException("", __FILE__, __LINE__);
 167         // WARNING : do not manage profil for the moment
 168         // if there is a profil,
 169         // 1. we should set mProfil to NO_PROFIL for this time_step
 170         // 2. we should extract data according to the profil
 171
 172         int nval = pSetIndices.size() * subset->mNGauss[itTimeStep];
 173         subset->mNVal.push_back(nval);
 174
 175         int sizeOfData = nval * mSizeOfType * mNumComponents;
 176         subset->mSizeOfData.push_back(sizeOfData);
 177
 178         unsigned char* data = new unsigned char[sizeOfData];
 179         unsigned char* dest = data;
 180         int sizeOfOneData = mSizeOfType * mNumComponents * subset->mNGauss[itTimeStep];
 181
 182         // for each element to extract
 183         for (set<med_int>::iterator itSet = pSetIndices.begin() ; itSet != pSetIndices.end() ; itSet++)
 184         {
 185             int indexElt = (*itSet);
 186
 187             // MED index start at 1.
 188             if (indexElt < 1) throw new IllegalArgumentException("", __FILE__, __LINE__);
 189
 190             unsigned char* src = mVal[itTimeStep] + (indexElt - 1) * sizeOfOneData;
 191             memcpy(dest, src, sizeOfOneData);
 192
 193             dest += sizeOfOneData;
 194         }
 195         subset->mVal.push_back(data);
 196     }
 197
 198     return subset;
 199 }
 200
 201
 202 Field* Field::merge(Field* pField)
 203 {
 204     Field* field = new Field();
 205
 206     if (strcmp(mName, pField->mName) != 0) throw IllegalStateException("incompatible field", __FILE__, __LINE__);
 207     memcpy(field->mName, mName, MED_TAILLE_NOM + 1);
 208
 209     if (mEntity != pField->mEntity) throw IllegalStateException("incompatible entity", __FILE__, __LINE__);
 210     field->mEntity        = mEntity;
 211
 212     if (mGeom != pField->mGeom) throw IllegalStateException("incompatible geom", __FILE__, __LINE__);
 213     field->mGeom          = mGeom;
 214
 215     if (mType != pField->mType) throw IllegalStateException("incompatible type", __FILE__, __LINE__);
 216     field->mType          = mType;
 217
 218     field->mSizeOfType    = mSizeOfType;
 219
 220     if (mNumComponents != pField->mNumComponents) throw IllegalStateException("incompatible #components", __FILE__, __LINE__);
 221     field->mNumComponents = mNumComponents;
 222
 223     field->mStrComponent  = mStrComponent;
 224     field->mStrUnit       = mStrUnit;
 225
 226     if (mNGauss.size() != pField->mNGauss.size()) throw IllegalStateException("incompatible #time stamps", __FILE__, __LINE__);
 227     field->mNGauss        = mNGauss;
 228
 229     field->mDT            = mDT;
 230     field->mNumDT         = mNumDT;
 231     field->mDTUnit        = mDTUnit;
 232     field->mNumO          = mNumO;
 233     field->mGaussLoc      = mGaussLoc;
 234     field->mProfil        = mProfil;
 235
 236     // for each time step
 237     for (unsigned itTimeStep = 0 ; itTimeStep < mNGauss.size() ; itTimeStep++)
 238     {
 239         int sizeOfData = mSizeOfData[itTimeStep] + pField->mSizeOfData[itTimeStep];
 240         field->mSizeOfData.push_back(sizeOfData);
 241
 242         int nVal = mNVal[itTimeStep] + pField->mNVal[itTimeStep];
 243         field->mNVal.push_back(nVal);
 244
 245         unsigned char* data = new unsigned char[sizeOfData];
 246         memcpy(data, mVal[itTimeStep], mSizeOfData[itTimeStep]);
 247         memcpy(data + mSizeOfData[itTimeStep], pField->mVal[itTimeStep], pField->mSizeOfData[itTimeStep]);
 248         field->mVal.push_back(data);
 249     }
 250
 251     return field;
 252 }
 253
 254
 255 Field* Field::merge(vector<Field*> pFields, int pFieldIt)
 256 {
 257     if (pFields.size() == 0) throw IllegalArgumentException("pElements should contain at least 1 element", __FILE__, __LINE__);
 258
 259     Field* field = new Field();
 260
 261     // check if fields are compatible
 262     for (unsigned i = 0 ; i < pFields.size() ; i++)
 263     {
 264         if (strcmp(mName, pFields[i]->mName) != 0) throw IllegalStateException("incompatible field", __FILE__, __LINE__);
 265         if (mEntity != pFields[i]->mEntity) throw IllegalStateException("incompatible entity", __FILE__, __LINE__);
 266         if (mGeom != pFields[i]->mGeom) throw IllegalStateException("incompatible geom", __FILE__, __LINE__);
 267         if (mType != pFields[i]->mType) throw IllegalStateException("incompatible type", __FILE__, __LINE__);
 268         if (mNumComponents != pFields[i]->mNumComponents) throw IllegalStateException("incompatible #components", __FILE__, __LINE__);
 269         if (mNGauss.size() != pFields[i]->mNGauss.size()) throw IllegalStateException("incompatible #time stamps", __FILE__, __LINE__);
 270     }
 271
 272     memcpy(field->mName, mName, MED_TAILLE_NOM + 1);
 273     field->mEntity        = mEntity;
 274     field->mGeom          = mGeom;
 275     field->mType          = mType;
 276     field->mSizeOfType    = mSizeOfType;
 277     field->mNumComponents = mNumComponents;
 278     field->mStrComponent  = mStrComponent;
 279     field->mStrUnit       = mStrUnit;
 280
 281     if (pFieldIt == -1) // merge all time step
 282     {
 283         field->mNGauss        = mNGauss;
 284         field->mDT            = mDT;
 285         field->mNumDT         = mNumDT;
 286         field->mDTUnit        = mDTUnit;
 287         field->mNumO          = mNumO;
 288         field->mGaussLoc      = mGaussLoc;
 289         field->mProfil        = mProfil;
 290
 291         // for each time step
 292         for (unsigned itTimeStep = 0 ; itTimeStep < mNGauss.size() ; itTimeStep++)
 293         {
 294             int sizeOfData = mSizeOfData[itTimeStep];
 295             int nVal = mNVal[itTimeStep];
 296
 297             for (unsigned i = 0 ; i < pFields.size() ; i++)
 298             {
 299                 sizeOfData += pFields[i]->mSizeOfData[itTimeStep];
 300                 nVal += pFields[i]->mNVal[itTimeStep];
 301             }
 302
 303             field->mSizeOfData.push_back(sizeOfData);
 304             field->mNVal.push_back(nVal);
 305
 306             unsigned char* data = new unsigned char[sizeOfData];
 307             memcpy(data, mVal[itTimeStep], mSizeOfData[itTimeStep]);
 308
 309             int offsetData = mSizeOfData[itTimeStep];
 310             for (unsigned i = 0 ; i < pFields.size() ; i++)
 311             {
 312                 memcpy(data + offsetData, pFields[i]->mVal[itTimeStep], pFields[i]->mSizeOfData[itTimeStep]);
 313                 offsetData += pFields[i]->mSizeOfData[itTimeStep];
 314             }
 315
 316             field->mVal.push_back(data);
 317         }
 318     }
 319     else // only merge the given time step
 320     {
 321         field->mNGauss.push_back(mNGauss[pFieldIt]);
 322         field->mDT.push_back(mDT[pFieldIt]);
 323         field->mNumDT.push_back(mNumDT[pFieldIt]);
 324         field->mDTUnit.push_back(mDTUnit[pFieldIt]);
 325         field->mNumO.push_back(mNumO[pFieldIt]);
 326
 327         if (mGaussLoc.size() != 0) field->mGaussLoc.push_back(mGaussLoc[pFieldIt]);
 328         if (mProfil.size() != 0) field->mProfil.push_back(mProfil[pFieldIt]);
 329
 330         // to finish
 331         throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 332     }
 333
 334     return field;
 335 }
 336
 337 void Field::getMinMax(float& pMin, float& pMax)
 338 {
 339     pMin = std::numeric_limits<med_float>::max();
 340     pMax = std::numeric_limits<med_float>::min();
 341
 342     for (unsigned itTimeStep = 0 ; itTimeStep < mNGauss.size() ; itTimeStep++)
 343     {
 344         switch (mType)
 345         {
 346         case MED_FLOAT64:
 347             {
 348                 med_float* src = reinterpret_cast<med_float*>(mVal[itTimeStep]);
 349                 for (int itVal = 0 ; itVal < mNVal[itTimeStep] * mNumComponents ; itVal++)
 350                 {
 351                     if (src[itVal] < pMin)
 352                     {
 353                         pMin = src[itVal];
 354                     }
 355                     if (src[itVal] > pMax)
 356                     {
 357                         pMax = src[itVal];
 358                     }
 359                 }
 360             }
 361             break;
 362         case MED_INT:
 363         case MED_INT32:
 364             {
 365                 med_int* src = reinterpret_cast<med_int*>(mVal[itTimeStep]);
 366                 for (int itVal = 0 ; itVal < mNVal[itTimeStep] * mNumComponents ; itVal++)
 367                 {
 368                     if (src[itVal] < pMin)
 369                     {
 370                         pMin = src[itVal];
 371                     }
 372                     if (src[itVal] > pMax)
 373                     {
 374                         pMax = src[itVal];
 375                     }
 376                 }
 377             }
 378             break;
 379         case MED_INT64:
 380             {
 381                 long* src = reinterpret_cast<long*>(mVal[itTimeStep]);
 382                 for (int itVal = 0 ; itVal < mNVal[itTimeStep] * mNumComponents ; itVal++)
 383                 {
 384                     if (src[itVal] < pMin)
 385                     {
 386                         pMin = src[itVal];
 387                     }
 388                     if (src[itVal] > pMax)
 389                     {
 390                         pMax = src[itVal];
 391                     }
 392                 }
 393             }
 394             break;
 395         default:
 396             // should not be there
 397             throw IllegalStateException("should not be there", __FILE__, __LINE__);
 398         }
 399     }
 400 }
 401
 402 void Field::getSetOfGaussLoc(set<string>& pSetOfGaussLoc) const
 403 {
 404     for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
 405     {
 406         const string& gaussLocName = mGaussLoc[itGaussLoc];
 407
 408         if (gaussLocName.length() != 0)
 409         {
 410             pSetOfGaussLoc.insert(gaussLocName);
 411         }
 412     }
 413 }
 414
 415
 416 void Field::readMED(med_idt pMEDfile, med_int pIndex, char* pMeshName, med_geometrie_element pGeom)
 417 {
 418     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 419     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 420     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 421     if (pIndex < 1) throw IllegalArgumentException("", __FILE__, __LINE__);
 422
 423     reset();
 424
 425     mNumComponents = MEDnChamp(pMEDfile, pIndex);
 426
 427     if (mNumComponents < 0) throw IOException("", __FILE__, __LINE__);
 428
 429     char* strComponent = new char[mNumComponents * MED_TAILLE_PNOM + 1];
 430     char* strUnit      = new char[mNumComponents * MED_TAILLE_PNOM + 1];
 431
 432     strComponent[0] = '\0';
 433     strUnit[0] = '\0';
 434
 435     med_err ret = MEDchampInfo(
 436         pMEDfile,
 437         pIndex,
 438         mName,
 439         &(mType),
 440         strComponent,
 441         strUnit,
 442         mNumComponents);
 443     mName[MED_TAILLE_NOM] = '\0';
 444     if (ret != 0) throw IOException("", __FILE__, __LINE__);
 445
 446     mStrComponent = strComponent;
 447     mStrUnit = strUnit;
 448
 449     delete[] strUnit;
 450     delete[] strComponent;
 451
 452     switch (mType)
 453     {
 454         case MED_FLOAT64: mSizeOfType = 8; break;
 455         case MED_INT64: mSizeOfType = 8; break;
 456         case MED_INT32: mSizeOfType = 4; break;
 457         case MED_INT: mSizeOfType = 4; break;
 458         default: throw IllegalStateException("Unknown field data type", __FILE__, __LINE__);
 459     }
 460
 461     //---------------------------------------------------------------------
 462     // Read fields over nodes
 463     //---------------------------------------------------------------------
 464     bool fieldOnNodes = false;
 465     {
 466         med_int numTimeStepNodes = MEDnPasdetemps(
 467             pMEDfile,
 468             mName,
 469             MED_NOEUD,
 470             (med_geometrie_element) 0);
 471
 472         if (numTimeStepNodes < 0) throw IOException("", __FILE__, __LINE__);
 473
 474         if (numTimeStepNodes > 0)
 475         {
 476             fieldOnNodes = true;
 477             mEntity = MED_NOEUD;
 478             mGeom = (med_geometrie_element) 0;
 479                         mGeomIdx = eMaxMedMesh;
 480             readMEDtimeSteps(pMEDfile, numTimeStepNodes, pMeshName);
 481         }
 482     }
 483
 484     //---------------------------------------------------------------------
 485     // Read fields over elements
 486     //---------------------------------------------------------------------
 487     {
 488         med_int numTimeStepElt = MEDnPasdetemps(
 489             pMEDfile,
 490             mName,
 491             MED_MAILLE,
 492             pGeom);
 493
 494                 for (int i = 0; i < eMaxMedMesh; ++i)
 495                 {
 496                         if (CELL_TYPES[i] == pGeom)
 497                         {
 498                                 mGeomIdx = (eMeshType)i;
 499                                 break;
 500                         }
 501                 }
 502         if (numTimeStepElt  < 0) throw IOException("", __FILE__, __LINE__);
 503
 504         if (numTimeStepElt > 0)
 505         {
 506             if (fieldOnNodes) throw IllegalStateException("", __FILE__, __LINE__);
 507             mEntity = MED_MAILLE;
 508             mGeom = pGeom;
 509             readMEDtimeSteps(pMEDfile, numTimeStepElt, pMeshName);
 510         }
 511     }
 512 }
 513
 514
 515 void Field::readMEDtimeSteps(med_idt pMEDfile, med_int pNumberOfTimeSteps, char* pMeshName)
 516 {
 517     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 518     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 519     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 520     if (pNumberOfTimeSteps < 0) throw IllegalArgumentException("", __FILE__, __LINE__);
 521
 522     char strEntity[8];
 523     switch (mEntity)
 524     {
 525         case MED_ARETE:
 526         case MED_FACE:
 527         case MED_MAILLE: strcpy(strEntity, "CELLS"); break;
 528         case MED_NOEUD: strcpy(strEntity, "NODES"); break;
 529         default: strcpy(strEntity, "UNKNOWN"); break;
 530     }
 531
 532     // iterates over time step
 533     for (int itTimeStep = 1 ; itTimeStep <= pNumberOfTimeSteps ; itTimeStep++)
 534     {
 535         med_int ngauss;
 536         med_int numdt;
 537         med_int numo;
 538         char dtunit[MED_TAILLE_PNOM + 1];
 539         med_float dt;
 540         char maa[MED_TAILLE_NOM + 1];
 541         med_booleen local;
 542         med_int nmaa;
 543
 544         med_err ret = MEDpasdetempsInfo(
 545             pMEDfile,
 546             mName,
 547             mEntity,
 548             mGeom,
 549             itTimeStep,
 550                         &ngauss,
 551             &numdt,
 552             &numo,
 553             dtunit,
 554             &dt,
 555             maa,
 556             &local,
 557             &nmaa);
 558
 559         if (ret != 0) throw IOException("i/o error while reading #timesteps in MED file", __FILE__, __LINE__);
 560
 561         // mesh must be local
 562         if (local != MED_VRAI) throw IllegalStateException("only local fields are currently supported", __FILE__, __LINE__);
 563
 564         // #mesh must be 1
 565         if (nmaa != 1) throw IllegalStateException("field shoud be associated with 1 mesh only", __FILE__, __LINE__);
 566
 567         // mesh must be pMeshName
 568         // if field does not apply on the current mesh, skip
 569         if (strcmp(maa, pMeshName) != 0)
 570         {
 571             continue;
 572         }
 573
 574         mNGauss.push_back(ngauss);
 575         mDT.push_back(dt);
 576         mNumDT.push_back(numdt);
 577         mDTUnit.push_back(dtunit);
 578         mNumO.push_back(numo);
 579
 580         med_int nval = MEDnVal(
 581             pMEDfile,
 582             mName,
 583             mEntity,
 584             mGeom,
 585             numdt,
 586             numo,
 587             pMeshName,
 588             MED_GLOBAL);
 589
 590         if (nval < 0) throw IOException("i/o error while reading field in MED file", __FILE__, __LINE__);
 591
 592         mNVal.push_back(nval);
 593
 594         char gaussLocName[MED_TAILLE_NOM + 1];
 595         char profilName[MED_TAILLE_NOM + 1];
 596         int sizeOfData = mSizeOfType * mNumComponents * nval;
 597         mSizeOfData.push_back(sizeOfData);
 598         unsigned char* fieldData = new unsigned char[sizeOfData];
 599         ret = MEDchampLire(
 600             pMEDfile,
 601             pMeshName,
 602             mName,
 603             fieldData,
 604             MED_FULL_INTERLACE,
 605             MED_ALL,
 606             gaussLocName,
 607             profilName,
 608             MED_COMPACT, // should be: MED_GLOBAL, MED_NO_PFLMOD, MED_COMPACT
 609             mEntity,
 610             mGeom,
 611             numdt,
 612             numo);
 613         if (ret != 0) throw IOException("i/o error while reading field in MED file", __FILE__, __LINE__);
 614
 615         mGaussLoc.push_back(gaussLocName);
 616         mProfil.push_back(profilName);
 617         mVal.push_back(fieldData);
 618     }
 619 }
 620
 621
 622 void Field::writeMED(med_idt pMEDfile, char* pMeshName, bool pCreateField)
 623 {
 624     if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 625     if (pMeshName == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 626     if (strlen(pMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("", __FILE__, __LINE__);
 627     if (mNumComponents < 1) throw IllegalStateException("", __FILE__, __LINE__);
 628     med_err ret;
 629
 630     if (pCreateField)
 631     {
 632         ret = MEDchampCr(
 633         pMEDfile,
 634         mName,                                        // name of the field
 635         mType,                                        // type of data (MED_FLOAT64, MED_INT32, etc.)
 636         const_cast<char*>(mStrComponent.c_str()),     // name of components
 637         const_cast<char*>(mStrUnit.c_str()),          // name of units
 638         mNumComponents);                              // number of components
 639
 640         if (ret != 0) throw IOException("i/o error while creating field in MED file", __FILE__, __LINE__);
 641     }
 642     // for each time step
 643     for (unsigned i = 0 ; i < mNGauss.size() ; i++)
 644     {
 645         // skip if no values
 646         if (mNVal[i] == 0) continue;
 647
 648         ret = MEDchampEcr(
 649             pMEDfile,
 650             pMeshName,          // name of the mesh (first call to MEDchampEcr => name of reference)
 651             mName,              // name of the field
 652             mVal[i],            // data (= values)
 653             MED_FULL_INTERLACE, // data organization
 654             mNVal[i],           // number of values
 655             const_cast<char*>(mGaussLoc[i].c_str()), // name of Gauss reference
 656             MED_ALL,            // components to be selected
 657             const_cast<char*>(mProfil[i].c_str()), // name of profil
 658             MED_GLOBAL,         // how to read data: MED_NO_PFLMOD,MED_COMPACT,MED_GLOBAL
 659             mEntity,            // type of entity (MED_NOEUD, MED_MAILLE, etc.)
 660             mGeom,              // type of geometry (TETRA10, etc.)
 661             mNumDT[i],          // time step iteration
 662             const_cast<char*>(mDTUnit[i].c_str()), // unit of time step
 663             mDT[i],             // time key
 664             mNumO[i]);          // order number
 665
 666         if (ret != 0) throw IOException("i/o error while writing field in MED file", __FILE__, __LINE__);
 667     }
 668 }
 669
 670 void    Field::writeMEDOptimized(std::vector<MeshDisPart*>* pParts, const char* pMeshName, GaussIndexList* pGaussList, int pGeomIdx, std::vector<med_int>& pFiles, std::map<std::string, Profil*>& pProfils)
 671 {
 672         med_err                                 ret;
 673         unsigned char*                  lValue = NULL;
 674         unsigned                                lGaussIdx = 0;
 675     bool                    lCreateField = false;
 676     med_int                 nbField;
 677     char                    lName[MED_TAILLE_NOM + 1];
 678     med_type_champ          type;
 679     char*                   comp = new char[mNumComponents * MED_TAILLE_PNOM + 1];
 680     char*                   unit = new char[mNumComponents * MED_TAILLE_PNOM + 1];
 681     set<med_int>*           idxSet;
 682     std::set< med_int>      lSetOfElt;
 683     bool                    completeProfil;
 684
 685     // For each time step.
 686         for (unsigned itTimeStep = 0;  itTimeStep < mNGauss.size(); ++itTimeStep)
 687         {
 688                 lGaussIdx = 0;
 689         // For each part.
 690                 for (unsigned itPart = 0 ; itPart < pParts->size() ; ++itPart)
 691                 {
 692                         // Skip empty part/mesh !
 693                         if ((*pParts)[itPart]->getMesh()->getNumberOfElements((eMeshType)pGeomIdx) == 0 &&
 694                                 !this->isFieldOnNodes())
 695                         {
 696                                 continue;
 697                         }
 698
 699                         // Get the index of the nodes or elements of this field. If the index is empty
 700             // the corresponding mesh was removed during domain split so we just skip it.
 701                         do
 702                         {
 703                                 if (this->isFieldOnNodes())
 704                                 {
 705                     // Field is on nodes.
 706                                         idxSet = &(*pGaussList)[lGaussIdx].second;
 707                                 }
 708                                 else
 709                                 {
 710                     // Field is on elements.
 711                                         idxSet = &(*pGaussList)[lGaussIdx].first[pGeomIdx];
 712                                 }
 713                                 lGaussIdx++;
 714                         }
 715                         while (idxSet->size() == 0 && lGaussIdx < pGaussList->size());
 716
 717             // Check if the gauss list is empty.
 718             if (lGaussIdx > pGaussList->size() || idxSet == 0)
 719                         {
 720                                 // We should never go here...
 721                                 throw IllegalStateException("Corrupted file !", __FILE__, __LINE__);
 722                         }
 723             // We check if the field already exists in this MED file.
 724             // I guess we could ignore this and always create the field but it seems
 725             // to create corrupted MED files.
 726             lCreateField = true;
 727             nbField = MEDnChamp(pFiles[itPart], 0);
 728             for (int i = 1; i <= nbField; ++i)
 729             {
 730                 ret = MEDchampInfo(pFiles[itPart], i, lName, &type, comp, unit, mNumComponents);
 731                 lName[MED_TAILLE_NOM] = '\0';
 732                 if (strncmp(mName, lName, MED_TAILLE_NOM) == 0)
 733                 {
 734                     lCreateField = false;
 735                     break;
 736                 }
 737             }
 738
 739             completeProfil = false;
 740             // If we have some profile.
 741             if (mProfil.size() > 0)
 742             {
 743                 // Try to find it.
 744                 if (pProfils.find(mProfil[itTimeStep]) != pProfils.end())
 745                 {
 746                     lSetOfElt.clear();
 747                     // Extract the subset of elements contained in the profil and the current part.
 748                     pProfils[mProfil[itTimeStep]]->filterSetOfElement(*idxSet, lSetOfElt);
 749                     // If this part has no nodes of elements in the profil we skip it.
 750                     if (lSetOfElt.size() == 0)
 751                     {
 752                         continue;
 753                     }
 754                     // Detect if the profil contains the whole field
 755                     if (lSetOfElt.size() == pProfils[mProfil[itTimeStep]]->getSet().size())
 756                     {
 757                         completeProfil = true;
 758                     }
 759                     // Replace the previous set.
 760                     idxSet = &lSetOfElt;
 761                 }
 762             }
 763
 764             // If the field doesn't exist we create it.
 765                         if (lCreateField == true)
 766             {
 767                 ret = MEDchampCr(
 768                     pFiles[itPart],
 769                     mName,                                                                       // name of the field
 770                     mType,                                                                       // type of data (MED_FLOAT64, MED_INT32, etc.)
 771                     const_cast<char*>(mStrComponent.c_str()),// name of components
 772                     const_cast<char*>(mStrUnit.c_str()),         // name of units
 773                     mNumComponents);                                             // number of components
 774             }
 775
 776             int nval;
 777             // If the profil contains the whole field, we dont need to copy it.
 778             if (completeProfil)
 779             {
 780                 nval = mNVal[itTimeStep];
 781                 lValue = mVal[itTimeStep];
 782             }
 783             else
 784             {
 785                 // Compute field number of values, size of one value etc.
 786                 nval = idxSet->size() * this->mNGauss[itTimeStep];
 787                 int sizeOfData = nval * mSizeOfType * mNumComponents;
 788                 int sizeOfOneData = mSizeOfType * mNumComponents * this->mNGauss[itTimeStep];
 789                 lValue = new unsigned char[sizeOfData + 1];
 790
 791                 // Copy the subset of values.
 792                 unsigned char* dest = lValue;
 793                 for (set<med_int>::iterator itSet = idxSet->begin() ; itSet != idxSet->end() ; ++itSet)
 794                 {
 795                     int indexElt = (*itSet);
 796                     unsigned char* src = mVal[itTimeStep] + (indexElt - 1) * sizeOfOneData;
 797                     memcpy(dest, src, sizeOfOneData);
 798                     dest += sizeOfOneData;
 799                 }
 800             }
 801
 802             // Write the field.
 803             ret =MEDchampEcr(
 804                 pFiles[itPart],
 805                 const_cast<char*>(pMeshName),// name of the mesh (first call to MEDchampEcr => name of reference)
 806                 mName,// name of the field
 807                 lValue,// data (= values)
 808                 MED_FULL_INTERLACE ,// data organization
 809                 nval,// number of values
 810                 const_cast<char*>(mGaussLoc[itTimeStep].c_str()),// name of Gauss reference
 811                 MED_ALL,// components to be selected
 812                 const_cast<char*>(mProfil[itTimeStep].c_str()),// name of profil
 813                 MED_COMPACT,// how to read data: MED_NO_PFLMOD,MED_COMPACT,MED_GLOBAL
 814                 mEntity,// type of entity (MED_NOEUD, MED_MAILLE, etc.)
 815                 mGeom,// type of geometry (TETRA10, etc.)
 816                 mNumDT[itTimeStep],// time step iteration
 817                 const_cast<char*>(mDTUnit[itTimeStep].c_str()), // unit of time step
 818                 mDT[itTimeStep],// time key
 819                 mNumO[itTimeStep]);// order number
 820
 821             // If we use the field's data we dont need to delete it.
 822             if (!completeProfil)
 823             {
 824                 delete[] lValue;
 825             }
 826
 827                         if (ret != 0)
 828                         {
 829                                 throw IOException("i/o error while writing field in MED file", __FILE__, __LINE__);
 830                         }
 831
 832                 }
 833         }
 834 }
 835
 836 ostream& operator<<(ostream& pOs, Field& pF)
 837 {
 838     char strEntity[16];
 839     switch (pF.mEntity)
 840     {
 841         case MED_MAILLE: strcpy(strEntity, "MED_MAILLE"); break;
 842         case MED_FACE:   strcpy(strEntity, "MED_FACE"); break;
 843         case MED_ARETE:  strcpy(strEntity, "MED_ARETE"); break;
 844         case MED_NOEUD:  strcpy(strEntity, "MED_NOEUD"); break;
 845         default:         strcpy(strEntity, "UNKNOWN"); break;
 846     }
 847
 848     char strType[16];
 849     switch (pF.mType)
 850     {
 851         case MED_FLOAT64: strcpy(strType, "MED_FLOAT64"); break;
 852         case MED_INT32:   strcpy(strType, "MED_INT32"); break;
 853         case MED_INT64:   strcpy(strType, "MED_INT64"); break;
 854         case MED_INT:     strcpy(strType, "MED_INT"); break;
 855         default:          strcpy(strType, "UNKNOWN"); break;
 856     }
 857
 858     pOs << "Field: " << endl;
 859     pOs << "    Name       =|" << pF.mName << "|" << endl;
 860     pOs << "    Entity     =" << strEntity << endl;
 861     pOs << "    Geom       =" << pF.mGeom << endl;
 862     pOs << "    Type       =" << strType << " (size=" << pF.mSizeOfType << ")" << endl;
 863     pOs << "    #Components=" << pF.mNumComponents << endl;
 864     pOs << "        Name component=|" << pF.mStrComponent << "|" << endl;
 865     pOs << "        Unit component=|" << pF.mStrUnit << "|" << endl;
 866     pOs << "    #Time steps=" << pF.mNGauss.size() << endl;
 867
 868     for (unsigned itTimeStep = 0 ; itTimeStep < pF.mNGauss.size() ; itTimeStep++)
 869     {
 870         pOs << "        Time=" << pF.mDT[itTimeStep];
 871         pOs << "  it=" << pF.mNumDT[itTimeStep];
 872         pOs << "  order=" << pF.mNumO[itTimeStep];
 873         pOs << "  #gauss=" << pF.mNGauss[itTimeStep];
 874         pOs << "  #val=" << pF.mNVal[itTimeStep];
 875         pOs << "  sizeof_val=" << pF.mSizeOfData[itTimeStep];
 876         pOs << "  gauss_loc=|" << ((pF.mGaussLoc[itTimeStep].size() == 0)?"NONE":pF.mGaussLoc[itTimeStep]) << "|  size=" << pF.mGaussLoc[itTimeStep].size();
 877         pOs << "  profil=|" << ((pF.mProfil[itTimeStep].size() == 0)?"NONE":pF.mProfil[itTimeStep]) << "|  size=" << pF.mProfil[itTimeStep].size() << endl;
 878
 879         if (pF.mFlagPrintAll)
 880         {
 881             cout << "    Values: ";
 882             switch (pF.mType)
 883             {
 884             case MED_FLOAT64:
 885                 {
 886                     med_float* src = reinterpret_cast<med_float*>(pF.mVal[itTimeStep]);
 887                     for (int itVal = 0 ; itVal < pF.mNVal[itTimeStep] * pF.mNumComponents ; itVal++)
 888                     {
 889                         cout << src[itVal] << " ";
 890                     }
 891                 }
 892                 break;
 893             case MED_INT:
 894             case MED_INT32:
 895                 {
 896                     med_int* src = reinterpret_cast<med_int*>(pF.mVal[itTimeStep]);
 897                     for (int itVal = 0 ; itVal < pF.mNVal[itTimeStep] * pF.mNumComponents ; itVal++)
 898                     {
 899                         cout << src[itVal] << " ";
 900                     }
 901                 }
 902                 break;
 903             case MED_INT64:
 904                 // not yet implemented
 905                 throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
 906             default:
 907                 // should not be there
 908                 throw IllegalStateException("should not be there", __FILE__, __LINE__);
 909             }
 910             cout << endl;
 911         }
 912
 913     }
 914
 915     return pOs;
 916 }
 917
 918
 919 } // namespace  multipr
 920
 921 // EOF