src/HDF5ParallelReader.cpp

Go to the documentation of this file.

#include "HDF5ParallelReader.h"
using namespace std;

HDF5ParallelReader::~HDF5ParallelReader() {
    delete[] _tetids;
    delete[] _matids;
    delete[] _tet_partition;
    delete[] _tet_global_id;
    delete[] _point_global_id;
    delete[] _point_coords;
}

HDF5ParallelReader::HDF5ParallelReader(const char* inputFile, MeshBuilder* builder, MPI_Comm comm)
    : _inputFile(inputFile), _builder(builder), _comm(comm) {
    MPI_Comm_size(_comm, &_mpi_size);
    MPI_Comm_rank(_comm, &_mpi_rank);  
    _info  = MPI_INFO_NULL;
}

double HDF5ParallelReader::seconds() {
    static long clk_tck = sysconf(_SC_CLK_TCK);
    struct tms buf;
    return (double)(times(&buf)) / clk_tck;
}

void HDF5ParallelReader::print_tet_global_id() {
    cout << "proc " << _mpi_rank << ": ";  
    for (int i = 0; i < _numtets_local; i ++) {
        if (_tet_global_id[i] < 10)
            cout << " ";
        cout << _tet_global_id[i] << " ";
    }
    cout << endl;
}

void HDF5ParallelReader::read_tets() {

    hsize_t     dimsf[2];
    hsize_t     count[2];
#if defined(USE_HDF5_HSSIZE)
    hssize_t    offset[2];
#else
    hsize_t     offset[2];
#endif
    herr_t      status;
   
    hid_t dset_id = H5Dopen(_file_id, "/tetids");
    assert(dset_id >= 0);

// Create the dataspace for the dataset.
    hid_t x = H5Dget_space(dset_id);
    H5Sget_simple_extent_dims(x, dimsf, NULL);
    hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
    _numtets_global = dimsf[0];

// Set offset[2] which describes where in the file to start reading
    _numtets_local = _numtets_global / _mpi_size;
    int remainder = _numtets_global % _mpi_size;

    offset[0] = _mpi_rank * _numtets_local;
    if (_mpi_rank < remainder) {
        offset[0] += _mpi_rank;
        _numtets_local += 1;
    } else 
        offset[0] += remainder;
    offset[1] = 0;

    count[0] = _numtets_local;
    count[1] = 4; //Each tet has 4 nodes
    hid_t memspace = H5Screate_simple(2, count, NULL);
    
// Select hyperslab in the file.
    filespace = H5Dget_space(dset_id);
    assert(filespace >= 0);

    H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
    assert(filespace >= 0);

    _tetids = new idxtype[_numtets_local * 4];

// Create property list for collective dataset write.
    _plist_id = H5Pcreate(H5P_DATASET_XFER);
    H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);

    status = H5Dread(dset_id, INTMEMTYPE, memspace, filespace, _plist_id, _tetids);
    assert(status >= 0);

// store global numbers of local tets
    _tet_global_id = new int[_numtets_local];
    for (int i = 0; i < _numtets_local; i ++)
        _tet_global_id[i] = offset[0] + i;
   
// Close HDF5 stuff
    H5Dclose(dset_id);
    H5Sclose(filespace);
    H5Sclose(memspace);
    
    rInfoAll("Number of tets: %d local, %d global", _numtets_local, _numtets_global);
}


