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