src/RateOfFlowThroughSection/plugin/RateOfFlowThroughSectionModule/VTKToMEDMem.cxx

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