void HDF5ParallelReader::read_mats() {

    hsize_t     dimsf[2];
    hsize_t     count[2];
#if defined(USE_HDF5_HSSIZE)
    hssize_t    offset[2];
#else
    hsize_t     offset[2];
#endif
    herr_t      status;
   
    hid_t dset_id = H5Dopen(_file_id, "/matids");
    assert(dset_id >= 0);

    // Create the dataspace for the dataset.
    hid_t x = H5Dget_space(dset_id);
    H5Sget_simple_extent_dims(x, dimsf, NULL);
    hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
    // _numtets_global = dimsf[0];                 OSCAR: As many materials as tetrahedra ...

    // Set offset[2] which describes where in the file to start reading
    _numtets_local = _numtets_global / _mpi_size;
    int remainder = _numtets_global % _mpi_size;

    offset[0] = _mpi_rank * _numtets_local;
    if (_mpi_rank < remainder) {
        offset[0] += _mpi_rank;
        _numtets_local += 1;
    } else 
        offset[0] += remainder;
    offset[1] = 0;

    count[0] = _numtets_local;
    count[1] = 1; //Each tet has 1 material id
    hid_t memspace = H5Screate_simple(2, count, NULL);
    
    // Select hyperslab in the file.
    filespace = H5Dget_space(dset_id);
    assert(filespace >= 0);

    H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
    assert(filespace >= 0);

    _matids = new idxtype[_numtets_local];

    // Create property list for collective dataset write.
    _plist_id = H5Pcreate(H5P_DATASET_XFER);
    H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);

    status = H5Dread(dset_id, INTMEMTYPE, memspace, filespace, _plist_id, _matids);
    assert(status >= 0);    

    rDebug("Number of material ids: %d", _numtets_local);

    // Close HDF5 stuff
    H5Dclose(dset_id);
    H5Sclose(filespace);
    H5Sclose(memspace);
}




void HDF5ParallelReader::partition() {
// All processors will build data to call ParMetis_V3_PartMeshKway

// elmdist specifies distribution of elements
    int remainder = _numtets_global % _mpi_size;
    idxtype* elmdist = new idxtype[_mpi_size + 1];
    int off = 0;
    for (int i = 0; i < _mpi_size + 1; i ++) {
        elmdist[i] = off;
        off += _numtets_global / _mpi_size;
        if (i < remainder) 
            off += 1;
    }

// eind specifies at what index of eptr the elements that make up node i starts
    idxtype* eptr = new idxtype[_numtets_local + 1];
    for (int i = 0; i < _numtets_local + 1; i ++) {
        eptr[i] = 4 * i;
    }
    
// with elmdist and eptr as above, eind is equal to data!

// More ParMetis arguments
    int wgtflag = 0; //No weights
    int numflag = 0; //C-style numbering starts from 0

    int ncon = 1; //number of weights on each vertex
    int ncommonnodes = 3; //no of common nodes to be considered connected
    int nparts = _mpi_size; //number of subdomains in output

    float* tpwgts = new float[nparts]; //fraction of nodes to each processor
    for (int i = 0; i < nparts; i ++)
        tpwgts[i] = 1.0 / (float) nparts;

    float ubvec = 1.05; //imbalance acceptance

    int options[10];
    options[0] = 0; //1 means debug
    options[1] = 3; //level of output between 1..127
    options[2] = 15; //random seed

//Output pointers
    _tet_partition = new idxtype[_numtets_local]; //output of ParMETIS
    int edgecut; //The number of edges cut

//Call ParMETIS
    ParMETIS_V3_PartMeshKway(elmdist, eptr, _tetids, //OK data may not be idx
                             NULL, &wgtflag, //OK
                             &numflag, &ncon, //OK
                             &ncommonnodes, &nparts, //OK
                             tpwgts, &ubvec, //OK
                             options, &edgecut, //OK
                             _tet_partition, &_comm); //OK

//Some output
//     if (_mpi_rank == 0)
//      cout << "edgecut: " << edgecut << endl;

//     if (edgecut < 100) {
//        cout << "partition of processor " << _mpi_rank << ": ";
//     for (int i = 0; i < _numtets_local; i ++)
//      cout << _tet_partition[i] << " ";
//     cout << endl;
//     }

    delete[] tpwgts;
    delete[] elmdist;
    delete[] eptr;
}

