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