CO2INJECTION: hdf5dealwriter.cpp Source File

00001 #include "hdf5dealwriter.h"
00002 #include <assert.h>
00003 #include <hdf5_hl.h>
00004 #include "globals.h"
00005 #include "dealbase.h"
00006 #include "vecdouble.h"
00007 
00008 HDF5DealWriter* HDF5DealWriter::ptr = NULL;
00009 
00010 HDF5DealWriter::HDF5DealWriter() 
00011 {
00012 }
00013 
00014 
00015 
00016 
00017 
00018 
00019 
00020 
00021 
00022 #define HDF5_BUFF_SIZE 1024
00023 void HDF5DealWriter::writeTriangulation(string triaName,Triangulation<2> &tria) 
00024 {
00025   assert(getFile() != -1);
00026   double buff[HDF5_BUFF_SIZE][2];
00027   int buff_int[HDF5_BUFF_SIZE][4];
00028   
00029   //Memorize a triangulacao.
00030   this->registerTriangulation(triaName,tria);
00031 
00032 
00033   //Criar Grupo
00034   char strPath[500];
00035   sprintf(strPath,"/Triangulations/%s",triaName.c_str());
00036   hid_t triaGroup=H5Gcreate1(getFile(),strPath,0);
00037   
00038   H5LTset_attribute_string(getFile(),strPath,"Type","GRID_2D");
00039 
00040   
00041   hsize_t dim[]={tria.n_vertices(),2};
00042   hid_t dataspace = H5Screate_simple(2,dim,NULL);
00043   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00044   H5Tset_order(datatype,H5T_ORDER_LE);
00045   hid_t dataset = H5Dcreate1(triaGroup, "vertices", datatype, dataspace, H5P_DEFAULT);
00046 
00047   int nSteps = tria.n_vertices()/HDF5_BUFF_SIZE;
00048   hsize_t count[]={HDF5_BUFF_SIZE,2};
00049   hsize_t offset[]={0,0};
00050   unsigned k=0;
00051   for (int i=0;i<=nSteps;i++)
00052   {
00053     int MaxJ = (i==nSteps) ? (tria.n_vertices() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00054     count[0]=MaxJ;
00055     for (int j=0;j<MaxJ;j++)
00056     {
00057       buff[j][0]= tria.get_vertices()[k](0);
00058       buff[j][1]= tria.get_vertices()[k](1);
00059       k++;
00060     }
00061     writeInChunks(dataset,offset,count,buff);
00062     offset[0]+=MaxJ;
00063   }
00064   H5Sclose(dataspace);
00065   H5Tclose(datatype);
00066   H5Dclose(dataset);
00067 
00068   dim[0]=tria.n_active_cells();
00069   dim[1]=4;
00070   dataspace = H5Screate_simple(2,dim,NULL);
00071   datatype  = H5Tcopy(H5T_NATIVE_INT);
00072   H5Tset_order(datatype,H5T_ORDER_LE);
00073   dataset = H5Dcreate1(triaGroup, "connections", datatype, dataspace, H5P_DEFAULT);
00074 
00075 
00076   nSteps = tria.n_active_cells()/HDF5_BUFF_SIZE;
00077   count[1]=4;
00078   offset[0]=0;
00079   Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
00080                                          endc = tria.end();
00081 
00082   for (int i=0;i<=nSteps;i++)
00083   {
00084     int MaxJ = (i==nSteps) ? (tria.n_active_cells() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00085     count[0]=MaxJ;
00086     for (int j=0;j<MaxJ;j++)
00087     {
00088       buff_int[j][0]= cell->vertex_index(BL_VERTEX);
00089       buff_int[j][1]= cell->vertex_index(BR_VERTEX);
00090       buff_int[j][2]= cell->vertex_index(UR_VERTEX);
00091       buff_int[j][3]= cell->vertex_index(UL_VERTEX);
00092       cell++;
00093     }
00094     writeInChunks(dataset,offset,count,buff_int);
00095     offset[0]+=MaxJ;
00096   }
00097   H5Sclose(dataspace);
00098   H5Tclose(datatype);
00099   H5Dclose(dataset);
00100   H5Gclose(triaGroup);
00101   
00102 }
00103 
00104 
00105 
00106 
00107 void HDF5DealWriter::writeTriangulation(string triaName,Triangulation<3> &tria) 
00108 {
00109   assert(getFile() != -1);
00110   double buff[HDF5_BUFF_SIZE][3];
00111   int buff_int[HDF5_BUFF_SIZE][8];
00112   
00113   this->registerTriangulation(triaName,tria);
00114   
00115   //Criar Grupo
00116   char strPath[500];
00117   sprintf(strPath,"/Triangulations/%s",triaName.c_str());
00118   hid_t triaGroup=H5Gcreate1(getFile(),strPath,0);
00119 
00120   H5LTset_attribute_string(getFile(),strPath,"Type","GRID_3D");
00121 
00122   
00123   
00124   hsize_t dim[]={tria.n_vertices(),3};
00125   hid_t dataspace = H5Screate_simple(2,dim,NULL);
00126   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00127   H5Tset_order(datatype,H5T_ORDER_LE);
00128   hid_t dataset = H5Dcreate1(triaGroup, "vertices", datatype, dataspace, H5P_DEFAULT);
00129 
00130   int nSteps = tria.n_vertices()/HDF5_BUFF_SIZE;
00131   hsize_t count[]={HDF5_BUFF_SIZE,3};
00132   hsize_t offset[]={0,0};
00133   unsigned k=0;
00134   for (int i=0;i<=nSteps;i++)
00135   {
00136     int MaxJ = (i==nSteps) ? (tria.n_vertices() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00137     count[0]=MaxJ;
00138     for (int j=0;j<MaxJ;j++)
00139     {
00140       buff[j][0]= tria.get_vertices()[k](0);
00141       buff[j][1]= tria.get_vertices()[k](1);
00142       buff[j][2]= tria.get_vertices()[k](2);
00143       k++;
00144     }
00145     writeInChunks(dataset,offset,count,buff);
00146     offset[0]+=MaxJ;
00147   }
00148   H5Sclose(dataspace);
00149   H5Tclose(datatype);
00150   H5Dclose(dataset);
00151 
00152   dim[0]=tria.n_active_cells();
00153   dim[1]=8;
00154   dataspace = H5Screate_simple(2,dim,NULL);
00155   datatype  = H5Tcopy(H5T_NATIVE_INT);
00156   H5Tset_order(datatype,H5T_ORDER_LE);
00157   dataset = H5Dcreate1(triaGroup, "connections", datatype, dataspace, H5P_DEFAULT);
00158 
00159 
00160   nSteps = tria.n_active_cells()/HDF5_BUFF_SIZE;
00161   count[1]=8;
00162   offset[0]=0;
00163   Triangulation<3>::active_cell_iterator cell = tria.begin_active(),
00164                                          endc = tria.end();
00165 
00166   for (int i=0;i<=nSteps;i++)
00167   {
00168     int MaxJ = (i==nSteps) ? (tria.n_active_cells() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00169     count[0]=MaxJ;
00170     //Write the cell points conectivity. The order of the cell vertex
00171     //is like vtkVOXEL 
00172     for (int j=0;j<MaxJ;j++)
00173     {
00174       buff_int[j][0]= cell->vertex_index(VERTEX_000);
00175       buff_int[j][1]= cell->vertex_index(VERTEX_100);
00176       buff_int[j][2]= cell->vertex_index(VERTEX_010);
00177       buff_int[j][3]= cell->vertex_index(VERTEX_110);
00178       buff_int[j][4]= cell->vertex_index(VERTEX_001);
00179       buff_int[j][5]= cell->vertex_index(VERTEX_101);
00180       buff_int[j][6]= cell->vertex_index(VERTEX_011);
00181       buff_int[j][7]= cell->vertex_index(VERTEX_111);
00182       cell++;
00183     }
00184     writeInChunks(dataset,offset,count,buff_int);
00185     offset[0]+=MaxJ;
00186   }
00187   H5Sclose(dataspace);
00188   H5Tclose(datatype);
00189   H5Dclose(dataset);
00190   H5Gclose(triaGroup);
00191   
00192 }
00193 
00194 
00195 /*
00196 void HDF5DealWriter::writeStructuredGrid(string triaName,int nX,int nY, int nZ, double X, double Y, double Z) 
00197 {
00198   assert(getFile() != -1);
00199   
00200   //Memorize a triangulacao.
00201   struct TriaInformation inf( (nX+1)*(nZ+1)*(nY+1), nX*nY*nZ,STRUCTURED_GRID);
00202   triaHeap[triaName]=inf;
00203 
00204   //Criar Grupo
00205   char strPath[500];
00206   sprintf(strPath,"/Triangulations/%s",triaName.c_str());
00207   hid_t triaGroup=H5Gcreate1(getFile(),strPath,0);
00208 
00209   H5LTset_attribute_string(getFile(),strPath,"Type","STRUCTURED_GRID");
00210 
00211   //Create the dataset
00212   hsize_t dim[]={6};
00213   hid_t dataspace = H5Screate_simple(1,dim,NULL);
00214   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00215   H5Tset_order(datatype,H5T_ORDER_LE);
00216   hid_t dataset = H5Dcreate1(triaGroup, "Geometry_Information", datatype, dataspace, H5P_DEFAULT);
00217 
00218   //Copy the data.
00219   double d[]={nX,nY,nZ,X,Y,Z};
00220   H5Dwrite(dataset,H5T_NATIVE_DOUBLE,H5S_ALL,H5S_ALL,H5P_DEFAULT, d);
00221 
00222   
00223   H5Sclose(dataspace);
00224   H5Tclose(datatype);
00225   H5Dclose(dataset);
00226   H5Gclose(triaGroup);
00227 
00228 }
00229 */
00230 
00231 
00232 
00233 
00234 
00235 
00236 
00237 HDF5DealWriter& HDF5DealWriter::getHDF5DealWriter() 
00238 {
00239    if (ptr == NULL)
00240      ptr = new HDF5DealWriter();
00241    return *ptr;
00242 }
00243 
00254 void HDF5DealWriter::writeScalarField_(const VecDouble &vec, std::string fieldName, std::string triaId) 
00255 {
00256  if (triaId == "")
00257  {
00258    triaId = this->getNameOfLastRegisteredTriangulation();
00259  }
00260   
00261     
00262   validateField(triaId,fieldName);
00263   assert(getFile() != -1);
00264   assert(getTriaInformation(triaId).n_vertices == vec.size() ||
00265          getTriaInformation(triaId).n_cells == vec.size());
00266 
00267   char str[2000];
00268   //Validation of the parameters.
00269   
00270   //Now get the counter of this field and increment it. 
00271   unsigned refCount = getFieldCount(fieldName);
00272   if (refCount == 0)
00273   {
00274      //This field is new, create a directory entry for the field name and initialize the counter
00275     sprintf(str,"/DataSets/%s",fieldName.c_str());
00276     H5Gclose(H5Gcreate1(getFile(),str,0));    
00277   }
00278     
00279     
00280     
00281     
00282   
00283   //Now write the dataset.
00284   //Create the dataset
00285   hsize_t dim[]={vec.size()};
00286   hid_t dataspace = H5Screate_simple(1,dim,NULL);
00287   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00288   H5Tset_order(datatype,H5T_ORDER_LE);
00289   sprintf(str,"/DataSets/%s/%s_%d",fieldName.c_str(),fieldName.c_str(),getFieldCount(fieldName));
00290   hid_t dataset = H5Dcreate1(getFile(), str, datatype, dataspace, H5P_DEFAULT);
00291 
00292   //Copy the data.
00293   H5Dwrite(dataset,H5T_NATIVE_DOUBLE,H5S_ALL,H5S_ALL,H5P_DEFAULT, vec.begin());
00294   
00295   //Set attributes
00296   H5LTset_attribute_string(getFile(),str,"Triangulation",triaId.c_str());
00297   this->writeContextVariables(str);
00298 
00299   //Close the resources.3
00300   H5Sclose(dataspace);
00301   H5Tclose(datatype);
00302   H5Dclose(dataset);
00303 
00304   //Increment the field count
00305   incFieldCount(fieldName);
00306 
00307 }
00308 
00309 
00310 
00311 
00312 
00313 
00314 
00315 
00316 
00321 void HDF5DealWriter::registerTriangulation(string triaId, Triangulation<3> &tria) 
00322 {
00323   struct TriaInformation inf(tria.n_vertices(),tria.n_active_cells());
00324   appendMeshInfo(triaId,inf);
00325 
00326 
00327 }
00328 
00329 
00334 void HDF5DealWriter::registerTriangulation(string triaId, Triangulation<2> &tria) 
00335 {
00336   struct TriaInformation inf(tria.n_vertices(),tria.n_active_cells());
00337   appendMeshInfo(triaId,inf);
00338 
00339 }
00340 
00341 
00342 
00343 
00344 HDF5DealWriter::~HDF5DealWriter() 
00345 {
00346 
00347 }
00348 
00349 
00350 
00351 
00352 
00353 
00363 void HDF5DealWriter::sortVerticesValuesToStructuredGrid(Triangulation<3> &tria, const VecDouble &vOrig,VecDouble &vDest) 
00364 {
00365   assert(tria.n_vertices() == vOrig.size());
00366   assert(vDest.size() == vOrig.size());
00367   
00368 
00369   Triangulation<3>::active_cell_iterator cellX,cellY,cellZ = tria.begin_active();
00370   //Walk first in X direction
00371   unsigned i=0;
00372   while (1)
00373   {
00374     cellY = cellZ;
00375     while(1)
00376     {
00377       //Run in X
00378       cellX = cellY;
00379       while(1)
00380       {
00381         vDest(i++)=vOrig(cellX->vertex_index(VERTEX_000));
00382         if (cellX->neighbor_index(RIGHT_CELL) == DealBase::INVALID_INDEX)
00383         {
00384           vDest(i++)=vOrig(cellX->vertex_index(VERTEX_100));
00385           break;
00386         }
00387         else
00388           cellX = cellX->neighbor(RIGHT_CELL);
00389       }
00390       if (cellY->neighbor_index(BACK_CELL) == DealBase::INVALID_INDEX)
00391       {
00392         cellX = cellY;
00393         while(1)
00394         {
00395           vDest(i++)=vOrig(cellX->vertex_index(VERTEX_010));
00396           if (cellX->neighbor_index(RIGHT_CELL) == DealBase::INVALID_INDEX)
00397           {
00398             vDest(i++)=vOrig(cellX->vertex_index(VERTEX_110));
00399             break;
00400           }
00401           else
00402             cellX = cellX->neighbor(RIGHT_CELL);
00403         }
00404         break;
00405       }
00406       else
00407         cellY=cellY->neighbor(BACK_CELL);
00408     }
00409     if (cellZ->neighbor_index(UP_CELL) == DealBase::INVALID_INDEX)
00410     {
00411       cellY = cellZ;
00412       while(1)
00413       {
00414         //Run in X
00415         cellX = cellY;
00416         while(1)
00417         {
00418           vDest(i++)=vOrig(cellX->vertex_index(VERTEX_001));
00419           if (cellX->neighbor_index(RIGHT_CELL) == DealBase::INVALID_INDEX)
00420           {
00421             vDest(i++)=vOrig(cellX->vertex_index(VERTEX_101));
00422             break;
00423           }
00424           else
00425             cellX = cellX->neighbor(RIGHT_CELL);
00426         }
00427         if (cellY->neighbor_index(BACK_CELL) == DealBase::INVALID_INDEX)
00428         {
00429           cellX = cellY;
00430           while(1)
00431           {
00432             vDest(i++)=vOrig(cellX->vertex_index(VERTEX_011));
00433             if (cellX->neighbor_index(RIGHT_CELL) == DealBase::INVALID_INDEX)
00434             {
00435               vDest(i++)=vOrig(cellX->vertex_index(VERTEX_111));
00436               break;
00437             }
00438             else
00439               cellX = cellX->neighbor(RIGHT_CELL);
00440           }
00441           break;
00442         }
00443         else
00444           cellY=cellY->neighbor(BACK_CELL);
00445       }
00446       break;
00447     }
00448     else
00449     {
00450       cellZ = cellZ->neighbor(UP_CELL);
00451     }
00452     
00453   }
00454 }
00455