void HDF5ParallelReader::sort_buffer(int* sendcnts) {

//This sorts the sendbuffer by the target proc.
//A bin/bucket sort where the sizes of the buckets is known

    int* next_pos = new int[_mpi_size];
    next_pos[0] = 0;
    for (int i = 1; i < _mpi_size; i ++) 
        next_pos[i] = next_pos[i - 1] + sendcnts[i - 1];

    int* sorted_buffer = new int[4 * _numtets_local];
    int* sorted_matbuffer = new int[_numtets_local];
    int* sorted_partition = new int[_numtets_local];
    int* sorted_tet_global_id = new int[_numtets_local];

    for (int i = 0; i < _numtets_local; i ++) {
      for (int j = 0; j < 4; j ++){
        sorted_buffer[next_pos[_tet_partition[i]] + j] = _tetids[4 * i + j];
      }
      sorted_matbuffer[next_pos[_tet_partition[i]]/4] = _matids[i];             // OSCAR: Everything you do with tets you do with mats (div 4)
        
      sorted_partition[next_pos[_tet_partition[i]] / 4] = _tet_partition[i];
      sorted_tet_global_id[next_pos[_tet_partition[i]] / 4] = _tet_global_id[i];
      next_pos[_tet_partition[i]] += 4;
    }

    delete[] _tetids;
    _tetids = sorted_buffer;
    delete[] _matids;
    _matids = sorted_matbuffer;
    delete[] _tet_global_id;
    _tet_global_id = sorted_tet_global_id;

//     if (_numtets_local < 20) {
//      cout << "sorted partition of processor " << _mpi_rank << ": ";
//      for (int i = 0; i < _numtets_local; i ++)
//          cout << sorted_partition[i] << " ";
//      cout << endl;
//     }

    delete[] next_pos;
    delete[] sorted_partition;    

}

void HDF5ParallelReader::distribute() {
/* distribute just tets, other things will be read after distribution to reduce communication
   there is no need for tet numbering as long as node/tet pairs remain intact */

//Count how many POINTS to send to each proc (sendcnts)
    int* sendcnts = new int[_mpi_size];
    for (int i = 0; i < _mpi_size; i ++)
        sendcnts[i]=0;
    for (int i = 0; i < _numtets_local; i++)
        sendcnts[_tet_partition[i]] += 4;

 //    cout << "Proc " << _mpi_rank << " sends ";
//     for (int i = 0; i < _mpi_size; i ++)
//      cout << sendcnts[i] << " ";
//     cout << " to the others" << endl;
   
//Use MPI_Alltoall to find no of points proc i will recieve from proc others (recvcnts)
    int* recvcnts = new int[_mpi_size];
    MPI_Alltoall (sendcnts, 1, MPI_INT, recvcnts, 1, MPI_INT, _comm);

 //    cout << "Proc " << _mpi_rank << " recieves ";
//     for (int i = 0; i < _mpi_size; i ++)
//      cout << recvcnts[i] << " ";
//     cout << " from the others" << endl;

//Construct send displacements (sdispls)
    int* sdispls = new int[_mpi_size];
    sdispls[0] = 0;
    for (int i = 1; i < _mpi_size; i ++)
        sdispls[i] = sdispls[i - 1] + sendcnts[i - 1];
    
//Construct recieve displacements (rdispls)
    int* rdispls = new int[_mpi_size];
    rdispls[0] = 0;
    for (int i = 1; i < _mpi_size; i ++)
        rdispls[i] = rdispls[i - 1] + recvcnts[i - 1];

//Problem: The sendbuffer is not sequential
//Must sort the sendbuffer.
    sort_buffer(sendcnts);
    //print_tet_global_id();

    int recvbuf_size = 0;
    for (int i = 0; i < _mpi_size; i ++)
        recvbuf_size += recvcnts[i];
    int* recvbuf = new int[recvbuf_size];

    
    MPI_Alltoallv (_tetids, sendcnts, sdispls, MPI_INT, recvbuf, recvcnts, rdispls, MPI_INT, _comm );
    delete[] _tetids;

    _tetids = recvbuf;
    _numtets_local = recvbuf_size / 4;
    
    // Send global materials and numbering _tet_global_id. Reuse sendcnts, recvcnts, sdispls, rdispls and recvbuf_size
    for (int i = 0; i < _mpi_size; i++) {
        sendcnts[i] /= 4;
        recvcnts[i] /= 4;
        sdispls[i] /= 4;
        rdispls[i] /= 4;
    }
    recvbuf_size /=4;
    
    recvbuf = new int[recvbuf_size];
    MPI_Alltoallv (_matids, sendcnts, sdispls, MPI_INT, recvbuf, recvcnts, rdispls, MPI_INT, _comm );
    delete[] _matids;
    _matids = recvbuf;


    int* distributed_tet_global_id = new int[recvbuf_size];
    MPI_Alltoallv (_tet_global_id, sendcnts, sdispls, MPI_INT, distributed_tet_global_id, recvcnts, rdispls, MPI_INT, _comm);
    delete[] _tet_global_id;
    _tet_global_id = distributed_tet_global_id;
    

   //  if (_mpi_rank == 0) {
//      for (int i = 0; i < _numtets_local; i ++) {
//          for (int j = 0; j < 4; j ++)
//              cout << _tetids[4 * i + j] << " ";
//          cout << endl;
//      }
//     }
    delete[] sendcnts;
    delete[] recvcnts;
    delete[] sdispls;
    delete[] rdispls;
}

