src/VISU_I/VISU_Convertor_impl.cxx

   1 using namespace std;
   2
   3 // File:        VISU_Convertor_impl.cxx
   4 // Created:     Fri Jan 10 11:44:54 2003
   5 // Author:      Alexey PETROV
   6 //              <apo@ivanox.nnov.matra-dtv.fr>
   7
   8
   9 #include "VISU_Convertor_impl.hxx"
  10
  11 #include <vtkUnstructuredGridReader.h>
  12 #include <qstring.h>
  13 #include <qfileinfo.h>
  14 #include <vtkCellType.h>
  15 #include <valarray>
  16 #include <memory>
  17 using namespace std;
  18
  19 #ifdef DEBUG
  20 static int MYDEBUG = 1;
  21 static int MYDEBUGWITHFILES = 0;
  22 #else
  23 static int MYDEBUG = 0;
  24 static int MYDEBUGWITHFILES = 0;
  25 #endif
  26 static int PRECISION = 7;
  27
  28 #define MED2VTK(MEDTYPE,VTKTYPE,VTKNBNODES) \
  29  {MEDTYPE,#MEDTYPE,getNbMedNodes(MEDTYPE),VTKTYPE,#VTKTYPE,VTKNBNODES}
  30 Med2vtk med2vtk[MED_NBR_GEOMETRIE_MAILLE] = {
  31   MED2VTK(MED_POINT1,VTK_VERTEX,1),
  32   MED2VTK(MED_SEG2,VTK_LINE,2),
  33   MED2VTK(MED_SEG3,VTK_LINE,2),
  34   MED2VTK(MED_TRIA3,VTK_TRIANGLE,3),
  35   MED2VTK(MED_TRIA6,VTK_TRIANGLE,3),
  36   MED2VTK(MED_QUAD4,VTK_QUAD,4),
  37   MED2VTK(MED_QUAD8,VTK_QUAD,4),
  38   MED2VTK(MED_TETRA4,VTK_TETRA,4),
  39   MED2VTK(MED_TETRA10,VTK_TETRA,4),
  40   MED2VTK(MED_HEXA8,VTK_HEXAHEDRON,8),
  41   MED2VTK(MED_HEXA20,VTK_HEXAHEDRON,8),
  42   MED2VTK(MED_PENTA6,VTK_WEDGE,6),
  43   MED2VTK(MED_PENTA15,VTK_WEDGE,6),
  44   MED2VTK(MED_PYRA5,VTK_PYRAMID,5),
  45   MED2VTK(MED_PYRA13,VTK_PYRAMID,5)
  46 };
  47 #undef MED2VTK
  48
  49 extern "C" {
  50   VISU_Convertor* CreateConvertor(const string& theFileName) throw(std::runtime_error&){
  51     if(QFileInfo(theFileName.c_str()).extension(false) == "med")
  52       return CreateMedConvertor(theFileName);
  53     else
  54       return CreateDatConvertor(theFileName);
  55   }
  56
  57   int getNbMedConnect(int theMedType, int theMedEntity, int theMeshDim){
  58     int anElemDim = theMedType / 100, nsup = 0;
  59     if(theMedEntity == VISU::CELL_ENTITY && anElemDim < theMeshDim) nsup = 1;
  60     return nsup + theMedType % 100;
  61   }
  62
  63   int getNbMedNodes(int geom){
  64     return geom % 100;
  65   }
  66
  67   int getIdMedType(int medType){
  68     for(int i = 0; i < MED_NBR_GEOMETRIE_MAILLE; i++)
  69       if(med2vtk[i].medType == medType) return i;
  70     return -1;
  71   }
  72 }
  73
  74 int VISU_Convertor_impl::ToFile(const string& theFileName) throw(std::runtime_error&){
  75   if(QFileInfo(theFileName.c_str()).extension(false) == "med")
  76     return ToMedFile(theFileName);
  77   else
  78     return ToDatFile(theFileName);
  79 }
  80
  81 VISU_Convertor_impl::VISU_Convertor_impl() : myIsDone(false) {}
  82
  83 VISU_Convertor_impl::~VISU_Convertor_impl() {}
  84
  85 VISU_Convertor::OutputType* VISU_Convertor_impl::GetMeshOnEntity(const string& theMeshName,
  86                                                                  const VISU::TEntity& theEntity,
  87                                                                  const string& theFamilyName)
  88   throw (std::runtime_error&)
  89 {
  90   if(MYDEBUG)
  91     MESSAGE("GetMeshOnEntity - theMeshName = '"<<theMeshName<<
  92             "'; theEntity = "<<theEntity<<"; theFamilyName = '"<<theFamilyName<<"'");
  93   //Cheching possibility do the query
  94   if(myMeshMap.find(theMeshName) == myMeshMap.end())
  95     throw std::runtime_error("MeshOnEntityToString >> There is no mesh with the name!!!");
  96   VISU::TMesh& aMesh = myMeshMap[theMeshName];
  97   VISU::TMeshOnEntityMap& aMeshOnEntityMap = aMesh.myMeshOnEntityMap;
  98   if(aMeshOnEntityMap.find(theEntity) == aMeshOnEntityMap.end())
  99     throw std::runtime_error("MeshOnEntityToString >> There is no mesh on the entity!!!");
 100   VISU::TMeshOnEntity& aMeshOnEntity = aMeshOnEntityMap[theEntity];
 101   VISU::TFamily* pFamily = VISU::GetFamily(aMeshOnEntity,theFamilyName);
 102   VISU::TVTKReader* pReader;
 103   if(pFamily != NULL)
 104     pReader = &(pFamily->myStorage);
 105   else
 106     pReader = &(aMeshOnEntity.myStorage);
 107   VISU::TVTKReader& aReader = *pReader;
 108   //Main part of code
 109   if(aReader.get() == NULL){
 110     LoadMeshOnEntity(aMeshOnEntity,theFamilyName);
 111     ostrstream strOut;
 112     strOut<<GetHead(theMeshName + dtos("-%d-",theEntity) + theFamilyName);
 113     strOut<<GetPoints(aMesh);
 114     pair<int,int> aCellsDim = aMeshOnEntity.GetCellsDims(theFamilyName);
 115     int aNbCells = aCellsDim.first, aCellsSize = aCellsDim.second;
 116     ostrstream strCellsOut, strTypesOut;
 117     GetCells(strCellsOut,strTypesOut,aMeshOnEntity,theFamilyName);
 118     strCellsOut<<ends;
 119     strTypesOut<<ends;
 120     strOut<<"CELLS "<<aNbCells<<"\t"<<aCellsSize<<endl;
 121     strOut<<strCellsOut.str()<<endl;
 122     strOut<<"CELL_TYPES "<<aNbCells<<endl;
 123     strOut<<strTypesOut.str()<<endl;
 124     strOut<<ends;
 125     aReader.reset(OutputType::New());
 126     //aReader = OutputType::New();
 127     //aReader->DebugOn();
 128     aReader->ReadFromInputStringOn();
 129     aReader->SetInputString(strOut.str());
 130     //aReader->Update();
 131     //aReader->Print(cout);
 132     if(MYDEBUGWITHFILES){
 133       string aMeshName = QString(theMeshName.c_str()).simplifyWhiteSpace().latin1();
 134       string aFamilyName = QString(theFamilyName.c_str()).simplifyWhiteSpace().latin1();
 135       string aFileName = string("/users/")+getenv("USER")+"/"+getenv("USER")+"-";
 136       aFileName += aMeshName + dtos("-%d-",theEntity) + aFamilyName + "-Conv.vtk";
 137       ofstream stmOut(aFileName.c_str(),ios::out);
 138       stmOut<<strOut.str();
 139     }
 140   }
 141   return aReader.get();
 142   //return aReader;
 143 }
 144
 145 VISU_Convertor::OutputType* VISU_Convertor_impl::GetMeshOnGroup(const string& theMeshName,
 146                                                                 const string& theGroupName)
 147      throw(std::runtime_error&)
 148 {
 149   if(MYDEBUG) MESSAGE("GetMeshOnGroup - theMeshName = '"<<theMeshName<<
 150                       "'; theGroupName = '"<<theGroupName<<"'");
 151   //Cheching possibility do the query
 152   if(myMeshMap.find(theMeshName) == myMeshMap.end())
 153     throw std::runtime_error("GetMeshOnGroup >> There is no mesh with the name!!!");
 154   VISU::TMesh& aMesh = myMeshMap[theMeshName];
 155   VISU::TGroupMap& aGroupMap = aMesh.myGroupMap;
 156   VISU::TGroupMap::iterator aGroupMapIter = aGroupMap.find(theGroupName);
 157   if(aGroupMapIter == aGroupMap.end())
 158     throw std::runtime_error("GetMeshOnGroup >> There is no the group in the mesh!!!");
 159   VISU::TGroup& aGroup = aGroupMapIter->second;
 160   const VISU::TFamilyAndEntitySet& aFamilyAndEntitySet = aGroup.myFamilyAndEntitySet;
 161   VISU::TVTKReader& aReader = aGroup.myStorage;
 162   if(aReader.get() == NULL){
 163     LoadMeshOnGroup(aMesh,aFamilyAndEntitySet);
 164     ostrstream strOut;
 165     strOut<<GetHead(theMeshName + "-" + theGroupName);
 166     strOut<<GetPoints(aMesh);
 167     VISU::TFamilyAndEntitySet::const_iterator aFamilyAndEntitySetIter = aFamilyAndEntitySet.begin();
 168     int aNbCells = 0, aCellsSize = 0;
 169     ostrstream strCellsOut, strTypesOut;
 170     for(; aFamilyAndEntitySetIter != aFamilyAndEntitySet.end(); aFamilyAndEntitySetIter++){
 171       const VISU::TFamilyAndEntity& aFamilyAndEntity = *aFamilyAndEntitySetIter;
 172       const string& aFamilyName = aFamilyAndEntity.first;
 173       const VISU::TEntity& anEntity = aFamilyAndEntity.second;
 174       VISU::TMeshOnEntity& aMeshOnEntity = aMesh.myMeshOnEntityMap[anEntity];
 175       pair<int,int> aCellsDim = aMeshOnEntity.GetCellsDims(aFamilyName);
 176       aNbCells += aCellsDim.first;
 177       aCellsSize += aCellsDim.second;
 178       GetCells(strCellsOut,strTypesOut,aMeshOnEntity,aFamilyName);
 179     }
 180     strCellsOut<<ends;
 181     strTypesOut<<ends;
 182     strOut<<"CELLS "<<aNbCells<<"\t"<<aCellsSize<<endl;
 183     strOut<<strCellsOut.str()<<endl;
 184     strOut<<"CELL_TYPES "<<aNbCells<<endl;
 185     strOut<<strTypesOut.str()<<endl;
 186     strOut<<ends;
 187     aReader.reset(OutputType::New());
 188     //aReader = OutputType::New();
 189     aReader->ReadFromInputStringOn();
 190     aReader->SetInputString(strOut.str());
 191     if(MYDEBUGWITHFILES){
 192       string aMeshName = QString(theMeshName.c_str()).simplifyWhiteSpace().latin1();
 193       string aGroupName = QString(theGroupName.c_str()).simplifyWhiteSpace().latin1();
 194       string aFileName = string("/users/")+getenv("USER")+"/"+getenv("USER")+"-";
 195       aFileName += aMeshName + "-" + aGroupName + "-Conv.vtk";
 196       ofstream stmOut(aFileName.c_str(),ios::out);
 197       stmOut<<strOut.str();
 198     }
 199   }
 200   return aReader.get();
 201   //return aReader;
 202 }
 203
 204 VISU_Convertor::OutputType* VISU_Convertor_impl::GetFieldOnMesh(const string& theMeshName,
 205                                                                 const VISU::TEntity& theEntity,
 206                                                                 const string& theFieldName,
 207                                                                 int theStampsNum)
 208      throw(std::runtime_error&)
 209 {
 210   if(MYDEBUG){
 211     MESSAGE("GetFieldOnMesh - theMeshName = '"<<theMeshName<<"; theEntity = "<<theEntity);
 212     MESSAGE("GetFieldOnMesh - theFieldName = '"<<theFieldName<<"'; theStampsNum = "<<theStampsNum);
 213   }
 214   //Cheching possibility do the query
 215   if(myMeshMap.find(theMeshName) == myMeshMap.end())
 216     throw std::runtime_error("GetFieldOnMesh >> There is no mesh with the name!!!");
 217   VISU::TMesh& aMesh = myMeshMap[theMeshName];
 218   VISU::TMeshOnEntityMap& aMeshOnEntityMap = aMesh.myMeshOnEntityMap;
 219   if(aMeshOnEntityMap.find(theEntity) == aMeshOnEntityMap.end())
 220     throw std::runtime_error("GetFieldOnMesh >> There is no mesh on the entity!!!");
 221   VISU::TMeshOnEntity& aMeshOnEntity = aMeshOnEntityMap[theEntity];
 222   VISU::TMeshOnEntity* pVtkMeshOnEntity = &aMeshOnEntity;
 223   if(theEntity == VISU::NODE_ENTITY){
 224     if(aMeshOnEntityMap.find(VISU::CELL_ENTITY) == aMeshOnEntityMap.end())
 225       throw std::runtime_error("GetFieldOnMesh >> There is no mesh on CELL_ENTITY!!!");
 226     pVtkMeshOnEntity = &aMeshOnEntityMap[VISU::CELL_ENTITY];
 227   }
 228   VISU::TMeshOnEntity& aVtkMeshOnEntity = *pVtkMeshOnEntity;
 229   VISU::TFieldMap& aFieldMap = aMeshOnEntity.myFieldMap;
 230   if(aFieldMap.find(theFieldName) == aFieldMap.end())
 231     throw std::runtime_error("GetFieldOnMesh >> There is no field on the mesh!!!");
 232   VISU::TField& aField = aFieldMap[theFieldName];
 233   VISU::TField::TValField& aValField = aField.myValField;
 234   if(aValField.find(theStampsNum) == aValField.end())
 235     throw std::runtime_error("GetFieldOnMesh >> There is no field with the timestamp!!!");
 236   VISU::TField::TValForTime& aValForTime = aValField[theStampsNum];
 237   VISU::TVTKReader& aReader = aValForTime.myStorage;
 238   //Main part of code
 239   if(aReader.get() == NULL){
 240     LoadMeshOnEntity(aVtkMeshOnEntity);
 241     LoadFieldOnMesh(aMesh,aMeshOnEntity,aField,aValForTime);
 242     ostrstream strOut;
 243     strOut<<GetHead(theMeshName + dtos("-%d-",theEntity) + theFieldName + dtos("-%d",theStampsNum));
 244     strOut<<GetPoints(aMesh);
 245     pair<int,int> aCellsDim = aVtkMeshOnEntity.GetCellsDims();
 246     int aNbCells = aCellsDim.first, aCellsSize = aCellsDim.second;
 247     ostrstream strCellsOut, strTypesOut;
 248     GetCells(strCellsOut,strTypesOut,aVtkMeshOnEntity);
 249     strCellsOut<<ends;
 250     strTypesOut<<ends;
 251     strOut<<"CELLS "<<aNbCells<<"\t"<<aCellsSize<<endl;
 252     strOut<<strCellsOut.str()<<endl;
 253     strOut<<"CELL_TYPES "<<aNbCells<<endl;
 254     strOut<<strTypesOut.str()<<endl;
 255     int aNbPoints = aMesh.myPointsCoord.size()/aMesh.myDim;
 256     strOut<<GetField(aField,aValForTime,aMesh.myDim,aNbPoints,aNbCells);
 257     aReader.reset(OutputType::New());
 258     //aReader = OutputType::New();
 259     aReader->ReadFromInputStringOn();
 260     aReader->SetInputString(strOut.str());
 261     if(MYDEBUGWITHFILES){
 262       string aMeshName = QString(theMeshName.c_str()).simplifyWhiteSpace().latin1();
 263       string aFieldName = QString(theFieldName.c_str()).simplifyWhiteSpace().latin1();
 264       string aFileName = string("/users/")+getenv("USER")+"/"+getenv("USER")+"-";
 265       aFileName += aMeshName + dtos("-%d-",theEntity) + aFieldName + dtos("-%d",theStampsNum) + "-Conv.vtk";
 266       ofstream stmOut(aFileName.c_str(),ios::out);
 267       stmOut<<strOut.str();
 268     }
 269   }
 270   return aReader.get();
 271   //return aReader;
 272 }
 273
 274 inline void PrintCells(ostrstream& strCellsOut, const VISU::TMeshOnEntity::TConnect& theVector,
 275                        ostrstream& strTypesOut, const int theVtkType)
 276 {
 277   strTypesOut<<theVtkType<<"\n";
 278   int kEnd = theVector.size();
 279   strCellsOut<<kEnd<<"\t";
 280   for(int k = 0; k < kEnd; k++)
 281     strCellsOut<<theVector[k]<<"\t";
 282   strCellsOut<<endl;
 283 }
 284
 285 int VISU_Convertor_impl::GetCells(ostrstream& strCellsOut, ostrstream& strTypesOut,
 286                                   const VISU::TMeshOnEntity& theMeshOnEntity,
 287                                   const string& theFamilyName) const
 288        throw (std::runtime_error&)
 289 {
 290   //Check on existing family
 291   const VISU::TFamily* pFamily = VISU::GetFamily(theMeshOnEntity,theFamilyName);
 292   bool isFamilyPresent = (pFamily != NULL);
 293   const VISU::TFamily& aFamily = *pFamily;
 294   //Main part of code
 295   if(MYDEBUG) MESSAGE("GetCells - isFamilyPresent = "<<isFamilyPresent);
 296   const VISU::TMeshOnEntity::TCellsConn &aCellsConn = theMeshOnEntity.myCellsConn;
 297   VISU::TMeshOnEntity::TCellsConn::const_iterator aCellsConnIter = aCellsConn.begin();
 298   for(; aCellsConnIter != aCellsConn.end(); aCellsConnIter++){
 299     const VISU::TMeshOnEntity::TConnForCellType& anArray = aCellsConnIter->second;
 300     if(MYDEBUG) MESSAGE("GetCells - anArray.size() = "<<anArray.size());
 301     int aVtkType = aCellsConnIter->first;
 302     if(!isFamilyPresent)
 303       for(int j = 0, jEnd = anArray.size(); j < jEnd; j++)
 304         PrintCells(strCellsOut,anArray[j],strTypesOut,aVtkType);
 305     else{
 306       const VISU::TFamily::TSubMesh& aSubMesh = aFamily.mySubMesh;
 307       VISU::TFamily::TSubMesh::const_iterator aSubMeshIter = aSubMesh.find(aVtkType);
 308       const VISU::TFamily::TSubMeshOnCellType& aSubMeshOnCellType = aSubMeshIter->second;
 309       VISU::TFamily::TSubMeshOnCellType::const_iterator aSubMeshOnCellTypeIter = aSubMeshOnCellType.begin();
 310       for(; aSubMeshOnCellTypeIter != aSubMeshOnCellType.end(); aSubMeshOnCellTypeIter++)
 311         PrintCells(strCellsOut,anArray[*aSubMeshOnCellTypeIter],strTypesOut,aVtkType);
 312     }
 313   }
 314   return 1;
 315 }
 316
 317 string VISU_Convertor_impl::GetHead(const string& theMeshName) const throw (std::runtime_error&){
 318   ostrstream strOut;
 319   strOut<<"# vtk DataFile Version 3.2\n";
 320   strOut<<theMeshName<<endl;
 321   strOut<<"ASCII\n";
 322   strOut<<"DATASET UNSTRUCTURED_GRID\n";
 323   strOut<<ends;
 324   return strOut.str();
 325 }
 326
 327 string VISU_Convertor_impl::GetPoints(const VISU::TMesh& theMesh) const
 328        throw (std::runtime_error&)
 329 {
 330   ostrstream strOut;
 331   strOut.setf(ios::scientific,ios::floatfield);  strOut.precision(PRECISION);
 332   const VISU::TMesh::TPointsCoord& anArray = theMesh.myPointsCoord;
 333   int iEnd = theMesh.myPointsCoord.size();
 334   int aNbPoints = iEnd/theMesh.myDim;
 335   strOut<<"POINTS "<<aNbPoints<<" float\n";
 336   if(MYDEBUG)
 337     MESSAGE("GetPoints - aNbPoints = "<<aNbPoints<<"; myDim = "<<theMesh.myDim);
 338   switch(theMesh.myDim) {
 339   case 1:
 340     for (int i = 0; i < iEnd;)
 341       strOut<<anArray[i++]<<"\n";
 342     break;
 343   case 2:
 344     for (int i = 0; i < iEnd;){
 345       strOut<<anArray[i++]<<"\t";
 346       strOut<<anArray[i++]<<"\t";
 347       strOut<<0.0<<"\n";
 348     }
 349     break;
 350   case 3:
 351     for (int i = 0; i < iEnd;){
 352       strOut<<anArray[i++]<<"\t";
 353       strOut<<anArray[i++]<<"\t";
 354       strOut<<anArray[i++]<<"\n";
 355     }
 356     break;
 357   }
 358   strOut<<ends;
 359   return strOut.str();
 360 }
 361
 362 string VISU_Convertor_impl::GetField(const VISU::TField& theField,
 363                                      const VISU::TField::TValForTime& theValForTime,
 364                                      int theDim, int theNbPoints, int theNbCells) const
 365        throw (std::runtime_error&)
 366 {
 367   ostrstream strOut;
 368   strOut.setf(ios::scientific,ios::floatfield);  strOut.precision(PRECISION);
 369   string aName = GenerateName(theField.myName,theValForTime.myId);
 370   if(MYDEBUG) MESSAGE("GetField - aTime = "<<theValForTime.myId<<"; aName = "<<aName);
 371   switch(theField.myEntity) {
 372   case VISU::NODE_ENTITY :
 373     strOut<<"POINT_DATA "<<theNbPoints<<"\n";
 374     break;
 375   default:
 376     strOut<<"CELL_DATA "<<theNbCells<<"\n";
 377   }
 378   switch(theField.myNbComp) {
 379   case 1:
 380     strOut<<"SCALARS "<<aName<<" float\n";
 381     strOut<<"LOOKUP_TABLE default\n";
 382     break;
 383   default:
 384     strOut<<"VECTORS "<<aName<<" float\n";
 385     break;
 386   }
 387   const VISU::TField::TValForCells& aValForCells = theValForTime.myValForCells;
 388   VISU::TField::TValForCells::const_iterator aValForCellsIter = aValForCells.begin();
 389   for(; aValForCellsIter != aValForCells.end(); aValForCellsIter++) {
 390     const VISU::TField::TValForCellsWithType& anArray = aValForCellsIter->second;
 391     int iEnd = anArray.size()/theField.myNbComp;
 392     int aVtkType = aValForCellsIter->first;
 393     if(MYDEBUG) MESSAGE("GetField -  iEnd = "<<iEnd<<"; aVtkType = "<<aVtkType<<"; theDim = "<<theDim);
 394     if(theField.myNbComp == 1){
 395       for (int i = 0; i < iEnd;) strOut<<anArray[i++]<<"\n";
 396     }else if(theDim == 2){
 397       for (int i = 0, ji = 0; i < iEnd; ji = (++i)*theField.myNbComp){
 398         strOut<<anArray[ji++]<<"\t";
 399         strOut<<anArray[ji++]<<"\t";
 400         strOut<<0.0<<"\n";
 401       }
 402     }else if(theDim == 3){
 403       for (int i = 0, ji = 0; i < iEnd; ji = (++i)*theField.myNbComp){
 404         strOut<<anArray[ji++]<<"\t";
 405         strOut<<anArray[ji++]<<"\t";
 406         strOut<<anArray[ji++]<<"\n";
 407       }
 408     }else
 409       throw std::runtime_error("GetField - There is algorithm for representation the field !!!");
 410   }
 411   strOut<<ends;
 412   return strOut.str();
 413 }