src/Plugins/MEDWriter/IO/VTKToMEDMem.cxx

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