void HDF5ParallelReader::generate_point_global_id(int* slab_offsets, int* slab_lengths, int slab_count) {
    _point_global_id = new int[_numpoints_local];
    int counter = 0;
    
    for(int i = 0; i < slab_count; i ++) {
        for(int j = 0; j < slab_lengths[i]; j ++) {
            _point_global_id[counter] = slab_offsets[i] + j;
            
            //cout << _point_global_id[counter] << " ";
            
            counter += 1;
            
            
        }
    }
}

void HDF5ParallelReader::read_points() {
    
    hid_t dset_id = H5Dopen(_file_id, "/coords");
    assert(dset_id >= 0);

// Create the dataspace for the dataset.
    hsize_t dimsf[2]; 
    herr_t status;
 
    hid_t x = H5Dget_space(dset_id);
    H5Sget_simple_extent_dims(x, dimsf, NULL);
//    hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
    _numpoints_global = dimsf[0];

   
//What nodes should proc i read
//Assume no of points known
    bool* read_this_node = new bool[_numpoints_global];
    _numpoints_local = 0;
    for (int i = 0; i < _numpoints_global; i ++)
        read_this_node[i] = 0;
    for (int i = 0; i < 4 * _numtets_local; i ++) {
        if (read_this_node[_tetids[i]] == 0)
            _numpoints_local += 1;
        read_this_node[_tetids[i]] = 1;
    }

//Use read_this_node[] to generate slabs
    int slab_count = 0;
    bool prev = 0;
    for (int i = 0; i < _numpoints_global; i ++) {
        if (read_this_node[i] != prev) {
            prev = read_this_node[i];
            if (read_this_node[i] == 1)
                slab_count += 1;
        }
    }

    int* slab_offsets = new int[slab_count];
    int* slab_lengths = new int[slab_count];
    int counter = 0;
    int start = 0;
    prev = 0;
//And go through read_this_node again
    for (int i = 0; i < _numpoints_global; i ++) {
        if (read_this_node[i] != prev) {
            prev = read_this_node[i];
            if (read_this_node[i] == 1) {
                start = i;
                slab_offsets[counter] = start;
                slab_lengths[counter] = _numpoints_global - start;
            } else {
                slab_lengths[counter] = i - start;
                counter ++;
            }
        }
    }
    
    //Generate global_local point index map
    generate_point_global_id(slab_offsets, slab_lengths, slab_count);

//     cout << "Proc " << _mpi_rank << " reads: ";
// //     for (int i = 0; i < _numpoints_global; i ++)
// //   cout << read_this_node[i];
// //     for (int i = 0; i < slab_count; i ++) 
// //   cout << " " << slab_offsets[i] << "-" << slab_offsets[i] + slab_lengths[i] - 1;
//     cout << " totally " << slab_count << " slabs" << endl;

 //    cout << "Number of points globally: " << _numpoints_global << endl;
//     cout << "Number of points locally: " << _numpoints_local << endl;
//     cout << "Number of slabs: " << slab_count << endl;

    hsize_t count[2];    
    count[0] = _numpoints_local;
    count[1] = 3;
    hid_t memspace = H5Screate_simple(2, count, NULL);
    assert (memspace >= 0);

//Construct hyperslab in file
    hid_t filespace = H5Dget_space(dset_id);
    assert (filespace >= 0);

#if defined(USE_HDF5_HSSIZE)
    hssize_t sstart[2];
#else
    hsize_t sstart[2];
#endif
    hsize_t sblock[2];
    hsize_t sstride[2];
    hsize_t scount[2];
    sblock[0] = 1; sblock[1] = 1;
    sstride[0] = 1; sstride[1] = 1;
    sstart[1] = 0;
    scount[1] = 3;
    for (int i = 0; i < slab_count; i ++) {
        sstart[0] = slab_offsets[i]; 
        scount[0] = slab_lengths[i];
        if (i == 0)
            status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, sstart, sstride, scount, sblock);
        else
            status = H5Sselect_hyperslab(filespace, H5S_SELECT_OR, sstart, sstride, scount, sblock);
        
        assert (status >= 0);
     }

    H5Sget_simple_extent_dims(filespace, dimsf, NULL);
  //   if (_mpi_rank == 0) {
