src/Plugins/MEDWriter/plugin/MEDWriterIO/VTKToMEDMem.cxx

   1 // Copyright (C) 2017-2020  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, or (at your option) any later version.
   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 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "VTKToMEDMem.h"
  22
  23 #include "vtkAdjacentVertexIterator.h"
  24 #include "vtkIntArray.h"
  25 #include "vtkLongArray.h"
  26 #include "vtkCellData.h"
  27 #include "vtkPointData.h"
  28 #include "vtkFloatArray.h"
  29 #include "vtkCellArray.h"
  30
  31 #include "vtkStreamingDemandDrivenPipeline.h"
  32 #include "vtkInformationDataObjectMetaDataKey.h"
  33 #include "vtkUnstructuredGrid.h"
  34 #include "vtkMultiBlockDataSet.h"
  35 #include "vtkRectilinearGrid.h"
  36 #include "vtkInformationStringKey.h"
  37 #include "vtkAlgorithmOutput.h"
  38 #include "vtkObjectFactory.h"
  39 #include "vtkMutableDirectedGraph.h"
  40 #include "vtkMultiBlockDataSet.h"
  41 #include "vtkPolyData.h"
  42 #include "vtkDataSet.h"
  43 #include "vtkInformationVector.h"
  44 #include "vtkInformation.h"
  45 #include "vtkDataArraySelection.h"
  46 #include "vtkTimeStamp.h"
  47 #include "vtkInEdgeIterator.h"
  48 #include "vtkInformationDataObjectKey.h"
  49 #include "vtkExecutive.h"
  50 #include "vtkVariantArray.h"
  51 #include "vtkStringArray.h"
  52 #include "vtkDoubleArray.h"
  53 #include "vtkCharArray.h"
  54 #include "vtkUnsignedCharArray.h"
  55 #include "vtkDataSetAttributes.h"
  56 #include "vtkDemandDrivenPipeline.h"
  57 #include "vtkDataObjectTreeIterator.h"
  58 #include "vtkWarpScalar.h"
  59
  60 #include <map>
  61 #include <deque>
  62 #include <sstream>
  63 #include <cstring>
  64
  65 using VTKToMEDMem::Grp;
  66 using VTKToMEDMem::Fam;
  67
  68 using MEDCoupling::MEDFileData;
  69 using MEDCoupling::MEDFileMesh;
  70 using MEDCoupling::MEDFileCMesh;
  71 using MEDCoupling::MEDFileUMesh;
  72 using MEDCoupling::MEDFileFields;
  73 using MEDCoupling::MEDFileMeshes;
  74
  75 using MEDCoupling::MEDFileInt32Field1TS;
  76 using MEDCoupling::MEDFileInt64Field1TS;
  77 using MEDCoupling::MEDFileField1TS;
  78 using MEDCoupling::MEDFileIntFieldMultiTS;
  79 using MEDCoupling::MEDFileFieldMultiTS;
  80 using MEDCoupling::MEDFileAnyTypeFieldMultiTS;
  81 using MEDCoupling::DataArray;
  82 using MEDCoupling::DataArrayInt32;
  83 using MEDCoupling::DataArrayInt64;
  84 using MEDCoupling::DataArrayFloat;
  85 using MEDCoupling::DataArrayDouble;
  86 using MEDCoupling::MEDCouplingMesh;
  87 using MEDCoupling::MEDCouplingUMesh;
  88 using MEDCoupling::MEDCouplingCMesh;
  89 using MEDCoupling::MEDCouplingFieldDouble;
  90 using MEDCoupling::MEDCouplingFieldFloat;
  91 using MEDCoupling::MEDCouplingFieldInt;
  92 using MEDCoupling::MCAuto;
  93 using MEDCoupling::Traits;
  94 using MEDCoupling::MLFieldTraits;
  95
  96 ///////////////////
  97
  98 Fam::Fam(const std::string& name)
  99 {
 100   static const char ZE_SEP[]="@@][@@";
 101   std::size_t pos(name.find(ZE_SEP));
 102   std::string name0(name.substr(0,pos)),name1(name.substr(pos+strlen(ZE_SEP)));
 103   std::istringstream iss(name1);
 104   iss >> _id;
 105   _name=name0;
 106 }
 107
 108 ///////////////////
 109
 110 #include "VTKMEDTraits.hxx"
 111
 112 std::map<int,int> ComputeMapOfType()
 113 {
 114   std::map<int,int> ret;
 115   int nbOfTypesInMC(sizeof(MEDCOUPLING2VTKTYPETRADUCER)/sizeof( decltype(MEDCOUPLING2VTKTYPETRADUCER[0]) ));
 116   for(int i=0;i<nbOfTypesInMC;i++)
 117     {
 118       auto vtkId(MEDCOUPLING2VTKTYPETRADUCER[i]);
 119       if(vtkId!=MEDCOUPLING2VTKTYPETRADUCER_NONE)
 120         ret[vtkId]=i;
 121     }
 122   return ret;
 123 }
 124
 125 std::string GetMeshNameWithContext(const std::vector<int>& context)
 126 {
 127   static const char DFT_MESH_NAME[]="Mesh";
 128   if(context.empty())
 129     return DFT_MESH_NAME;
 130   std::ostringstream oss; oss << DFT_MESH_NAME;
 131   for(std::vector<int>::const_iterator it=context.begin();it!=context.end();it++)
 132     oss << "_" << *it;
 133   return oss.str();
 134 }
 135
 136 DataArrayIdType *ConvertVTKArrayToMCArrayInt(vtkDataArray *data)
 137 {
 138   if(!data)
 139     throw MZCException("ConvertVTKArrayToMCArrayInt : internal error !");
 140   int nbTuples(data->GetNumberOfTuples()),nbComp(data->GetNumberOfComponents());
 141   std::size_t nbElts(nbTuples*nbComp);
 142   MCAuto<DataArrayIdType> ret(DataArrayIdType::New());
 143   ret->alloc(nbTuples,nbComp);
 144   for(int i=0;i<nbComp;i++)
 145     {
 146       const char *comp(data->GetComponentName(i));
 147       if(comp)
 148         ret->setInfoOnComponent(i,comp);
 149     }
 150   mcIdType *ptOut(ret->getPointer());
 151   vtkIntArray *d0(vtkIntArray::SafeDownCast(data));
 152   if(d0)
 153     {
 154       const int *pt(d0->GetPointer(0));
 155       std::copy(pt,pt+nbElts,ptOut);
 156       return ret.retn();
 157     }
 158   vtkLongArray *d1(vtkLongArray::SafeDownCast(data));
 159   if(d1)
 160     {
 161       const long *pt(d1->GetPointer(0));
 162       std::copy(pt,pt+nbElts,ptOut);
 163       return ret.retn();
 164     }
 165   vtkIdTypeArray *d3(vtkIdTypeArray::SafeDownCast(data));
 166   if(d3)
 167     {
 168       const vtkIdType *pt(d3->GetPointer(0));
 169       std::copy(pt,pt+nbElts,ptOut);
 170       return ret.retn();
 171     }
 172   vtkUnsignedCharArray *d2(vtkUnsignedCharArray::SafeDownCast(data));
 173   if(d2)
 174     {
 175       const unsigned char *pt(d2->GetPointer(0));
 176       std::copy(pt,pt+nbElts,ptOut);
 177       return ret.retn();
 178     }
 179   std::ostringstream oss;
 180   oss << "ConvertVTKArrayToMCArrayInt : unrecognized array \"" << typeid(*data).name() << "\" type !";
 181   throw MZCException(oss.str());
 182 }
 183
 184 template<class T>
 185 typename Traits<T>::ArrayType *ConvertVTKArrayToMCArrayDouble(vtkDataArray *data)
 186 {
 187   if(!data)
 188     throw MZCException("ConvertVTKArrayToMCArrayDouble : internal error !");
 189   int nbTuples(data->GetNumberOfTuples()),nbComp(data->GetNumberOfComponents());
 190   std::size_t nbElts(nbTuples*nbComp);
 191   MCAuto< typename Traits<T>::ArrayType > ret(Traits<T>::ArrayType::New());
 192   ret->alloc(nbTuples,nbComp);
 193   for(int i=0;i<nbComp;i++)
 194     {
 195       const char *comp(data->GetComponentName(i));
 196       if(comp)
 197         ret->setInfoOnComponent(i,comp);
 198       else
 199         {
 200           if(nbComp>1 && nbComp<=3)
 201             {
 202               char tmp[2];
 203               tmp[0]=(char)('X'+i); tmp[1]='\0';
 204               ret->setInfoOnComponent(i,tmp);
 205             }
 206         }
 207     }
 208   T *ptOut(ret->getPointer());
 209   typename MEDFileVTKTraits<T>::VtkType *d0(MEDFileVTKTraits<T>::VtkType::SafeDownCast(data));
 210   if(d0)
 211     {
 212       const T *pt(d0->GetPointer(0));
 213       for(std::size_t i=0;i<nbElts;i++)
 214         ptOut[i]=(T)pt[i];
 215       return ret.retn();
 216     }
 217   std::ostringstream oss;
 218   oss << "ConvertVTKArrayToMCArrayDouble : unrecognized array \"" << data->GetClassName() << "\" type !";
 219   throw MZCException(oss.str());
 220 }
 221
 222 DataArrayDouble *ConvertVTKArrayToMCArrayDoubleForced(vtkDataArray *data)
 223 {
 224   if(!data)
 225     throw MZCException("ConvertVTKArrayToMCArrayDoubleForced : internal error 0 !");
 226   vtkFloatArray *d0(vtkFloatArray::SafeDownCast(data));
 227   if(d0)
 228     {
 229       MCAuto<DataArrayFloat> ret(ConvertVTKArrayToMCArrayDouble<float>(data));
 230       MCAuto<DataArrayDouble> ret2(ret->convertToDblArr());
 231       return ret2.retn();
 232     }
 233   vtkDoubleArray *d1(vtkDoubleArray::SafeDownCast(data));
 234   if(d1)
 235     return ConvertVTKArrayToMCArrayDouble<double>(data);
 236   throw MZCException("ConvertVTKArrayToMCArrayDoubleForced : unrecognized type of data for double !");
 237 }
 238
 239 DataArray *ConvertVTKArrayToMCArray(vtkDataArray *data)
 240 {
 241   if(!data)
 242     throw MZCException("ConvertVTKArrayToMCArray : internal error !");
 243   vtkFloatArray *d0(vtkFloatArray::SafeDownCast(data));
 244   if(d0)
 245     return ConvertVTKArrayToMCArrayDouble<float>(data);
 246   vtkDoubleArray *d1(vtkDoubleArray::SafeDownCast(data));
 247   if(d1)
 248     return ConvertVTKArrayToMCArrayDouble<double>(data);
 249   vtkIntArray *d2(vtkIntArray::SafeDownCast(data));
 250   vtkLongArray *d3(vtkLongArray::SafeDownCast(data));
 251   vtkUnsignedCharArray *d4(vtkUnsignedCharArray::SafeDownCast(data));
 252   vtkIdTypeArray *d5(vtkIdTypeArray::SafeDownCast(data));
 253   if(d2 || d3 || d4 || d5)
 254     return ConvertVTKArrayToMCArrayInt(data);
 255   std::ostringstream oss;
 256   oss << "ConvertVTKArrayToMCArray : unrecognized array \"" << typeid(*data).name() << "\" type !";
 257   throw MZCException(oss.str());
 258 }
 259
 260 MEDCouplingUMesh *BuildMeshFromCellArray(vtkCellArray *ca, DataArrayDouble *coords, int meshDim, INTERP_KERNEL::NormalizedCellType type)
 261 {
 262   MCAuto<MEDCouplingUMesh> subMesh(MEDCouplingUMesh::New("",meshDim));
 263   subMesh->setCoords(coords); subMesh->allocateCells();
 264   int nbCells(ca->GetNumberOfCells());
 265   if(nbCells==0)
 266     return 0;
 267   vtkIdType nbEntries(ca->GetNumberOfConnectivityEntries());
 268   const vtkIdType *conn(ca->GetData()->GetPointer(0));
 269   for(int i=0;i<nbCells;i++)
 270     {
 271       mcIdType sz(ToIdType(*conn++));
 272       std::vector<mcIdType> conn2(sz);
 273       for(int jj=0;jj<sz;jj++)
 274         conn2[jj]=ToIdType(conn[jj]);
 275       subMesh->insertNextCell(type,sz,&conn2[0]);
 276       conn+=sz;
 277     }
 278   return subMesh.retn();
 279 }
 280
 281 MEDCouplingUMesh *BuildMeshFromCellArrayTriangleStrip(vtkCellArray *ca, DataArrayDouble *coords, MCAuto<DataArrayIdType>& ids)
 282 {
 283   MCAuto<MEDCouplingUMesh> subMesh(MEDCouplingUMesh::New("",2));
 284   subMesh->setCoords(coords); subMesh->allocateCells();
 285   int nbCells(ca->GetNumberOfCells());
 286   if(nbCells==0)
 287     return 0;
 288   vtkIdType nbEntries(ca->GetNumberOfConnectivityEntries());
 289   const vtkIdType *conn(ca->GetData()->GetPointer(0));
 290   ids=DataArrayIdType::New() ; ids->alloc(0,1);
 291   for(int i=0;i<nbCells;i++)
 292     {
 293       int sz(*conn++);
 294       int nbTri(sz-2);
 295       if(nbTri>0)
 296         {
 297           for(int j=0;j<nbTri;j++,conn++)
 298             {
 299               mcIdType conn2[3]; conn2[0]=conn[0] ; conn2[1]=conn[1] ; conn2[2]=conn[2];
 300               subMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn2);
 301               ids->pushBackSilent(i);
 302             }
 303         }
 304       else
 305         {
 306           std::ostringstream oss; oss << "BuildMeshFromCellArrayTriangleStrip : on cell #" << i << " the triangle stip looks bab !";
 307           throw MZCException(oss.str());
 308         }
 309       conn+=sz;
 310     }
 311   return subMesh.retn();
 312 }
 313
 314 class MicroField
 315 {
 316 public:
 317   MicroField(const MCAuto<MEDCouplingUMesh>& m, const std::vector<MCAuto<DataArray> >& cellFs):_m(m),_cellFs(cellFs) { }
 318   MicroField(const std::vector< MicroField >& vs);
 319   void setNodeFields(const std::vector<MCAuto<DataArray> >& nf) { _nodeFs=nf; }
 320   MCAuto<MEDCouplingUMesh> getMesh() const { return _m; }
 321   std::vector<MCAuto<DataArray> > getCellFields() const { return _cellFs; }
 322 private:
 323   MCAuto<MEDCouplingUMesh> _m;
 324   std::vector<MCAuto<DataArray> > _cellFs;
 325   std::vector<MCAuto<DataArray> > _nodeFs;
 326 };
 327
 328 MicroField::MicroField(const std::vector< MicroField >& vs)
 329 {
 330   std::size_t sz(vs.size());
 331   std::vector<const MEDCouplingUMesh *> vs2(sz);
 332   std::vector< std::vector< MCAuto<DataArray> > > arrs2(sz);
 333   int nbElts(-1);
 334   for(std::size_t ii=0;ii<sz;ii++)
 335     {
 336       vs2[ii]=vs[ii].getMesh();
 337       arrs2[ii]=vs[ii].getCellFields();
 338       if(nbElts<0)
 339         nbElts=(int)arrs2[ii].size();
 340       else
 341         if((int)arrs2[ii].size()!=nbElts)
 342           throw MZCException("MicroField cstor : internal error !");
 343     }
 344   _cellFs.resize(nbElts);
 345   for(int ii=0;ii<nbElts;ii++)
 346     {
 347       std::vector<const DataArray *> arrsTmp(sz);
 348       for(std::size_t jj=0;jj<sz;jj++)
 349         {
 350           arrsTmp[jj]=arrs2[jj][ii];
 351         }
 352       _cellFs[ii]=DataArray::Aggregate(arrsTmp);
 353     }
 354   _m=MEDCouplingUMesh::MergeUMeshesOnSameCoords(vs2);
 355 }
 356
 357 template<class T>
 358 void AppendToFields(MEDCoupling::TypeOfField tf, MEDCouplingMesh *mesh, const DataArrayIdType *n2oPtr, typename MEDFileVTKTraits<T>::MCType *dadPtr, MEDFileFields *fs, double timeStep, int tsId)
 359 {
 360   std::string fieldName(dadPtr->getName());
 361   MCAuto< typename Traits<T>::FieldType > f(Traits<T>::FieldType::New(tf));
 362   f->setTime(timeStep,tsId,0);
 363   {
 364     std::string fieldNameForChuckNorris(MEDCoupling::MEDFileAnyTypeField1TSWithoutSDA::FieldNameToMEDFileConvention(fieldName));
 365     f->setName(fieldNameForChuckNorris);
 366   }
 367   if(!n2oPtr)
 368     f->setArray(dadPtr);
 369   else
 370     {
 371       MCAuto< typename Traits<T>::ArrayType > dad2(dadPtr->selectByTupleId(n2oPtr->begin(),n2oPtr->end()));
 372       f->setArray(dad2);
 373     }
 374   f->setMesh(mesh);
 375   MCAuto< typename MLFieldTraits<T>::FMTSType > fmts(MLFieldTraits<T>::FMTSType::New());
 376   MCAuto< typename MLFieldTraits<T>::F1TSType > f1ts(MLFieldTraits<T>::F1TSType::New());
 377   f1ts->setFieldNoProfileSBT(f);
 378   fmts->pushBackTimeStep(f1ts);
 379   fs->pushField(fmts);
 380 }
 381
 382 void AppendMCFieldFrom(MEDCoupling::TypeOfField tf, MEDCouplingMesh *mesh, MEDFileData *mfd, MCAuto<DataArray> da, const DataArrayIdType *n2oPtr, double timeStep, int tsId)
 383 {
 384   static const char FAMFIELD_FOR_CELLS[]="FamilyIdCell";
 385   static const char FAMFIELD_FOR_NODES[]="FamilyIdNode";
 386   if(!da || !mesh || !mfd)
 387     throw MZCException("AppendMCFieldFrom : internal error !");
 388   MEDFileFields *fs(mfd->getFields());
 389   MEDFileMeshes *ms(mfd->getMeshes());
 390   if(!fs || !ms)
 391     throw MZCException("AppendMCFieldFrom : internal error 2 !");
 392   MCAuto<DataArrayDouble> dad(MEDCoupling::DynamicCast<DataArray,DataArrayDouble>(da));
 393   if(dad.isNotNull())
 394     {
 395       AppendToFields<double>(tf,mesh,n2oPtr,dad,fs,timeStep,tsId);
 396       return ;
 397     }
 398   MCAuto<DataArrayFloat> daf(MEDCoupling::DynamicCast<DataArray,DataArrayFloat>(da));
 399   if(daf.isNotNull())
 400     {
 401       AppendToFields<float>(tf,mesh,n2oPtr,daf,fs,timeStep,tsId);
 402       return ;
 403     }
 404   MCAuto<DataArrayInt> dai(MEDCoupling::DynamicCast<DataArray,DataArrayInt>(da));
 405   MCAuto<DataArrayIdType> daId(MEDCoupling::DynamicCast<DataArray,DataArrayIdType>(da));
 406   if(dai.isNotNull() || daId.isNotNull())
 407     {
 408       std::string fieldName(da->getName());
 409       if((fieldName!=FAMFIELD_FOR_CELLS || tf!=MEDCoupling::ON_CELLS) && (fieldName!=FAMFIELD_FOR_NODES || tf!=MEDCoupling::ON_NODES))
 410         {
 411           if(!dai)
 412             throw MZCException("AppendMCFieldFrom : internal error 3 (not int32) !");
 413           AppendToFields<int>(tf,mesh,n2oPtr,dai,fs,timeStep,tsId);
 414           return ;
 415         }
 416       else if(fieldName==FAMFIELD_FOR_CELLS && tf==MEDCoupling::ON_CELLS)
 417         {
 418           MEDFileMesh *mm(ms->getMeshWithName(mesh->getName()));
 419           if(!mm)
 420             throw MZCException("AppendMCFieldFrom : internal error 3 !");
 421           if(!daId)
 422             throw MZCException("AppendMCFieldFrom : internal error 3 (not mcIdType) !");
 423           if(!n2oPtr)
 424             mm->setFamilyFieldArr(mesh->getMeshDimension()-mm->getMeshDimension(),daId);
 425           else
 426             {
 427               MCAuto<DataArrayIdType> dai2(daId->selectByTupleId(n2oPtr->begin(),n2oPtr->end()));
 428               mm->setFamilyFieldArr(mesh->getMeshDimension()-mm->getMeshDimension(),dai2);
 429             }
 430         }
 431       else if(fieldName==FAMFIELD_FOR_NODES || tf==MEDCoupling::ON_NODES)
 432         {
 433           MEDFileMesh *mm(ms->getMeshWithName(mesh->getName()));
 434           if(!mm)
 435             throw MZCException("AppendMCFieldFrom : internal error 4 !");
 436           if(!daId)
 437             throw MZCException("AppendMCFieldFrom : internal error 4 (not mcIdType) !");
 438           if(!n2oPtr)
 439             mm->setFamilyFieldArr(1,daId);
 440           else
 441             {
 442               MCAuto<DataArrayIdType> dai2(daId->selectByTupleId(n2oPtr->begin(),n2oPtr->end()));
 443               mm->setFamilyFieldArr(1,dai2);
 444             }
 445         }
 446     }
 447 }
 448
 449 void PutAtLevelDealOrder(MEDFileData *mfd, int meshDimRel, const MicroField& mf, double timeStep, int tsId)
 450 {
 451   if(!mfd)
 452     throw MZCException("PutAtLevelDealOrder : internal error !");
 453   MEDFileMesh *mm(mfd->getMeshes()->getMeshAtPos(0));
 454   MEDFileUMesh *mmu(dynamic_cast<MEDFileUMesh *>(mm));
 455   if(!mmu)
 456     throw MZCException("PutAtLevelDealOrder : internal error 2 !");
 457   MCAuto<MEDCouplingUMesh> mesh(mf.getMesh());
 458   mesh->setName(mfd->getMeshes()->getMeshAtPos(0)->getName());
 459   MCAuto<DataArrayIdType> o2n(mesh->sortCellsInMEDFileFrmt());
 460   const DataArrayIdType *o2nPtr(o2n);
 461   MCAuto<DataArrayIdType> n2o;
 462   mmu->setMeshAtLevel(meshDimRel,mesh);
 463   const DataArrayIdType *n2oPtr(0);
 464   if(o2n)
 465     {
 466       n2o=o2n->invertArrayO2N2N2O(mesh->getNumberOfCells());
 467       n2oPtr=n2o;
 468       if(n2oPtr && n2oPtr->isIota(mesh->getNumberOfCells()))
 469         n2oPtr=0;
 470       if(n2oPtr)
 471         mm->setRenumFieldArr(meshDimRel,n2o);
 472     }
 473   //
 474   std::vector<MCAuto<DataArray> > cells(mf.getCellFields());
 475   for(std::vector<MCAuto<DataArray> >::const_iterator it=cells.begin();it!=cells.end();it++)
 476     {
 477       MCAuto<DataArray> da(*it);
 478       AppendMCFieldFrom(MEDCoupling::ON_CELLS,mesh,mfd,da,n2oPtr,timeStep,tsId);
 479     }
 480 }
 481
 482 void AssignSingleGTMeshes(MEDFileData *mfd, const std::vector< MicroField >& ms, double timeStep, int tsId)
 483 {
 484   if(!mfd)
 485     throw MZCException("AssignSingleGTMeshes : internal error !");
 486   MEDFileMesh *mm0(mfd->getMeshes()->getMeshAtPos(0));
 487   MEDFileUMesh *mm(dynamic_cast<MEDFileUMesh *>(mm0));
 488   if(!mm)
 489     throw MZCException("AssignSingleGTMeshes : internal error 2 !");
 490   int meshDim(-std::numeric_limits<int>::max());
 491   std::map<int, std::vector< MicroField > > ms2;
 492   for(std::vector< MicroField >::const_iterator it=ms.begin();it!=ms.end();it++)
 493     {
 494       const MEDCouplingUMesh *elt((*it).getMesh());
 495       if(elt)
 496         {
 497           int myMeshDim(elt->getMeshDimension());
 498           meshDim=std::max(meshDim,myMeshDim);
 499           ms2[myMeshDim].push_back(*it);
 500         }
 501     }
 502   if(ms2.empty())
 503     return ;
 504   for(std::map<int, std::vector< MicroField > >::const_iterator it=ms2.begin();it!=ms2.end();it++)
 505     {
 506       const std::vector< MicroField >& vs((*it).second);
 507       if(vs.size()==1)
 508         {
 509           PutAtLevelDealOrder(mfd,(*it).first-meshDim,vs[0],timeStep,tsId);
 510         }
 511       else
 512         {
 513           MicroField merge(vs);
 514           PutAtLevelDealOrder(mfd,(*it).first-meshDim,merge,timeStep,tsId);
 515         }
 516     }
 517 }
 518
 519 DataArrayDouble *BuildCoordsFrom(vtkPointSet *ds)
 520 {
 521   if(!ds)
 522     throw MZCException("BuildCoordsFrom : internal error !");
 523   vtkPoints *pts(ds->GetPoints());
 524   if(!pts)
 525     throw MZCException("BuildCoordsFrom : internal error 2 !");
 526   vtkDataArray *data(pts->GetData());
 527   if(!data)
 528     throw MZCException("BuildCoordsFrom : internal error 3 !");
 529   return ConvertVTKArrayToMCArrayDoubleForced(data);
 530 }
 531
 532 void AddNodeFields(MEDFileData *mfd, vtkDataSetAttributes *dsa, double timeStep, int tsId)
 533 {
 534   if(!mfd || !dsa)
 535     throw MZCException("AddNodeFields : internal error !");
 536   MEDFileMesh *mm(mfd->getMeshes()->getMeshAtPos(0));
 537   MEDFileUMesh *mmu(dynamic_cast<MEDFileUMesh *>(mm));
 538   if(!mmu)
 539     throw MZCException("AddNodeFields : internal error 2 !");
 540   MCAuto<MEDCouplingUMesh> mesh;
 541   if(!mmu->getNonEmptyLevels().empty())
 542     mesh=mmu->getMeshAtLevel(0);
 543   else
 544     {
 545       mesh=MEDCouplingUMesh::Build0DMeshFromCoords(mmu->getCoords());
 546       mesh->setName(mmu->getName());
 547     }
 548   int nba(dsa->GetNumberOfArrays());
 549   for(int i=0;i<nba;i++)
 550     {
 551       vtkDataArray *arr(dsa->GetArray(i));
 552       const char *name(arr->GetName());
 553       if(!arr)
 554         continue;
 555       MCAuto<DataArray> da(ConvertVTKArrayToMCArray(arr));
 556       da->setName(name);
 557       AppendMCFieldFrom(MEDCoupling::ON_NODES,mesh,mfd,da,NULL,timeStep,tsId);
 558     }
 559 }
 560
 561 std::vector<MCAuto<DataArray> > AddPartFields(const DataArrayIdType *part, vtkDataSetAttributes *dsa)
 562 {
 563   std::vector< MCAuto<DataArray> > ret;
 564   if(!dsa)
 565     return ret;
 566   int nba(dsa->GetNumberOfArrays());
 567   for(int i=0;i<nba;i++)
 568     {
 569       vtkDataArray *arr(dsa->GetArray(i));
 570       if(!arr)
 571         continue;
 572       const char *name(arr->GetName());
 573       int nbCompo(arr->GetNumberOfComponents());
 574       vtkIdType nbTuples(arr->GetNumberOfTuples());
 575       MCAuto<DataArray> mcarr(ConvertVTKArrayToMCArray(arr));
 576       if(part)
 577         mcarr=mcarr->selectByTupleId(part->begin(),part->end());
 578       mcarr->setName(name);
 579       ret.push_back(mcarr);
 580     }
 581   return ret;
 582 }
 583
 584 std::vector<MCAuto<DataArray> > AddPartFields2(int bg, int end, vtkDataSetAttributes *dsa)
 585 {
 586   std::vector< MCAuto<DataArray> > ret;
 587   if(!dsa)
 588     return ret;
 589   int nba(dsa->GetNumberOfArrays());
 590   for(int i=0;i<nba;i++)
 591     {
 592       vtkDataArray *arr(dsa->GetArray(i));
 593       if(!arr)
 594         continue;
 595       const char *name(arr->GetName());
 596       int nbCompo(arr->GetNumberOfComponents());
 597       vtkIdType nbTuples(arr->GetNumberOfTuples());
 598       MCAuto<DataArray> mcarr(ConvertVTKArrayToMCArray(arr));
 599       mcarr=mcarr->selectByTupleIdSafeSlice(bg,end,1);
 600       mcarr->setName(name);
 601       ret.push_back(mcarr);
 602     }
 603   return ret;
 604 }
 605
 606 void ConvertFromRectilinearGrid(MEDFileData *ret, vtkRectilinearGrid *ds, const std::vector<int>& context, double timeStep, int tsId)
 607 {
 608   if(!ds || !ret)
 609     throw MZCException("ConvertFromRectilinearGrid : internal error !");
 610   //
 611   MCAuto<MEDFileMeshes> meshes(MEDFileMeshes::New());
 612   ret->setMeshes(meshes);
 613   MCAuto<MEDFileFields> fields(MEDFileFields::New());
 614   ret->setFields(fields);
 615   //
 616   MCAuto<MEDFileCMesh> cmesh(MEDFileCMesh::New());
 617   meshes->pushMesh(cmesh);
 618   MCAuto<MEDCouplingCMesh> cmeshmc(MEDCouplingCMesh::New());
 619   vtkDataArray *cx(ds->GetXCoordinates()),*cy(ds->GetYCoordinates()),*cz(ds->GetZCoordinates());
 620   if(cx)
 621     {
 622       MCAuto<DataArrayDouble> arr(ConvertVTKArrayToMCArrayDoubleForced(cx));
 623       cmeshmc->setCoordsAt(0,arr);
 624     }
 625   if(cy)
 626     {
 627       MCAuto<DataArrayDouble> arr(ConvertVTKArrayToMCArrayDoubleForced(cy));
 628       cmeshmc->setCoordsAt(1,arr);
 629     }
 630   if(cz)
 631     {
 632       MCAuto<DataArrayDouble> arr(ConvertVTKArrayToMCArrayDoubleForced(cz));
 633       cmeshmc->setCoordsAt(2,arr);
 634     }
 635   std::string meshName(GetMeshNameWithContext(context));
 636   cmeshmc->setName(meshName);
 637   cmesh->setMesh(cmeshmc);
 638   std::vector<MCAuto<DataArray> > cellFs(AddPartFields(0,ds->GetCellData()));
 639   for(std::vector<MCAuto<DataArray> >::const_iterator it=cellFs.begin();it!=cellFs.end();it++)
 640     {
 641       MCAuto<DataArray> da(*it);
 642       AppendMCFieldFrom(MEDCoupling::ON_CELLS,cmeshmc,ret,da,NULL,timeStep,tsId);
 643     }
 644   std::vector<MCAuto<DataArray> > nodeFs(AddPartFields(0,ds->GetPointData()));
 645   for(std::vector<MCAuto<DataArray> >::const_iterator it=nodeFs.begin();it!=nodeFs.end();it++)
 646     {
 647       MCAuto<DataArray> da(*it);
 648       AppendMCFieldFrom(MEDCoupling::ON_NODES,cmeshmc,ret,da,NULL,timeStep,tsId);
 649     }
 650 }
 651
 652 void ConvertFromPolyData(MEDFileData *ret, vtkPolyData *ds, const std::vector<int>& context, double timeStep, int tsId)
 653 {
 654   if(!ds || !ret)
 655     throw MZCException("ConvertFromPolyData : internal error !");
 656   //
 657   MCAuto<MEDFileMeshes> meshes(MEDFileMeshes::New());
 658   ret->setMeshes(meshes);
 659   MCAuto<MEDFileFields> fields(MEDFileFields::New());
 660   ret->setFields(fields);
 661   //
 662   MCAuto<MEDFileUMesh> umesh(MEDFileUMesh::New());
 663   meshes->pushMesh(umesh);
 664   MCAuto<DataArrayDouble> coords(BuildCoordsFrom(ds));
 665   umesh->setCoords(coords);
 666   umesh->setName(GetMeshNameWithContext(context));
 667   //
 668   int offset(0);
 669   std::vector< MicroField > ms;
 670   vtkCellArray *cd(ds->GetVerts());
 671   if(cd)
 672     {
 673       MCAuto<MEDCouplingUMesh> subMesh(BuildMeshFromCellArray(cd,coords,0,INTERP_KERNEL::NORM_POINT1));
 674       if((const MEDCouplingUMesh *)subMesh)
 675         {
 676           std::vector<MCAuto<DataArray> > cellFs(AddPartFields2(offset,offset+subMesh->getNumberOfCells(),ds->GetCellData()));
 677           offset+=subMesh->getNumberOfCells();
 678           ms.push_back(MicroField(subMesh,cellFs));
 679         }
 680     }
 681   vtkCellArray *cc(ds->GetLines());
 682   if(cc)
 683     {
 684       MCAuto<MEDCouplingUMesh> subMesh;
 685       try
 686         {
 687           subMesh=BuildMeshFromCellArray(cc,coords,1,INTERP_KERNEL::NORM_SEG2);
 688         }
 689       catch(INTERP_KERNEL::Exception& e)
 690         {
 691           std::ostringstream oss; oss << "MEDWriter does not manage polyline cell type because MED file format does not support it ! Maybe it is the source of the problem ? The cause of this exception was " << e.what() << std::endl;
 692           throw INTERP_KERNEL::Exception(oss.str());
 693         }
 694       if((const MEDCouplingUMesh *)subMesh)
 695         {
 696           std::vector<MCAuto<DataArray> > cellFs(AddPartFields2(offset,offset+subMesh->getNumberOfCells(),ds->GetCellData()));
 697           offset+=subMesh->getNumberOfCells();
 698           ms.push_back(MicroField(subMesh,cellFs));
 699         }
 700     }
 701   vtkCellArray *cb(ds->GetPolys());
 702   if(cb)
 703     {
 704       MCAuto<MEDCouplingUMesh> subMesh(BuildMeshFromCellArray(cb,coords,2,INTERP_KERNEL::NORM_POLYGON));
 705       if((const MEDCouplingUMesh *)subMesh)
 706         {
 707           std::vector<MCAuto<DataArray> > cellFs(AddPartFields2(offset,offset+subMesh->getNumberOfCells(),ds->GetCellData()));
 708           offset+=subMesh->getNumberOfCells();
 709           ms.push_back(MicroField(subMesh,cellFs));
 710         }
 711     }
 712   vtkCellArray *ca(ds->GetStrips());
 713   if(ca)
 714     {
 715       MCAuto<DataArrayIdType> ids;
 716       MCAuto<MEDCouplingUMesh> subMesh(BuildMeshFromCellArrayTriangleStrip(ca,coords,ids));
 717       if((const MEDCouplingUMesh *)subMesh)
 718         {
 719           std::vector<MCAuto<DataArray> > cellFs(AddPartFields(ids,ds->GetCellData()));
 720           offset+=subMesh->getNumberOfCells();
 721           ms.push_back(MicroField(subMesh,cellFs));
 722         }
 723     }
 724   AssignSingleGTMeshes(ret,ms,timeStep,tsId);
 725   AddNodeFields(ret,ds->GetPointData(),timeStep,tsId);
 726 }
 727
 728 void ConvertFromUnstructuredGrid(MEDFileData *ret, vtkUnstructuredGrid *ds, const std::vector<int>& context, double timeStep, int tsId)
 729 {
 730   if(!ds || !ret)
 731     throw MZCException("ConvertFromUnstructuredGrid : internal error !");
 732   //
 733   MCAuto<MEDFileMeshes> meshes(MEDFileMeshes::New());
 734   ret->setMeshes(meshes);
 735   MCAuto<MEDFileFields> fields(MEDFileFields::New());
 736   ret->setFields(fields);
 737   //
 738   MCAuto<MEDFileUMesh> umesh(MEDFileUMesh::New());
 739   meshes->pushMesh(umesh);
 740   MCAuto<DataArrayDouble> coords(BuildCoordsFrom(ds));
 741   umesh->setCoords(coords);
 742   umesh->setName(GetMeshNameWithContext(context));
 743   vtkIdType nbCells(ds->GetNumberOfCells());
 744   vtkCellArray *ca(ds->GetCells());
 745   if(!ca)
 746     return ;
 747   vtkIdType nbEnt(ca->GetNumberOfConnectivityEntries());
 748   vtkIdType *caPtr(ca->GetData()->GetPointer(0));
 749   vtkUnsignedCharArray *ct(ds->GetCellTypesArray());
 750   if(!ct)
 751     throw MZCException("ConvertFromUnstructuredGrid : internal error");
 752   vtkIdTypeArray *cla(ds->GetCellLocationsArray());
 753   const vtkIdType *claPtr(cla->GetPointer(0));
 754   if(!cla)
 755     throw MZCException("ConvertFromUnstructuredGrid : internal error 2");
 756   const unsigned char *ctPtr(ct->GetPointer(0));
 757   std::map<int,int> m(ComputeMapOfType());
 758   MCAuto<DataArrayInt> lev(DataArrayInt::New()) ;  lev->alloc(nbCells,1);
 759   int *levPtr(lev->getPointer());
 760   for(vtkIdType i=0;i<nbCells;i++)
 761     {
 762       std::map<int,int>::iterator it(m.find(ctPtr[i]));
 763       if(it!=m.end())
 764         {
 765           const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)(*it).second));
 766           levPtr[i]=cm.getDimension();
 767         }
 768       else
 769         {
 770           if(ctPtr[i]==VTK_POLY_VERTEX)
 771             {
 772               const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_POINT1));
 773               levPtr[i]=cm.getDimension();
 774             }
 775           else
 776             {
 777               std::ostringstream oss; oss << "ConvertFromUnstructuredGrid : at pos #" << i << " unrecognized VTK cell with type =" << ctPtr[i];
 778               throw MZCException(oss.str());
 779             }
 780         }
 781     }
 782   int dummy(0);
 783   MCAuto<DataArrayInt> levs(lev->getDifferentValues());
 784   std::vector< MicroField > ms;
 785   vtkIdTypeArray *faces(ds->GetFaces()),*faceLoc(ds->GetFaceLocations());
 786   for(const int *curLev=levs->begin();curLev!=levs->end();curLev++)
 787     {
 788       MCAuto<MEDCouplingUMesh> m0(MEDCouplingUMesh::New("",*curLev));
 789       m0->setCoords(coords); m0->allocateCells();
 790       MCAuto<DataArrayIdType> cellIdsCurLev(lev->findIdsEqual(*curLev));
 791       for(const mcIdType *cellId=cellIdsCurLev->begin();cellId!=cellIdsCurLev->end();cellId++)
 792         {
 793           int vtkType(ds->GetCellType(*cellId));
 794           std::map<int,int>::iterator it(m.find(vtkType));
 795           INTERP_KERNEL::NormalizedCellType ct=it!=m.end()?(INTERP_KERNEL::NormalizedCellType)((*it).second):INTERP_KERNEL::NORM_POINT1;
 796           if(ct!=INTERP_KERNEL::NORM_POLYHED && vtkType!=VTK_POLY_VERTEX)
 797             {
 798               vtkIdType sz(0);
 799               const vtkIdType *pts(nullptr);
 800               ds->GetCellPoints(*cellId, sz, pts);
 801               std::vector<mcIdType> conn2(pts,pts+sz);
 802               m0->insertNextCell(ct,sz,conn2.data());
 803             }
 804           else if(ct==INTERP_KERNEL::NORM_POLYHED)
 805             {
 806               // # de faces du polyèdre
 807               vtkIdType nbOfFaces(0);
 808               // connectivé des faces (numFace0Pts, id1, id2, id3, numFace1Pts,id1, id2, id3, ...)
 809               const vtkIdType *facPtr(nullptr);
 810               ds->GetFaceStream(*cellId, nbOfFaces, facPtr);
 811               std::vector<mcIdType> conn;
 812               for(vtkIdType k=0;k<nbOfFaces;k++)
 813                 {
 814                   vtkIdType nbOfNodesInFace(*facPtr++);
 815                   std::copy(facPtr,facPtr+nbOfNodesInFace,std::back_inserter(conn));
 816                   if(k<nbOfFaces-1)
 817                     conn.push_back(-1);
 818                   facPtr+=nbOfNodesInFace;
 819                 }
 820               m0->insertNextCell(ct,ToIdType(conn.size()),&conn[0]);
 821             }
 822           else
 823             {
 824               vtkIdType sz(0);
 825               const vtkIdType *pts(nullptr);
 826               ds->GetCellPoints(*cellId, sz, pts);
 827               if(sz!=1)
 828                 throw MZCException("ConvertFromUnstructuredGrid : non single poly vertex not managed by MED !");
 829               m0->insertNextCell(ct,1,(const mcIdType*)pts);
 830             }
 831         }
 832       std::vector<MCAuto<DataArray> > cellFs(AddPartFields(cellIdsCurLev,ds->GetCellData()));
 833       ms.push_back(MicroField(m0,cellFs));
 834     }
 835   AssignSingleGTMeshes(ret,ms,timeStep,tsId);
 836   AddNodeFields(ret,ds->GetPointData(),timeStep,tsId);
 837 }
 838
 839 ///////////////////
 840
 841 void WriteMEDFileFromVTKDataSet(MEDFileData *mfd, vtkDataSet *ds, const std::vector<int>& context, double timeStep, int tsId)
 842 {
 843   if(!ds || !mfd)
 844     throw MZCException("Internal error in WriteMEDFileFromVTKDataSet.");
 845   vtkPolyData *ds2(vtkPolyData::SafeDownCast(ds));
 846   vtkUnstructuredGrid *ds3(vtkUnstructuredGrid::SafeDownCast(ds));
 847   vtkRectilinearGrid *ds4(vtkRectilinearGrid::SafeDownCast(ds));
 848   if(ds2)
 849     {
 850       ConvertFromPolyData(mfd,ds2,context,timeStep,tsId);
 851     }
 852   else if(ds3)
 853     {
 854       ConvertFromUnstructuredGrid(mfd,ds3,context,timeStep,tsId);
 855     }
 856   else if(ds4)
 857     {
 858       ConvertFromRectilinearGrid(mfd,ds4,context,timeStep,tsId);
 859     }
 860   else
 861     throw MZCException("Unrecognized vtkDataSet ! Sorry ! Try to convert it to UnstructuredGrid to be able to write it !");
 862 }
 863
 864 void WriteMEDFileFromVTKMultiBlock(MEDFileData *mfd, vtkMultiBlockDataSet *ds, const std::vector<int>& context, double timeStep, int tsId)
 865 {
 866   if(!ds || !mfd)
 867     throw MZCException("Internal error in WriteMEDFileFromVTKMultiBlock.");
 868   int nbBlocks(ds->GetNumberOfBlocks());
 869   if(nbBlocks==1 && context.empty())
 870     {
 871       vtkDataObject *uniqueElt(ds->GetBlock(0));
 872       if(!uniqueElt)
 873         throw MZCException("Unique elt in multiblock is NULL !");
 874       vtkDataSet *uniqueEltc(vtkDataSet::SafeDownCast(uniqueElt));
 875       if(uniqueEltc)
 876         {
 877           WriteMEDFileFromVTKDataSet(mfd,uniqueEltc,context,timeStep,tsId);
 878           return ;
 879         }
 880     }
 881   for(int i=0;i<nbBlocks;i++)
 882     {
 883       vtkDataObject *elt(ds->GetBlock(i));
 884       std::vector<int> context2;
 885       context2.push_back(i);
 886       if(!elt)
 887         {
 888           std::ostringstream oss; oss << "In context ";
 889           std::copy(context.begin(),context.end(),std::ostream_iterator<int>(oss," "));
 890           oss << " at pos #" << i << " elt is NULL !";
 891           throw MZCException(oss.str());
 892         }
 893       vtkDataSet *elt1(vtkDataSet::SafeDownCast(elt));
 894       if(elt1)
 895         {
 896           WriteMEDFileFromVTKDataSet(mfd,elt1,context,timeStep,tsId);
 897           continue;
 898         }
 899       vtkMultiBlockDataSet *elt2(vtkMultiBlockDataSet::SafeDownCast(elt));
 900       if(elt2)
 901         {
 902           WriteMEDFileFromVTKMultiBlock(mfd,elt2,context,timeStep,tsId);
 903           continue;
 904         }
 905       std::ostringstream oss; oss << "In context ";
 906       std::copy(context.begin(),context.end(),std::ostream_iterator<int>(oss," "));
 907       oss << " at pos #" << i << " elt not recognized data type !";
 908       throw MZCException(oss.str());
 909     }
 910 }
 911
 912 void WriteMEDFileFromVTKGDS(MEDFileData *mfd, vtkDataObject *input, double timeStep, int tsId)
 913 {
 914   if(!input || !mfd)
 915     throw MZCException("WriteMEDFileFromVTKGDS : internal error !");
 916   std::vector<int> context;
 917   vtkDataSet *input1(vtkDataSet::SafeDownCast(input));
 918   if(input1)
 919     {
 920       WriteMEDFileFromVTKDataSet(mfd,input1,context,timeStep,tsId);
 921       return ;
 922     }
 923   vtkMultiBlockDataSet *input2(vtkMultiBlockDataSet::SafeDownCast(input));
 924   if(input2)
 925     {
 926       WriteMEDFileFromVTKMultiBlock(mfd,input2,context,timeStep,tsId);
 927       return ;
 928     }
 929   throw MZCException("WriteMEDFileFromVTKGDS : not recognized data type !");
 930 }
 931
 932 void PutFamGrpInfoIfAny(MEDFileData *mfd, const std::string& meshName, const std::vector<Grp>& groups, const std::vector<Fam>& fams)
 933 {
 934   if(!mfd)
 935     return ;
 936   if(meshName.empty())
 937     return ;
 938   MEDFileMeshes *meshes(mfd->getMeshes());
 939   if(!meshes)
 940     return ;
 941   if(meshes->getNumberOfMeshes()!=1)
 942     return ;
 943   MEDFileMesh *mm(meshes->getMeshAtPos(0));
 944   if(!mm)
 945     return ;
 946   mm->setName(meshName);
 947   for(std::vector<Fam>::const_iterator it=fams.begin();it!=fams.end();it++)
 948     mm->setFamilyId((*it).getName(),(*it).getID());
 949   for(std::vector<Grp>::const_iterator it=groups.begin();it!=groups.end();it++)
 950     mm->setFamiliesOnGroup((*it).getName(),(*it).getFamilies());
 951   MEDFileFields *fields(mfd->getFields());
 952   if(!fields)
 953     return ;
 954   for(int i=0;i<fields->getNumberOfFields();i++)
 955     {
 956       MEDFileAnyTypeFieldMultiTS *fmts(fields->getFieldAtPos(i));
 957       if(!fmts)
 958         continue;
 959       fmts->setMeshName(meshName);
 960     }
 961 }