src/Plugins/MEDWriter/IO/VTKToMEDMem.cxx

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