//      cout << "Size of filespace: " << dimsf[0] << "x" << dimsf[1] << endl;
//      cout << "Slab union # elems: " << H5Sget_select_npoints(filespace) << endl;
//     }


    H5Sget_simple_extent_dims(memspace, dimsf, NULL);
  //   if (_mpi_rank == 0)
//      cout << "Size of memspace: " << dimsf[0] << "x" << dimsf[1] << endl;


//Create property list for collective dataset read.
    _plist_id = H5Pcreate(H5P_DATASET_XFER);
    H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);

    _point_coords = new double[3 * _numpoints_local];
    status = H5Dread(dset_id, REALMEMTYPE, memspace, filespace, _plist_id, _point_coords);
    assert(status >= 0);

// Close HDF5 stuff
    H5Dclose(dset_id);
    H5Sclose(filespace);
    H5Sclose(memspace);

    delete[] slab_offsets;
    delete[] slab_lengths;
    delete [] read_this_node;
}

void HDF5ParallelReader::read_boundary() {

//So first read in all boundary points
    hsize_t     dimsf[2];                 /* dataset dimensions */
    herr_t      status;
   
    hid_t dset_id = H5Dopen(_file_id, "/surface");
    assert(dset_id >= 0);

// Create the dataspace for the dataset.
    hid_t x = H5Dget_space(dset_id);
    H5Sget_simple_extent_dims(x, dimsf, NULL);
    hid_t filespace = H5Screate_simple(2, dimsf, NULL); 
    _numbfaces_global = dimsf[0];

    int* allbfaces = new idxtype[_numbfaces_global * 4];

// Create property list for collective dataset write.
    _plist_id = H5Pcreate(H5P_DATASET_XFER);
    H5Pset_dxpl_mpio(_plist_id, H5FD_MPIO_COLLECTIVE);

    status = H5Dread(dset_id, INTMEMTYPE, H5S_ALL, filespace, _plist_id, allbfaces);
    assert(status >= 0);
        
// Store local boundary faces, discard others.
   _numbfaces_local = 0;
   int nof_sym = 0; //Count number of symmetry planes.

    for (int i = 0; i < _numbfaces_global; i ++) {
        try {
            _builder->set_bc(allbfaces[4*i+1], allbfaces[4*i+2], allbfaces[4*i+3], allbfaces[4*i]);
            _numbfaces_local += 1;
        } catch (std::runtime_error err) {}//Face was not local
        
        if (allbfaces[4 * i] > nof_sym)
            nof_sym = allbfaces[4 * i];
    }
    
    //cout << "Proc " << _mpi_rank <<": " << _numbfaces_local << " faces, " << nof_sym << " symmetry planes" << endl;
    
    _builder->finalize_bc(nof_sym);

// Close HDF5 stuff
    H5Dclose(dset_id);
    H5Sclose(filespace);
    
    delete[] allbfaces;
}

