CO2INJECTION: hdf5orthowriter.cpp Source File

00001 #include "hdf5orthowriter.h"
00002 #include <hdf5_hl.h>
00003 #include <hdf5.h>
00004 #include "numericmethods.h"
00005 #include "xdmfwriter.h"
00006 HDF5OrthoWriter* HDF5OrthoWriter::ptr = NULL;
00007 
00008 
00009  HDF5OrthoWriter::HDF5OrthoWriter() 
00010 {
00011   this->_supress_output=false;
00012 }
00013 
00014 
00015 
00016  HDF5OrthoWriter::~HDF5OrthoWriter() 
00017 {
00018   std::string file=getFileName();
00019   if (!file.empty())
00020   {
00021     XdmfWriter xdmf;
00022     xdmf.writeXdmf(file);
00023   }
00024   
00025 }
00026 
00027 
00028 
00029 
00030 
00031 
00032 
00033 HDF5OrthoWriter& HDF5OrthoWriter::getHDF5OrthoWriter() 
00034 {
00035   if (ptr == NULL)
00036     ptr=new HDF5OrthoWriter();
00037   return *ptr;
00038     
00039 }
00040 
00041 
00042 
00043 #define HDF5_BUFF_SIZE 1024
00044 void HDF5OrthoWriter::writeMesh(string meshName,OrthoMesh &mesh) 
00045 {
00046   if (_supress_output)
00047     return;
00048 
00049   if (_defaultMesh.empty())
00050     _defaultMesh=meshName;
00051   assert(getFile() != -1);
00052   double buff[HDF5_BUFF_SIZE][3];
00053   int buff_int[HDF5_BUFF_SIZE][8];
00054   
00055   putMeshIntoDictionary(meshName,mesh.numVertices(),mesh.numCells());
00056   
00057   //Criar Grupo
00058   char strPath[500];
00059   sprintf(strPath,"/Triangulations/%s",meshName.c_str());
00060   hid_t triaGroup=H5Gcreate1(getFile(),strPath,0);
00061 
00062   H5LTset_attribute_string(getFile(),strPath,"Type","GRID_3D");
00063   HDF5OrthoInfo info(mesh.getP(),mesh.getQ(),mesh.numElemX(),mesh.numElemY(),mesh.numElemZ(),mesh.getWells());
00064   info.writeOrthoMeshInfo(triaGroup);
00065 
00066   
00067   
00068   hsize_t dim[]={mesh.numVertices(),3};
00069   hid_t dataspace = H5Screate_simple(2,dim,NULL);
00070   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00071   H5Tset_order(datatype,H5T_ORDER_LE);
00072   hid_t dataset = H5Dcreate1(triaGroup, "vertices", datatype, dataspace, H5P_DEFAULT);
00073 
00074   int nSteps = mesh.numVertices()/HDF5_BUFF_SIZE;
00075   hsize_t count[]={HDF5_BUFF_SIZE,3};
00076   hsize_t offset[]={0,0};
00077   OrthoMesh::Vertex_It vIt = mesh.begin_vertice();
00078   for (int i=0;i<=nSteps;i++)
00079   {
00080     int MaxJ = (i==nSteps) ? (mesh.numVertices() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00081     count[0]=MaxJ;
00082     for (int j=0;j<MaxJ;j++)
00083     {
00084       Point3D p =vIt->getPoint();
00085       buff[j][0]= p[0];
00086       buff[j][1]= p[1];
00087       buff[j][2]= p[2];
00088       vIt++;
00089     }
00090     writeInChunks(dataset,offset,count,buff);
00091     offset[0]+=MaxJ;
00092   }
00093   H5Sclose(dataspace);
00094   H5Tclose(datatype);
00095   H5Dclose(dataset);
00096 
00097   dim[0]=mesh.numCells();
00098   dim[1]=8;
00099   dataspace = H5Screate_simple(2,dim,NULL);
00100   datatype  = H5Tcopy(H5T_NATIVE_INT);
00101   H5Tset_order(datatype,H5T_ORDER_LE);
00102   dataset = H5Dcreate1(triaGroup, "connections", datatype, dataspace, H5P_DEFAULT);
00103 
00104 
00105   nSteps = mesh.numCells()/HDF5_BUFF_SIZE;
00106   count[1]=8;
00107   offset[0]=0;
00108   OrthoMesh::Cell_It cell = mesh.begin_cell(),
00109     endc = mesh.end_cell();
00110 
00111   for (int i=0;i<=nSteps;i++)
00112   {
00113     int MaxJ = (i==nSteps) ? (mesh.numCells() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00114     count[0]=MaxJ;
00115     //Write the cell points conectivity. The order of the cell vertex
00116     //is like vtkVOXEL 
00117     for (int j=0;j<MaxJ;j++)
00118     {
00119       buff_int[j][0]= cell->vertex_index(VERTEX_000);
00120       buff_int[j][1]= cell->vertex_index(VERTEX_100);
00121       buff_int[j][2]= cell->vertex_index(VERTEX_010);
00122       buff_int[j][3]= cell->vertex_index(VERTEX_110);
00123       buff_int[j][4]= cell->vertex_index(VERTEX_001);
00124       buff_int[j][5]= cell->vertex_index(VERTEX_101);
00125       buff_int[j][6]= cell->vertex_index(VERTEX_011);
00126       buff_int[j][7]= cell->vertex_index(VERTEX_111);
00127       cell++;
00128     }
00129     writeInChunks(dataset,offset,count,buff_int);
00130     offset[0]+=MaxJ;
00131   }
00132   H5Sclose(dataspace);
00133   H5Tclose(datatype);
00134   H5Dclose(dataset);
00135   H5Gclose(triaGroup);
00136 
00137 }
00138 
00139 
00140 void HDF5OrthoWriter::writeScalarField(const VecDouble &vec, std::string sFieldName, std::string triaId)
00141 {
00142   if (_supress_output)
00143     return;
00144 
00145   //Set the dataset name
00146   char datasetName[1000];
00147   char datasetGroup[1000];
00148   char attValue[10];
00149 
00150   //If the user did not give the name of the mesh
00151   //use the last registered mesh
00152  if (triaId == "")
00153  {
00154    triaId = _defaultMesh;
00155  }
00156 
00157   //Validation of the parameters
00158  validateField(triaId,sFieldName);
00159 
00160   //Now get the counter of this field and increment it. 
00161   unsigned refCount = getFieldCount(sFieldName);
00162   if (refCount == 0)
00163   {
00164      //This field is new, create a directory entry for the field name and initialize the counter
00165     sprintf(datasetName,"/DataSets/%s",sFieldName.c_str());
00166     H5Gclose(H5Gcreate1(getFile(),datasetName,0));    
00167   }
00168   double min,max;
00169   NumericMethods::vectorMinMaxValues(vec,min,max);
00170   getFieldInfo(sFieldName).updateMinMaxValues(min,max);
00171   
00172     
00173   sprintf(datasetName,"/DataSets/%s/%s_%d",sFieldName.c_str(),sFieldName.c_str(),getFieldCount(sFieldName));
00174 
00175 
00176 
00177   //write it
00178   writeDataSet(vec,datasetName,triaId);
00179 
00180 
00181   sprintf(datasetGroup,"/DataSets/%s",sFieldName.c_str());
00182   sprintf(attValue,"%d",refCount+1);
00183   setAtt(datasetGroup,"NFields",attValue); 
00184   setAtt(datasetGroup,"Max",getFieldInfo(sFieldName).max);
00185   setAtt(datasetGroup,"Min",getFieldInfo(sFieldName).min);
00186   
00187   //Increment the field count
00188   incFieldCount(sFieldName);
00189 }
00190 
00191 herr_t exitA(hid_t t,void *b)
00192 {
00193   H5Eprint(t,stdout);
00194   abort();
00195 }
00196 
00197 void HDF5OrthoWriter::writeDataSet(const VecDouble &vec, std::string datasetName, std::string triaId) 
00198 {
00199   if (_supress_output)
00200     return;
00201   H5Eset_auto2(H5E_DEFAULT,(herr_t (*)(hid_t,void*))   &exitA,stdout);
00202 
00203   //If the user did not give the name of the mesh
00204  //use the last registered mesh
00205  if (triaId == "")
00206  {
00207    triaId = this->getNameOfLastRegisteredTriangulation();
00208  }
00209 
00210  if (!(getTriaInformation(triaId).n_vertices == vec.size() ||
00211        getTriaInformation(triaId).n_cells == vec.size()))
00212  {
00213    printf("Error %d != %d || %d\n",vec.size(),getTriaInformation(triaId).n_vertices,getTriaInformation(triaId).n_cells);
00214  }
00215  assert(getTriaInformation(triaId).n_vertices == vec.size() ||
00216         getTriaInformation(triaId).n_cells == vec.size());
00217  assert(getFile() != -1);
00218 
00219 
00220   
00221   //Now write the dataset.
00222   //Create the dataset
00223   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00224   H5Tset_order(datatype,H5T_ORDER_LE);
00225 
00226   struct TriaInformation &triaInfo = getTriaInformation(triaId);
00227   if (!triaInfo.isCollapsed())  //pV==NULL means that we are writting the original mesh
00228   {
00229 
00230     hsize_t dim[]={vec.size()};
00231     hid_t dataspace = H5Screate_simple(1,dim,NULL);
00232     hid_t dataset = H5Dcreate1(getFile(), datasetName.c_str(), datatype, dataspace, H5P_DEFAULT);
00233     //Copy the data.
00234     H5Dwrite(dataset,H5T_NATIVE_DOUBLE,H5S_ALL,H5S_ALL,H5P_DEFAULT, vec.begin());
00235 
00236     //Close the resources.3
00237     H5Sclose(dataspace);
00238     H5Tclose(datatype);
00239     H5Dclose(dataset);
00240   }
00241   else
00242   {
00243     hsize_t dim[]={triaInfo.getActiveVertices().size()};
00244     hid_t dataspace = H5Screate_simple(1,dim,NULL);
00245     hid_t dataset = H5Dcreate1(getFile(), datasetName.c_str(), datatype, dataspace, H5P_DEFAULT);
00246     hsize_t memDim[]={vec.size()};
00247     hid_t memDataSpace = H5Screate_simple(1,memDim,NULL);
00248     H5Sselect_elements(memDataSpace,H5S_SELECT_SET,
00249                        triaInfo.getActiveVertices().size(),
00250                        &(triaInfo.getActiveVertices()[0]));
00251     //Copy the data.
00252     H5Dwrite(dataset,H5T_NATIVE_DOUBLE,memDataSpace,H5S_ALL,H5P_DEFAULT, vec.begin());
00253 
00254     //Close the resources.3
00255     H5Sclose(dataspace);
00256     H5Tclose(datatype);
00257     H5Sclose(memDataSpace);
00258     H5Dclose(dataset);
00259 
00260   }
00261   
00262   //Set attributes
00263   H5LTset_attribute_string(getFile(),datasetName.c_str(),"Triangulation",triaId.c_str());
00264   this->writeContextVariables(datasetName);
00265 
00266 
00267   
00268 }
00269 
00279 void HDF5OrthoWriter::collapseCellBlock(OrthoMesh::Cell_It &cell,unsigned blX,unsigned blY,unsigned blZ,StackVector<8,unsigned> &verts) 
00280  {
00281    
00282    blX--;
00283    blY--;
00284    blZ--;
00285    unsigned r1,r2,r3,r4,r5,r6,r7;
00286    OrthoMesh::Cell_It c1=cell;
00287   verts[VERTEX_000]=cell->vertex_index(VERTEX_000);
00288 
00289   r1=c1->advance(RIGHT_CELL,blX);
00290   verts[VERTEX_100]=c1->vertex_index(VERTEX_100);
00291 
00292   r2=c1->advance(BACK_CELL,blY);
00293   verts[VERTEX_110]=c1->vertex_index(VERTEX_110);
00294 
00295   r3=c1->advance(LEFT_CELL,blX);
00296   verts[VERTEX_010]=c1->vertex_index(VERTEX_010);
00297 
00298   c1=cell;
00299   r4=c1->advance(UP_CELL,blZ);
00300   verts[VERTEX_001]=c1->vertex_index(VERTEX_001);
00301 
00302   r5=c1->advance(RIGHT_CELL,blX);
00303   verts[VERTEX_101]=c1->vertex_index(VERTEX_101);
00304 
00305   r6=c1->advance(BACK_CELL,blY);
00306   verts[VERTEX_111]=c1->vertex_index(VERTEX_111);
00307 
00308   r7=c1->advance(LEFT_CELL,blX);
00309   verts[VERTEX_011]=c1->vertex_index(VERTEX_011);
00310 
00311   if (r1!=blX || r2 != blY || r3 != blX || r4 != blZ ||
00312       r5!=blX || r6 != blY || r7 != blX)
00313   {
00314     throw new Exception("HDF5DealWriter: Method reach the boundary while collapsing a block cell"); 
00315   }
00316 }
00317 
00318 
00319 
00320 typedef StackVector<8,unsigned> Verts;
00321 void HDF5OrthoWriter::writeCoarseMesh(string triaName,OrthoMesh &mesh,unsigned blkX,unsigned blkY,unsigned blkZ) 
00322 {
00323   if (blkX == 1 && blkY == 1 && blkZ == 1)
00324     return writeMesh(triaName,mesh);
00325 
00326   if (!mesh.getWells().empty())
00327   {
00328     return writeCoarseMeshWithWells(triaName,mesh,blkX,blkY,blkZ);
00329   }
00330   
00331   putMeshIntoDictionary(triaName,mesh.numVertices(),mesh.numCells());
00332   
00333   unsigned offX = blkX-1;
00334   unsigned offY = blkY-1;
00335   unsigned offZ = blkZ-1;
00336   unsigned dX,dZ,dY;
00337   unsigned nX,nY,nZ;
00338   nX = mesh.numElemX()/blkX + ((mesh.numElemX()%blkX == 0) ? 0 : 1);
00339   nY = mesh.numElemY()/blkY + ((mesh.numElemY()%blkY == 0) ? 0 : 1);
00340   nZ = mesh.numElemZ()/blkZ + ((mesh.numElemZ()%blkZ == 0) ? 0 : 1);
00341 
00342   
00343   std::vector<Verts> cells;
00344   cells.reserve(nX*nY*nZ);
00345   Verts verts;
00346   
00347   
00348   OrthoMesh::Cell_It cellZ = mesh.begin_cell();
00349   while(1) {
00350     OrthoMesh::Cell_It aux = cellZ;
00351     dZ = aux->advance(UP_CELL,offZ);
00352     OrthoMesh::Cell_It cellY=cellZ;
00353     while(1)
00354     {
00355       aux=cellY;
00356       dY=aux->advance(BACK_CELL,offY);
00357       OrthoMesh::Cell_It cellX=cellY;
00358       while(1)
00359       {
00360         aux=cellX;
00361         dX=aux->advance(RIGHT_CELL,offX);
00362         collapseCellBlock(cellX,dX+1,dY+1,dZ+1,verts);
00363         cells.push_back(verts);
00364         cellX->advance(RIGHT_CELL,dX);
00365         if (!cellX->advance(RIGHT_CELL))
00366           break;
00367       }
00368       cellY->advance(BACK_CELL,dY);
00369       if (!cellY->advance(BACK_CELL))
00370           break;
00371     }
00372     cellY->advance(UP_CELL,dZ);
00373     if (!cellZ->advance(UP_CELL))
00374       break;
00375   }
00376 
00377 
00378 
00379   
00380   std::vector<unsigned> vI(mesh.numVertices(),UINT_MAX);
00381   for(unsigned i=0;i<cells.size();i++)
00382   {
00383     for(unsigned j=0;j<8;j++)
00384       vI[cells[i][j]]=1;
00385   }
00386   unsigned acum=0;
00387   for (unsigned i=0;i<vI.size();i++)
00388   {
00389     if (vI[i] != UINT_MAX)
00390       vI[i]=acum++;
00391   }
00392   for(unsigned i=0;i<cells.size();i++)
00393     for(unsigned j=0;j<8;j++)
00394     {
00395       cells[i][j]=vI[cells[i][j]];
00396     }
00397   std::vector<hsize_t> vActiveVert;
00398   vActiveVert.reserve(acum);
00399   for(unsigned i=0;i<vI.size();i++)
00400   {
00401     if (vI[i] != UINT_MAX)
00402       vActiveVert.push_back(i);
00403   }
00404   getTriaInformation(triaName).setActiveVertices(vActiveVert);
00405 
00406   /*
00407   for(unsigned i=0;i<cells.size();i++)
00408   {
00409     printf("%d) ",i);
00410     for(unsigned j=0;j<8;j++)
00411       printf("%d ",cells[i][j]);
00412     printf("\n");
00413   }
00414   */
00415 
00416   //Criar Grupo
00417   char strPath[500];
00418   sprintf(strPath,"/Triangulations/%s",triaName.c_str());
00419   hid_t triaGroup=H5Gcreate1(getFile(),strPath,0);
00420   H5LTset_attribute_string(getFile(),strPath,"Type","GRID_3D");
00421 
00422   
00423   hsize_t dim[]={vActiveVert.size(),3};
00424   hid_t dataspace = H5Screate_simple(2,dim,NULL);
00425   hid_t datatype  = H5Tcopy(H5T_NATIVE_DOUBLE);
00426   H5Tset_order(datatype,H5T_ORDER_LE);
00427   hid_t dataset = H5Dcreate1(triaGroup, "vertices", datatype, dataspace, H5P_DEFAULT);
00428 
00429   int nSteps = vActiveVert.size()/HDF5_BUFF_SIZE;
00430   hsize_t count[]={HDF5_BUFF_SIZE,3};
00431   hsize_t offset[]={0,0};
00432   double buff[HDF5_BUFF_SIZE][3];
00433   unsigned k=0;
00434   for (int i=0;i<=nSteps;i++)
00435   {
00436     int MaxJ = (i==nSteps) ? (vActiveVert.size() % HDF5_BUFF_SIZE) : HDF5_BUFF_SIZE;
00437     count[0]=MaxJ;
00438     for (int j=0;j<MaxJ;j++)
00439     {
00440       Point3D P = mesh.get_vertice(vActiveVert[k])->getPoint();
00441       buff[j][0]= P[0];
00442       buff[j][1]= P[1];
00443       buff[j][2]= P[2];
00444       k++;
00445     }
00446     writeInChunks(dataset,offset,count,buff);
00447     offset[0]+=MaxJ;
00448   }
00449   H5Sclose(dataspace);
00450   H5Tclose(datatype);
00451   H5Dclose(dataset);
00452 
00453   dim[0]=cells.size();
00454   dim[1]=8;
00455   dataspace = H5Screate_simple(2,dim,NULL);
00456   datatype  = H5Tcopy(H5T_NATIVE_UINT);
00457   H5Tset_order(datatype,H5T_ORDER_LE);
00458   dataset = H5Dcreate1(triaGroup, "connections", datatype, dataspace, H5P_DEFAULT);
00459   H5Dwrite(dataset,H5T_NATIVE_UINT,H5S_ALL,H5S_ALL,H5P_DEFAULT, &(cells[0][0]));
00460 
00461   H5Sclose(dataspace);
00462   H5Tclose(datatype);
00463   H5Dclose(dataset);
00464   H5Gclose(triaGroup);
00465   
00466     
00467 
00468 }
00469 
00470 
00471 
00472 
00473 void HDF5OrthoWriter::createGroup(std::string str) 
00474 {
00475   hid_t aux;
00476   H5Eset_auto2(H5E_DEFAULT,NULL,stdout);
00477   aux = H5Gopen1(getFile(),str.c_str());
00478   H5Eset_auto2(H5E_DEFAULT,(herr_t (*)(hid_t,void*))   &H5Eprint,stdout);
00479   if (aux >= 0)
00480   {
00481     H5Gclose(aux);
00482     return;
00483   }
00484   else
00485   {
00486     H5Gclose(H5Gcreate1(getFile(),str.c_str(),H5P_DEFAULT));
00487   }
00488 }
00489 
00490 
00491 
00492 
00493 void HDF5OrthoWriter::writeCoarseMeshWithWells(string triaName,OrthoMesh &mesh,unsigned blkX,unsigned blkY,unsigned blkZ) 
00494 {
00495   if (blkX == 1 && blkY == 1 && blkZ == 1)
00496     return writeMesh(triaName,mesh);
00497 
00498   unsigned nX = mesh.numElemX();
00499   unsigned nY = mesh.numElemY();
00500   unsigned nZ = mesh.numElemZ();
00501 
00502   blkX=adjustBlockSize(nX,blkX); 
00503   blkY=adjustBlockSize(nY,blkY);
00504   blkZ=adjustBlockSize(nZ,blkZ);
00505   printf("Adjusting blocks steps for mesh %s\n",triaName.c_str());
00506   printf("Block X: %u\n",blkX);
00507   printf("Block Y: %u\n",blkY);
00508   printf("Block Z: %u\n",blkZ);
00509  
00510   OrthoMesh newMesh(mesh.getP(),mesh.getQ(),nX/blkX,nY/blkY,nZ/blkZ);
00511   VecWellInfo wells=mesh.getWells();
00512   for (unsigned i=0;i<wells.size();i++)
00513   wells[i].adjustBoundaryWithGrid(newMesh.getDX());
00514   newMesh.putWells(wells);
00515   writeMesh(triaName,newMesh);
00516 
00517   getTriaInformation(triaName).n_vertices=mesh.numVertices();
00518   getTriaInformation(triaName).n_cells=mesh.numCells();
00519   std::vector<hsize_t> map(newMesh.numVertices(),UINT_MAX);
00520   for (OrthoMesh::Vertex_It vIt = newMesh.begin_vertice();vIt != newMesh.end_vertice();vIt++)
00521   {
00522     map[vIt->index()]=mesh.getNearestVertice(vIt->getPoint())->index();
00523   }
00524   getTriaInformation(triaName).setActiveVertices(map);
00525 }
00526 
00527 
00528 unsigned HDF5OrthoWriter::adjustBlockSize(unsigned nElems,unsigned blkSize) 
00529 {
00530   for (unsigned Itry=blkSize;Itry<=nElems;Itry++)
00531   {
00532     if (nElems%Itry == 0)
00533       return Itry;
00534   }
00535   return nElems;
00536 }
00537 
00538 
00539 void HDF5OrthoWriter::close() 
00540 {
00541   std::string file=getFileName();
00542   HDF5Writer::close();
00543   printf("File:   %s\n",file.c_str());
00544   if (!file.empty())
00545   {
00546     printf("Writing xdmf\n");
00547     XdmfWriter xdmf;
00548     xdmf.writeXdmf(file);
00549   }
00550  
00551 }
00552 
00553 
00554 void HDF5OrthoWriter::supressOutput() 
00555 {
00556   _supress_output=true;
00557 }