src/HDF5ParallelReader.cpp

Go to the documentation of this file.
00001 #include "HDF5ParallelReader.h"
00002 using namespace std;
00003 
00004 HDF5ParallelReader::~HDF5ParallelReader() {
00005     delete[] _tetids;
00006     delete[] _matids;
00007     delete[] _tet_partition;
00008     delete[] _tet_global_id;
00009     delete[] _point_global_id;
00010     delete[] _point_coords;
00011 }
00012 
00013 HDF5ParallelReader::HDF5ParallelReader(const char* inputFile, MeshBuilder* builder, MPI_Comm comm)
00014     : _inputFile(inputFile), _builder(builder), _comm(comm) {
00015     MPI_Comm_size(_comm, &_mpi_size);
00016     MPI_Comm_rank(_comm, &_mpi_rank);  
00017     _info  = MPI_INFO_NULL;
00018 }
00019 
00020 double HDF5ParallelReader::seconds() {
00021     static long clk_tck = sysconf(_SC_CLK_TCK);
00022     struct tms buf;
00023     return (double)(times(&buf)) / clk_tck;
00024 }
00025 
00026 void HDF5ParallelReader::print_tet_global_id() {
00027     cout << "proc " << _mpi_rank << ": ";  
00028     for (int i = 0; i < _numtets_local; i ++) {
00029         if (_tet_global_id[i] < 10)
00030             cout << " ";
00031         cout << _tet_global_id[i] << " ";
00032     }
00033     cout << endl;
00034 }
00035 
00036 void HDF5ParallelReader::read_tets() {
00037 
00038     hsize_t     dimsf[2];
00039     hsize_t     count[2];
00040 #if defined(USE_HDF5_HSSIZE)
00041     hssize_t    offset[2];
00042 #else
00043     hsize_t     offset[2];
00044 #endif
00045     herr_t      status;
00046    
00047     hid_t dset_id = H5Dopen(_file_id, "/tetids");
00048     assert(dset_id >= 0);
00049 
00050 // Create the dataspace for the dataset.
00051     hid_t x = H5Dget_space(dset_id);
00052     H5Sget_simple_extent_dims(x, dimsf, NULL);
00053     hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
00054     _numtets_global = dimsf[0];
00055 
00056 // Set offset[2] which describes where in the file to start reading
00057     _numtets_local = _numtets_global / _mpi_size;
00058     int remainder = _numtets_global % _mpi_size;
00059 
00060     offset[0] = _mpi_rank * _numtets_local;
00061     if (_mpi_rank < remainder) {
00062         offset[0] += _mpi_rank;
00063         _numtets_local += 1;
00064     } else 
00065         offset[0] += remainder;
00066     offset[1] = 0;
00067 
00068     count[0] = _numtets_local;
00069     count[1] = 4; //Each tet has 4 nodes
00070     hid_t memspace = H5Screate_simple(2, count, NULL);
00071     
00072 // Select hyperslab in the file.
00073     filespace = H5Dget_space(dset_id);
00074     assert(filespace >= 0);
00075 
00076     H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
00077     assert(filespace >= 0);
00078 
00079     _tetids = new idxtype[_numtets_local * 4];
00080 
00081 // Create property list for collective dataset write.
00082     _plist_id = H5Pcreate(H5P_DATASET_XFER);
00083     H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);
00084 
00085     status = H5Dread(dset_id, INTMEMTYPE, memspace, filespace, _plist_id, _tetids);
00086     assert(status >= 0);
00087 
00088 // store global numbers of local tets
00089     _tet_global_id = new int[_numtets_local];
00090     for (int i = 0; i < _numtets_local; i ++)
00091         _tet_global_id[i] = offset[0] + i;
00092    
00093 // Close HDF5 stuff
00094     H5Dclose(dset_id);
00095     H5Sclose(filespace);
00096     H5Sclose(memspace);
00097     
00098     rInfoAll("Number of tets: %d local, %d global", _numtets_local, _numtets_global);
00099 }
00100 
00101 
00102 void HDF5ParallelReader::read_mats() {
00103 
00104     hsize_t     dimsf[2];
00105     hsize_t     count[2];
00106 #if defined(USE_HDF5_HSSIZE)
00107     hssize_t    offset[2];
00108 #else
00109     hsize_t     offset[2];
00110 #endif
00111     herr_t      status;
00112    
00113     hid_t dset_id = H5Dopen(_file_id, "/matids");
00114     assert(dset_id >= 0);
00115 
00116     // Create the dataspace for the dataset.
00117     hid_t x = H5Dget_space(dset_id);
00118     H5Sget_simple_extent_dims(x, dimsf, NULL);
00119     hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
00120     // _numtets_global = dimsf[0];                 OSCAR: As many materials as tetrahedra ...
00121 
00122     // Set offset[2] which describes where in the file to start reading
00123     _numtets_local = _numtets_global / _mpi_size;
00124     int remainder = _numtets_global % _mpi_size;
00125 
00126     offset[0] = _mpi_rank * _numtets_local;
00127     if (_mpi_rank < remainder) {
00128         offset[0] += _mpi_rank;
00129         _numtets_local += 1;
00130     } else 
00131         offset[0] += remainder;
00132     offset[1] = 0;
00133 
00134     count[0] = _numtets_local;
00135     count[1] = 1; //Each tet has 1 material id
00136     hid_t memspace = H5Screate_simple(2, count, NULL);
00137     
00138     // Select hyperslab in the file.
00139     filespace = H5Dget_space(dset_id);
00140     assert(filespace >= 0);
00141 
00142     H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
00143     assert(filespace >= 0);
00144 
00145     _matids = new idxtype[_numtets_local];
00146 
00147     // Create property list for collective dataset write.
00148     _plist_id = H5Pcreate(H5P_DATASET_XFER);
00149     H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);
00150 
00151     status = H5Dread(dset_id, INTMEMTYPE, memspace, filespace, _plist_id, _matids);
00152     assert(status >= 0);    
00153 
00154     rDebug("Number of material ids: %d", _numtets_local);
00155 
00156     // Close HDF5 stuff
00157     H5Dclose(dset_id);
00158     H5Sclose(filespace);
00159     H5Sclose(memspace);
00160 }
00161 
00162 
00163 
00164 
00165 void HDF5ParallelReader::partition() {
00166 // All processors will build data to call ParMetis_V3_PartMeshKway
00167 
00168 // elmdist specifies distribution of elements
00169     int remainder = _numtets_global % _mpi_size;
00170     idxtype* elmdist = new idxtype[_mpi_size + 1];
00171     int off = 0;
00172     for (int i = 0; i < _mpi_size + 1; i ++) {
00173         elmdist[i] = off;
00174         off += _numtets_global / _mpi_size;
00175         if (i < remainder) 
00176             off += 1;
00177     }
00178 
00179 // eind specifies at what index of eptr the elements that make up node i starts
00180     idxtype* eptr = new idxtype[_numtets_local + 1];
00181     for (int i = 0; i < _numtets_local + 1; i ++) {
00182         eptr[i] = 4 * i;
00183     }
00184     
00185 // with elmdist and eptr as above, eind is equal to data!
00186 
00187 // More ParMetis arguments
00188     int wgtflag = 0; //No weights
00189     int numflag = 0; //C-style numbering starts from 0
00190 
00191     int ncon = 1; //number of weights on each vertex
00192     int ncommonnodes = 3; //no of common nodes to be considered connected
00193     int nparts = _mpi_size; //number of subdomains in output
00194 
00195     float* tpwgts = new float[nparts]; //fraction of nodes to each processor
00196     for (int i = 0; i < nparts; i ++)
00197         tpwgts[i] = 1.0 / (float) nparts;
00198 
00199     float ubvec = 1.05; //imbalance acceptance
00200 
00201     int options[10];
00202     options[0] = 0; //1 means debug
00203     options[1] = 3; //level of output between 1..127
00204     options[2] = 15; //random seed
00205 
00206 //Output pointers
00207     _tet_partition = new idxtype[_numtets_local]; //output of ParMETIS
00208     int edgecut; //The number of edges cut
00209 
00210 //Call ParMETIS
00211     ParMETIS_V3_PartMeshKway(elmdist, eptr, _tetids, //OK data may not be idx
00212                              NULL, &wgtflag, //OK
00213                              &numflag, &ncon, //OK
00214                              &ncommonnodes, &nparts, //OK
00215                              tpwgts, &ubvec, //OK
00216                              options, &edgecut, //OK
00217                              _tet_partition, &_comm); //OK
00218 
00219 //Some output
00220 //     if (_mpi_rank == 0)
00221 //      cout << "edgecut: " << edgecut << endl;
00222 
00223 //     if (edgecut < 100) {
00224 //        cout << "partition of processor " << _mpi_rank << ": ";
00225 //     for (int i = 0; i < _numtets_local; i ++)
00226 //      cout << _tet_partition[i] << " ";
00227 //     cout << endl;
00228 //     }
00229 
00230     delete[] tpwgts;
00231     delete[] elmdist;
00232     delete[] eptr;
00233 }
00234 
00235 void HDF5ParallelReader::sort_buffer(int* sendcnts) {
00236 
00237 //This sorts the sendbuffer by the target proc.
00238 //A bin/bucket sort where the sizes of the buckets is known
00239 
00240     int* next_pos = new int[_mpi_size];
00241     next_pos[0] = 0;
00242     for (int i = 1; i < _mpi_size; i ++) 
00243         next_pos[i] = next_pos[i - 1] + sendcnts[i - 1];
00244 
00245     int* sorted_buffer = new int[4 * _numtets_local];
00246     int* sorted_matbuffer = new int[_numtets_local];
00247     int* sorted_partition = new int[_numtets_local];
00248     int* sorted_tet_global_id = new int[_numtets_local];
00249 
00250     for (int i = 0; i < _numtets_local; i ++) {
00251       for (int j = 0; j < 4; j ++){
00252         sorted_buffer[next_pos[_tet_partition[i]] + j] = _tetids[4 * i + j];
00253       }
00254       sorted_matbuffer[next_pos[_tet_partition[i]]/4] = _matids[i];             // OSCAR: Everything you do with tets you do with mats (div 4)
00255         
00256       sorted_partition[next_pos[_tet_partition[i]] / 4] = _tet_partition[i];
00257       sorted_tet_global_id[next_pos[_tet_partition[i]] / 4] = _tet_global_id[i];
00258       next_pos[_tet_partition[i]] += 4;
00259     }
00260 
00261     delete[] _tetids;
00262     _tetids = sorted_buffer;
00263     delete[] _matids;
00264     _matids = sorted_matbuffer;
00265     delete[] _tet_global_id;
00266     _tet_global_id = sorted_tet_global_id;
00267 
00268 //     if (_numtets_local < 20) {
00269 //      cout << "sorted partition of processor " << _mpi_rank << ": ";
00270 //      for (int i = 0; i < _numtets_local; i ++)
00271 //          cout << sorted_partition[i] << " ";
00272 //      cout << endl;
00273 //     }
00274 
00275     delete[] next_pos;
00276     delete[] sorted_partition;    
00277 
00278 }
00279 
00280 void HDF5ParallelReader::distribute() {
00281 /* distribute just tets, other things will be read after distribution to reduce communication
00282    there is no need for tet numbering as long as node/tet pairs remain intact */
00283 
00284 //Count how many POINTS to send to each proc (sendcnts)
00285     int* sendcnts = new int[_mpi_size];
00286     for (int i = 0; i < _mpi_size; i ++)
00287         sendcnts[i]=0;
00288     for (int i = 0; i < _numtets_local; i++)
00289         sendcnts[_tet_partition[i]] += 4;
00290 
00291  //    cout << "Proc " << _mpi_rank << " sends ";
00292 //     for (int i = 0; i < _mpi_size; i ++)
00293 //      cout << sendcnts[i] << " ";
00294 //     cout << " to the others" << endl;
00295    
00296 //Use MPI_Alltoall to find no of points proc i will recieve from proc others (recvcnts)
00297     int* recvcnts = new int[_mpi_size];
00298     MPI_Alltoall (sendcnts, 1, MPI_INT, recvcnts, 1, MPI_INT, _comm);
00299 
00300  //    cout << "Proc " << _mpi_rank << " recieves ";
00301 //     for (int i = 0; i < _mpi_size; i ++)
00302 //      cout << recvcnts[i] << " ";
00303 //     cout << " from the others" << endl;
00304 
00305 //Construct send displacements (sdispls)
00306     int* sdispls = new int[_mpi_size];
00307     sdispls[0] = 0;
00308     for (int i = 1; i < _mpi_size; i ++)
00309         sdispls[i] = sdispls[i - 1] + sendcnts[i - 1];
00310     
00311 //Construct recieve displacements (rdispls)
00312     int* rdispls = new int[_mpi_size];
00313     rdispls[0] = 0;
00314     for (int i = 1; i < _mpi_size; i ++)
00315         rdispls[i] = rdispls[i - 1] + recvcnts[i - 1];
00316 
00317 //Problem: The sendbuffer is not sequential
00318 //Must sort the sendbuffer.
00319     sort_buffer(sendcnts);
00320     //print_tet_global_id();
00321 
00322     int recvbuf_size = 0;
00323     for (int i = 0; i < _mpi_size; i ++)
00324         recvbuf_size += recvcnts[i];
00325     int* recvbuf = new int[recvbuf_size];
00326 
00327     
00328     MPI_Alltoallv (_tetids, sendcnts, sdispls, MPI_INT, recvbuf, recvcnts, rdispls, MPI_INT, _comm );
00329     delete[] _tetids;
00330 
00331     _tetids = recvbuf;
00332     _numtets_local = recvbuf_size / 4;
00333     
00334     // Send global materials and numbering _tet_global_id. Reuse sendcnts, recvcnts, sdispls, rdispls and recvbuf_size
00335     for (int i = 0; i < _mpi_size; i++) {
00336         sendcnts[i] /= 4;
00337         recvcnts[i] /= 4;
00338         sdispls[i] /= 4;
00339         rdispls[i] /= 4;
00340     }
00341     recvbuf_size /=4;
00342     
00343     recvbuf = new int[recvbuf_size];
00344     MPI_Alltoallv (_matids, sendcnts, sdispls, MPI_INT, recvbuf, recvcnts, rdispls, MPI_INT, _comm );
00345     delete[] _matids;
00346     _matids = recvbuf;
00347 
00348 
00349     int* distributed_tet_global_id = new int[recvbuf_size];
00350     MPI_Alltoallv (_tet_global_id, sendcnts, sdispls, MPI_INT, distributed_tet_global_id, recvcnts, rdispls, MPI_INT, _comm);
00351     delete[] _tet_global_id;
00352     _tet_global_id = distributed_tet_global_id;
00353     
00354 
00355    //  if (_mpi_rank == 0) {
00356 //      for (int i = 0; i < _numtets_local; i ++) {
00357 //          for (int j = 0; j < 4; j ++)
00358 //              cout << _tetids[4 * i + j] << " ";
00359 //          cout << endl;
00360 //      }
00361 //     }
00362     delete[] sendcnts;
00363     delete[] recvcnts;
00364     delete[] sdispls;
00365     delete[] rdispls;
00366 }
00367 
00368 void HDF5ParallelReader::generate_point_global_id(int* slab_offsets, int* slab_lengths, int slab_count) {
00369     _point_global_id = new int[_numpoints_local];
00370     int counter = 0;
00371     
00372     for(int i = 0; i < slab_count; i ++) {
00373         for(int j = 0; j < slab_lengths[i]; j ++) {
00374             _point_global_id[counter] = slab_offsets[i] + j;
00375             
00376             //cout << _point_global_id[counter] << " ";
00377             
00378             counter += 1;
00379             
00380             
00381         }
00382     }
00383 }
00384 
00385 void HDF5ParallelReader::read_points() {
00386     
00387     hid_t dset_id = H5Dopen(_file_id, "/coords");
00388     assert(dset_id >= 0);
00389 
00390 // Create the dataspace for the dataset.
00391     hsize_t dimsf[2]; 
00392     herr_t status;
00393  
00394     hid_t x = H5Dget_space(dset_id);
00395     H5Sget_simple_extent_dims(x, dimsf, NULL);
00396 //    hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
00397     _numpoints_global = dimsf[0];
00398 
00399    
00400 //What nodes should proc i read
00401 //Assume no of points known
00402     bool* read_this_node = new bool[_numpoints_global];
00403     _numpoints_local = 0;
00404     for (int i = 0; i < _numpoints_global; i ++)
00405         read_this_node[i] = 0;
00406     for (int i = 0; i < 4 * _numtets_local; i ++) {
00407         if (read_this_node[_tetids[i]] == 0)
00408             _numpoints_local += 1;
00409         read_this_node[_tetids[i]] = 1;
00410     }
00411 
00412 //Use read_this_node[] to generate slabs
00413     int slab_count = 0;
00414     bool prev = 0;
00415     for (int i = 0; i < _numpoints_global; i ++) {
00416         if (read_this_node[i] != prev) {
00417             prev = read_this_node[i];
00418             if (read_this_node[i] == 1)
00419                 slab_count += 1;
00420         }
00421     }
00422 
00423     int* slab_offsets = new int[slab_count];
00424     int* slab_lengths = new int[slab_count];
00425     int counter = 0;
00426     int start = 0;
00427     prev = 0;
00428 //And go through read_this_node again
00429     for (int i = 0; i < _numpoints_global; i ++) {
00430         if (read_this_node[i] != prev) {
00431             prev = read_this_node[i];
00432             if (read_this_node[i] == 1) {
00433                 start = i;
00434                 slab_offsets[counter] = start;
00435                 slab_lengths[counter] = _numpoints_global - start;
00436             } else {
00437                 slab_lengths[counter] = i - start;
00438                 counter ++;
00439             }
00440         }
00441     }
00442     
00443     //Generate global_local point index map
00444     generate_point_global_id(slab_offsets, slab_lengths, slab_count);
00445 
00446 //     cout << "Proc " << _mpi_rank << " reads: ";
00447 // //     for (int i = 0; i < _numpoints_global; i ++)
00448 // //   cout << read_this_node[i];
00449 // //     for (int i = 0; i < slab_count; i ++) 
00450 // //   cout << " " << slab_offsets[i] << "-" << slab_offsets[i] + slab_lengths[i] - 1;
00451 //     cout << " totally " << slab_count << " slabs" << endl;
00452 
00453  //    cout << "Number of points globally: " << _numpoints_global << endl;
00454 //     cout << "Number of points locally: " << _numpoints_local << endl;
00455 //     cout << "Number of slabs: " << slab_count << endl;
00456 
00457     hsize_t count[2];    
00458     count[0] = _numpoints_local;
00459     count[1] = 3;
00460     hid_t memspace = H5Screate_simple(2, count, NULL);
00461     assert (memspace >= 0);
00462 
00463 //Construct hyperslab in file
00464     hid_t filespace = H5Dget_space(dset_id);
00465     assert (filespace >= 0);
00466 
00467 #if defined(USE_HDF5_HSSIZE)
00468     hssize_t sstart[2];
00469 #else
00470     hsize_t sstart[2];
00471 #endif
00472     hsize_t sblock[2];
00473     hsize_t sstride[2];
00474     hsize_t scount[2];
00475     sblock[0] = 1; sblock[1] = 1;
00476     sstride[0] = 1; sstride[1] = 1;
00477     sstart[1] = 0;
00478     scount[1] = 3;
00479     for (int i = 0; i < slab_count; i ++) {
00480         sstart[0] = slab_offsets[i]; 
00481         scount[0] = slab_lengths[i];
00482         if (i == 0)
00483             status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, sstart, sstride, scount, sblock);
00484         else
00485             status = H5Sselect_hyperslab(filespace, H5S_SELECT_OR, sstart, sstride, scount, sblock);
00486         
00487         assert (status >= 0);
00488      }
00489 
00490     H5Sget_simple_extent_dims(filespace, dimsf, NULL);
00491   //   if (_mpi_rank == 0) {
00492 //      cout << "Size of filespace: " << dimsf[0] << "x" << dimsf[1] << endl;
00493 //      cout << "Slab union # elems: " << H5Sget_select_npoints(filespace) << endl;
00494 //     }
00495 
00496 
00497     H5Sget_simple_extent_dims(memspace, dimsf, NULL);
00498   //   if (_mpi_rank == 0)
00499 //      cout << "Size of memspace: " << dimsf[0] << "x" << dimsf[1] << endl;
00500 
00501 
00502 //Create property list for collective dataset read.
00503     _plist_id = H5Pcreate(H5P_DATASET_XFER);
00504     H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);
00505 
00506     _point_coords = new double[3 * _numpoints_local];
00507     status = H5Dread(dset_id, REALMEMTYPE, memspace, filespace, _plist_id, _point_coords);
00508     assert(status >= 0);
00509 
00510 // Close HDF5 stuff
00511     H5Dclose(dset_id);
00512     H5Sclose(filespace);
00513     H5Sclose(memspace);
00514 
00515     delete[] slab_offsets;
00516     delete[] slab_lengths;
00517     delete [] read_this_node;
00518 }
00519 
00520 void HDF5ParallelReader::read_boundary() {
00521 
00522 //So first read in all boundary points
00523     hsize_t     dimsf[2];                 /* dataset dimensions */
00524     herr_t      status;
00525    
00526     hid_t dset_id = H5Dopen(_file_id, "/surface");
00527     assert(dset_id >= 0);
00528 
00529 // Create the dataspace for the dataset.
00530     hid_t x = H5Dget_space(dset_id);
00531     H5Sget_simple_extent_dims(x, dimsf, NULL);
00532     hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
00533     _numbfaces_global = dimsf[0];
00534 
00535     int* allbfaces = new idxtype[_numbfaces_global * 4];
00536 
00537 // Create property list for collective dataset write.
00538     _plist_id = H5Pcreate(H5P_DATASET_XFER);
00539     H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);
00540 
00541     status = H5Dread(dset_id, INTMEMTYPE, H5S_ALL, filespace, _plist_id, allbfaces);
00542     assert(status >= 0);
00543         
00544 // Store local boundary faces, discard others.
00545    _numbfaces_local = 0;
00546    int nof_sym = 0; //Count number of symmetry planes.
00547 
00548     for (int i = 0; i < _numbfaces_global; i ++) {
00549         try {
00550             _builder->set_bc(allbfaces[4*i+1], allbfaces[4*i+2], allbfaces[4*i+3], allbfaces[4*i]);
00551             _numbfaces_local += 1;
00552         } catch (std::runtime_error err) {}//Face was not local
00553         
00554         if (allbfaces[4 * i] > nof_sym)
00555             nof_sym = allbfaces[4 * i];
00556     }
00557     
00558     //cout << "Proc " << _mpi_rank <<": " << _numbfaces_local << " faces, " << nof_sym << " symmetry planes" << endl;
00559     
00560     _builder->finalize_bc(nof_sym);
00561 
00562 // Close HDF5 stuff
00563     H5Dclose(dset_id);
00564     H5Sclose(filespace);
00565     
00566     delete[] allbfaces;
00567 }
00568 
00569 void HDF5ParallelReader::export_dist(char* file) {
00570 
00571     std::ofstream of;
00572     int* recvcnts;
00573 
00574 //Proc 0 must read all points from file then output them
00575     if (_mpi_rank == 0) {
00576         of.open (file);
00577         assert (of.is_open());
00578         of.precision (6);
00579 
00580         of << "# vtk DataFile Version 2.0" << endl;                      // first line of a vtk file
00581         of << "generated using hdf5_parallel_read::export_dist" << endl;   // header
00582         of << "ASCII" << endl << endl;                                   // file format
00583         of << "DATASET UNSTRUCTURED_GRID" << endl;
00584         of << "POINTS " << _numpoints_global << " float" << endl;
00585         of << scientific; //????????
00586 
00587         hid_t dataset = H5Dopen(_file_id, "/coords");
00588         assert (dataset >= 0);
00589 
00590         double* nodes = new double[3 * _numpoints_global];
00591         herr_t status = H5Dread(dataset, REALMEMTYPE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nodes);
00592         assert (status >= 0);
00593 
00594         for (int i = 0; i < _numpoints_global; i ++)
00595             of << nodes[3 * i] << " " << nodes[3 * i + 1] << " " << nodes[3 * i + 2] << endl;
00596         of << endl;
00597 
00598         H5Dclose (dataset);
00599         delete[] nodes;
00600 
00601 //Read all tets from file and output them
00602         of << "CELLS " << _numtets_global << " " << 5 * _numtets_global << endl;
00603     }
00604 
00605     recvcnts = new int[_mpi_size];
00606     MPI_Gather (&_numtets_local, 1, MPI_INT, 
00607                 recvcnts, 1, MPI_INT, 
00608                 0, _comm);
00609     int* displs = new int[_mpi_size];
00610     displs[0] = 0;
00611     recvcnts[0] *= 4;
00612     for (int i = 1; i < _mpi_size; i ++) {
00613         recvcnts[i] *= 4;
00614         displs[i] = displs[i - 1] + recvcnts[i - 1];
00615     }
00616 //     for (int i = 0; i < _mpi_size; i ++) {
00617 //      if (_mpi_rank == 0)
00618 //          cout << recvcnts[i] << " "<< displs[i] << endl;
00619 //     }
00620 
00621     int* total_partition = new int[4 * _numtets_global];
00622     MPI_Gatherv (_tetids, 4 * _numtets_local, MPI_INT, 
00623                  total_partition, recvcnts, displs,
00624                  MPI_INT, 0, _comm);
00625 
00626     delete[] displs;
00627 
00628     if (_mpi_rank == 0) {
00629 
00630         for (int i = 0; i < _numtets_global; i ++) {
00631             of << "4 ";
00632             for (int j = 0; j < 4; j ++)
00633                 of << total_partition[4 * i + j] << " ";
00634             of << endl;
00635         }
00636         of << endl;
00637 
00638         of << "CELL_TYPES " << _numtets_global << endl;
00639         for (int i = 0; i < _numtets_global; i ++)
00640             of << 10 << endl;
00641         of << endl;
00642         of << endl;
00643     }
00644 
00645 // 
00646 
00647 // Proc 0 writes colormap
00648     if (_mpi_rank == 0) {
00649 
00650         of << "CELL_DATA " << _numtets_global << endl << endl;
00651         of << "SCALARS process_no int 1" << endl;
00652         of << "LOOKUP_TABLE default" << endl;
00653 
00654         for (int i = 0; i < _mpi_size; i ++)
00655             for (int j = 0; j < recvcnts[i] / 4; j ++)
00656                 of << i << endl;
00657     }
00658 
00659 
00660     
00661 }
00662 //
00663 void HDF5ParallelReader::open_file() {
00664 // Set up file access property list with parallel I/O access
00665     _plist_id = H5Pcreate(H5P_FILE_ACCESS); //Property list identifier
00666     H5Pset_fapl_mpio(_plist_id, _comm, _info);   
00667 // Create a new file collectively and release property list identifier.
00668     _file_id = H5Fopen(_inputFile, H5F_ACC_RDONLY, _plist_id);
00669     assert(_file_id >= 0);
00670     H5Pclose(_plist_id);
00671 }
00672 
00673 void HDF5ParallelReader::close_file() {
00674     H5Pclose(_plist_id);
00675     H5Fclose(_file_id);
00676 }
00677     
00678 void HDF5ParallelReader::read() {
00679     open_file();
00680     
00681     read_tets();
00682 
00683     read_mats();
00684     
00685     partition();
00686     
00687     distribute();
00688     
00689     read_points();
00690     
00691     _builder->init_coord(_numpoints_local);
00692     for (int i = 0; i < _numpoints_local; i ++)
00693         _builder->set_coord(_point_global_id[i], _point_coords[3 * i], _point_coords[3 * i + 1], _point_coords[3 * i + 2]);
00694         
00695     _builder->finalize_coord();
00696         
00697     _builder->init_tet(_numtets_local);       
00698     for (int i = 0; i < _numtets_local; i ++) {
00699       _builder->set_tet(i, _tetids[4 * i], _tetids[4 * i + 1], _tetids[4 * i + 2], _tetids[4 * i + 3], _matids[i]);
00700     }
00701 
00702     _builder->finalize_tet();
00703     
00704     read_boundary();
00705     
00706     close_file();
00707 }
00708 
00709 int HDF5ParallelReader::main (int argc, char **argv) {
00710 
00711 #ifdef __GNUC__
00712     std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
00713 #endif
00714 
00715     MPI_Init(&argc, &argv);    
00716     
00717     char* inputFile;
00718     
00719     if (argc != 2) {
00720         if (_mpi_rank == 0) 
00721             cout << "Usage: hdf5_parallel_read inputfile.h5" << endl;
00722         return 0;
00723     } else {
00724         inputFile = argv[1];
00725     }
00726     
00727     
00728     MPI_Comm comm = MPI_COMM_WORLD;
00729     
00730     TetMeshBuilder* builder = new TetMeshBuilder(comm);
00731     HDF5ParallelReader* reader = new HDF5ParallelReader(inputFile, builder, comm);
00732     
00733     double time;
00734     
00735 
00736     
00737 
00738     //time = seconds();
00739     open_file(); 
00740     //cout << "Proc: " << _mpi_rank << " open_file: " << seconds() - time << endl;
00741 
00742     time = seconds();
00743     read_tets();
00744     read_mats();
00745     //cout << "Proc: " << _mpi_rank << " read_tets: " << seconds() - time << endl;
00746     //print_tet_global_id();  
00747     
00748     time = seconds();
00749     partition(); 
00750     //cout << "Proc: " << _mpi_rank << " partition: " << seconds() - time << endl;
00751 
00752     time = seconds();
00753     reader->distribute();
00754     //cout << "Proc: " << _mpi_rank << " distribute: " << seconds() - time << endl;
00755     //print_tet_global_id();  
00756  
00757     time = seconds();
00758     reader->read_points();
00759     //cout << "Proc: " << _mpi_rank << " read_points: " << seconds() - time << endl;
00760 
00761     
00762     _builder->init_coord(_numpoints_local);
00763     
00764      for (int i = 0; i < _numpoints_local; i ++)
00765          _builder->set_coord(_point_global_id[i], _point_coords[3 * i], _point_coords[3 * i + 1], _point_coords[3 * i + 2]);
00766      
00767     _builder->finalize_coord();
00768     //
00769     
00770     _builder->init_tet(_numtets_local);       
00771         
00772     for (int i = 0; i < _numtets_local; i ++) {
00773         //cout << _mpi_rank << ":" << i;
00774         //cout << ": (" << _tetids[4 * i] << "," << _tetids[4 * i + 1] << "," << _tetids[4 * i + 2] << "," << _tetids[4 * i + 3] << ") " << endl;
00775       _builder->set_tet(i, _tetids[4 * i], _tetids[4 * i + 1], _tetids[4 * i + 2], _tetids[4 * i + 3], _matids[i]);
00776     }
00777 
00778     _builder->finalize_tet();
00779     
00780      time = seconds();
00781     read_boundary();
00782     //cout << "Proc: " << _mpi_rank << " read_boundary: " << seconds() - time << endl;
00783 
00784     time = seconds();
00785     export_dist("out.vtk");
00786     //cout << "Proc: " << _mpi_rank << " export_dist: " << seconds() - time << endl;
00787     
00788     
00789     
00790     
00791 //Close/release resources.
00792     close_file();
00793     
00794 
00795     mesh::ParallelTetMesh* mesh = dynamic_cast<mesh::ParallelTetMesh*> (builder->get_mesh());
00796     
00797     
00798     
00799     delete _builder;
00800     
00801     
00802     delete mesh;
00803     
00804     delete _tetids;
00805     delete _tet_partition;
00806 
00807     MPI_Finalize();
00808     
00809     return 0;
00810 }

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7