void HDF5ParallelReader::export_dist(char* file) {

    std::ofstream of;
    int* recvcnts;

//Proc 0 must read all points from file then output them
    if (_mpi_rank == 0) {
        of.open (file);
        assert (of.is_open());
        of.precision (6);

        of << "# vtk DataFile Version 2.0" << endl;                      // first line of a vtk file
        of << "generated using hdf5_parallel_read::export_dist" << endl;   // header
        of << "ASCII" << endl << endl;                                   // file format
        of << "DATASET UNSTRUCTURED_GRID" << endl;
        of << "POINTS " << _numpoints_global << " float" << endl;
        of << scientific; //????????

        hid_t dataset = H5Dopen(_file_id, "/coords");
        assert (dataset >= 0);

        double* nodes = new double[3 * _numpoints_global];
        herr_t status = H5Dread(dataset, REALMEMTYPE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nodes);
        assert (status >= 0);

        for (int i = 0; i < _numpoints_global; i ++)
            of << nodes[3 * i] << " " << nodes[3 * i + 1] << " " << nodes[3 * i + 2] << endl;
        of << endl;

        H5Dclose (dataset);
        delete[] nodes;

//Read all tets from file and output them
        of << "CELLS " << _numtets_global << " " << 5 * _numtets_global << endl;
    }

    recvcnts = new int[_mpi_size];
    MPI_Gather (&_numtets_local, 1, MPI_INT, 
                recvcnts, 1, MPI_INT, 
                0, _comm);
    int* displs = new int[_mpi_size];
    displs[0] = 0;
    recvcnts[0] *= 4;
    for (int i = 1; i < _mpi_size; i ++) {
        recvcnts[i] *= 4;
        displs[i] = displs[i - 1] + recvcnts[i - 1];
    }
//     for (int i = 0; i < _mpi_size; i ++) {
//      if (_mpi_rank == 0)
//          cout << recvcnts[i] << " "<< displs[i] << endl;
//     }

    int* total_partition = new int[4 * _numtets_global];
    MPI_Gatherv (_tetids, 4 * _numtets_local, MPI_INT, 
                 total_partition, recvcnts, displs,
                 MPI_INT, 0, _comm);

    delete[] displs;

    if (_mpi_rank == 0) {

        for (int i = 0; i < _numtets_global; i ++) {
            of << "4 ";
            for (int j = 0; j < 4; j ++)
                of << total_partition[4 * i + j] << " ";
            of << endl;
        }
        of << endl;

        of << "CELL_TYPES " << _numtets_global << endl;
        for (int i = 0; i < _numtets_global; i ++)
            of << 10 << endl;
        of << endl;
        of << endl;
    }

// 

// Proc 0 writes colormap
    if (_mpi_rank == 0) {

        of << "CELL_DATA " << _numtets_global << endl << endl;
        of << "SCALARS process_no int 1" << endl;
        of << "LOOKUP_TABLE default" << endl;

        for (int i = 0; i < _mpi_size; i ++)
            for (int j = 0; j < recvcnts[i] / 4; j ++)
                of << i << endl;
    }


    
}
//
void HDF5ParallelReader::open_file() {
// Set up file access property list with parallel I/O access
    _plist_id = H5Pcreate(H5P_FILE_ACCESS); //Property list identifier
    H5Pset_fapl_mpio(_plist_id, _comm, _info);   
// Create a new file collectively and release property list identifier.
    _file_id = H5Fopen(_inputFile, H5F_ACC_RDONLY, _plist_id);
    assert(_file_id >= 0);
    H5Pclose(_plist_id);
}

void HDF5ParallelReader::close_file() {
    H5Pclose(_plist_id);
    H5Fclose(_file_id);
}
    
void HDF5ParallelReader::read() {
    open_file();
    
    read_tets();

    read_mats();
    
    partition();
    
    distribute();
    
    read_points();
    
    _builder->init_coord(_numpoints_local);
    for (int i = 0; i < _numpoints_local; i ++)
        _builder->set_coord(_point_global_id[i], _point_coords[3 * i], _point_coords[3 * i + 1], _point_coords[3 * i + 2]);
        
    _builder->finalize_coord();
        
    _builder->init_tet(_numtets_local);       
    for (int i = 0; i < _numtets_local; i ++) {
      _builder->set_tet(i, _tetids[4 * i], _tetids[4 * i + 1], _tetids[4 * i + 2], _tetids[4 * i + 3], _matids[i]);
    }

    _builder->finalize_tet();
    
    read_boundary();
    
    close_file();
}

int HDF5ParallelReader::main (int argc, char **argv) {

#ifdef __GNUC__
    std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
#endif

    MPI_Init(&argc, &argv);    
    
    char* inputFile;
    
    if (argc != 2) {
        if (_mpi_rank == 0) 
            cout << "Usage: hdf5_parallel_read inputfile.h5" << endl;
        return 0;
    } else {
        inputFile = argv[1];
    }
    
    
    MPI_Comm comm = MPI_COMM_WORLD;
    
    TetMeshBuilder* builder = new TetMeshBuilder(comm);
    HDF5ParallelReader* reader = new HDF5ParallelReader(inputFile, builder, comm);
    
    double time;
    

    

    //time = seconds();
    open_file(); 
    //cout << "Proc: " << _mpi_rank << " open_file: " << seconds() - time << endl;

    time = seconds();
    read_tets();
    read_mats();
    //cout << "Proc: " << _mpi_rank << " read_tets: " << seconds() - time << endl;
    //print_tet_global_id();  
    
    time = seconds();
    partition(); 
    //cout << "Proc: " << _mpi_rank << " partition: " << seconds() - time << endl;

    time = seconds();
    reader->distribute();
    //cout << "Proc: " << _mpi_rank << " distribute: " << seconds() - time << endl;
    //print_tet_global_id();  
 
    time = seconds();
    reader->read_points();
    //cout << "Proc: " << _mpi_rank << " read_points: " << seconds() - time << endl;

    
    _builder->init_coord(_numpoints_local);
    
     for (int i = 0; i < _numpoints_local; i ++)
         _builder->set_coord(_point_global_id[i], _point_coords[3 * i], _point_coords[3 * i + 1], _point_coords[3 * i + 2]);
     
    _builder->finalize_coord();
    //
    
    _builder->init_tet(_numtets_local);       
        
    for (int i = 0; i < _numtets_local; i ++) {
        //cout << _mpi_rank << ":" << i;
        //cout << ": (" << _tetids[4 * i] << "," << _tetids[4 * i + 1] << "," << _tetids[4 * i + 2] << "," << _tetids[4 * i + 3] << ") " << endl;
      _builder->set_tet(i, _tetids[4 * i], _tetids[4 * i + 1], _tetids[4 * i + 2], _tetids[4 * i + 3], _matids[i]);
    }

    _builder->finalize_tet();
    
     time = seconds();
    read_boundary();
    //cout << "Proc: " << _mpi_rank << " read_boundary: " << seconds() - time << endl;

    time = seconds();
    export_dist("out.vtk");
    //cout << "Proc: " << _mpi_rank << " export_dist: " << seconds() - time << endl;
    
    
    
    
//Close/release resources.
    close_file();
    

    mesh::ParallelTetMesh* mesh = dynamic_cast<mesh::ParallelTetMesh*> (builder->get_mesh());
    
    
    
    delete _builder;
    
    
    delete mesh;
    
    delete _tetids;
    delete _tet_partition;

    MPI_Finalize();
    
    return 0;
}

---