src/h5_tools.cpp

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <ostream>
#include "mpi.h"
#include "hdf5.h"
#include <Epetra_Comm.h>
#include <Epetra_MpiComm.h>
#include <Epetra_Map.h>
#include <Epetra_Vector.h>
#include <Epetra_MultiVector.h>
#include <Epetra_SerialDenseVector.h>
#include <Epetra_Import.h>
#include <Epetra_Export.h>
#include <Teuchos_ParameterList.hpp>
#include "rlog/rlog.h"
#include "colarray/matrix.h"
#include "colarray/vector.h"
#include "tetmesh.h"
#include "femaxmesh.h"
#include "nedelecelement.h"
#include "nedelecelement1.h"
#include "nedelecelement2.h"
#include "nedelecmesh.h"
#include "paralleltetmesh.h"
#include "vector3.h"
#include "vector4.h"
#include "fmxxtopology.h"
#include "fmxxphysicomath.h"
#include "nonsciconst.h"
#include "dbgsupport.h"

Epetra_Map* h5_write_eigenfield_retrieve_DoF_robust(const Epetra_MultiVector& Q,  const Epetra_Comm& comm);

Epetra_Map* h5_write_eigenfeield_retrieve_DoF_efficient(const std::string& file_name,  const Epetra_Comm& comm,
                                                                                                        mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                                                                                                        const Epetra_MultiVector& Q, double* lambda);

void h5_write_eigenfield_write_hdf5e(const std::string& file_name,
                                                                const Epetra_Comm& comm,
                                                                int lntet,
                                                                int gntet,
                                                                int mthmode,
                                                                double* ev);

void h5_write_eigenfield_write_hdf5h(const std::string& file_name,
                                                                const Epetra_Comm& comm,
                                                                int lntet,
                                                                int gntet,
                                                                int mthmode,
                                                                double* hv);

void h5_write_eigenfield_writehdf5sloc(const std::string& file_name,
                                                                const Epetra_Comm& comm,
                                                                int ntet,
                                                                int gntet,
                                                                double* sloc);

void h5_compute_eigenquality(const std::string& file_name,  const Epetra_Comm& comm,
                                                                                mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                                                                                const Epetra_MultiVector& Q, double* lambda, double* &quality);


void h5_cartesian_sampling(const std::string& file_name,  const Epetra_Comm& comm,
                            mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                            const Epetra_MultiVector& Q, double* lambda,
                            double xmin, double xmax, int nx,
                            double ymin, double ymax, int ny,
                            double zmin, double zmax, int nz);

void h5_cartesian_sampling_write_hdf5(const Epetra_Comm& comm, const std::string& file_name, unsigned int nmode,
                    std::vector<std::vector<std::vector<unsigned int> > >& exists,
                    std::vector<std::vector<std::vector<std::vector<mesh::Vector3> > > >& efield,
                    std::vector<std::vector<std::vector<std::vector<mesh::Vector3> > > >& hfield,
                    double xmin, double xmax, const int nx,
                    double ymin, double ymax, const int ny,
                    double zmin, double zmax, const int nz);






void h5_create_empty_file(const std::string& file_name) {
    hid_t file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    if (file_id < 0)
        rError("Unable to create hdf5 file \"%s\".", file_name.c_str());
    H5Fclose(file_id);
}

static void h5_write_param_list_recursive(const Teuchos::ParameterList& params, hid_t group_id) {
    Teuchos::ParameterList::ConstIterator it = params.begin();
    for (; it != params.end(); ++ it) {
        std::string key(params.name(it));
        if (params.isSublist(key)) {
            //
            // Sublist
            //
            
            // Create subgroup for sublist.
            hid_t child_group_id = H5Gcreate(group_id, key.c_str(), 0);
            h5_write_param_list_recursive(params.sublist(key), child_group_id);
            H5Gclose(child_group_id);
        } else {
            //
            // Regular parameter
            //
            
            // Create dataspace/dataset.
            herr_t status;
            hsize_t one = 1;
            hid_t dataspace_id, dataset_id;

            // Write the dataset.
            if (params.isType<int>(key)) {
                int value = params.get<int>(key);
                dataspace_id = H5Screate_simple(1, &one, NULL);
                dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_INT, dataspace_id, H5P_DEFAULT);
                status = H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, 
                                  H5P_DEFAULT, &value);
                status = H5Dclose(dataset_id);
                status = H5Sclose(dataspace_id);
            } else if (params.isType<double>(key)) {
                double value = params.get<double>(key);
                dataspace_id = H5Screate_simple(1, &one, NULL);
                dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT);
                status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, 
                                  H5P_DEFAULT, &value);
                status = H5Dclose(dataset_id);
                status = H5Sclose(dataspace_id);
            } else if (params.isType<std::string>(key)) {
                std::string value = params.get<std::string>(key);
                hsize_t len = value.size() + 1;
                dataspace_id = H5Screate_simple(1, &len, NULL);
                dataset_id = H5Dcreate(group_id, key.c_str(), H5T_C_S1, dataspace_id, H5P_DEFAULT);
                status = H5Dwrite(dataset_id, H5T_C_S1, H5S_ALL, H5S_ALL, 
                                  H5P_DEFAULT, value.c_str());
                status = H5Dclose(dataset_id);
                status = H5Sclose(dataspace_id);
            } else if (params.isType<bool>(key)) {
                // Use H5T_NATIVE_USHORT to store a bool value
                unsigned short value = params.get<bool>(key) ? 1 : 0;
                dataspace_id = H5Screate_simple(1, &one, NULL);
                dataset_id = H5Dcreate(group_id, key.c_str(), H5T_NATIVE_USHORT, dataspace_id, H5P_DEFAULT);
                status = H5Dwrite(dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, 
                                  H5P_DEFAULT, &value);
                status = H5Dclose(dataset_id);
                status = H5Sclose(dataspace_id);
            } else {
                rError("Unable to export parameter \"%s\" to hdf5 file: Unsupported type.", key.c_str());        
            }
        }
    }
}

void h5_write_param_list(const std::string& file_name, const Teuchos::ParameterList& params) {
    // Open hdf5 file
    hid_t       file_id, group_id;  /* identifiers */
    herr_t      status;
    // file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);

    // Create group "parameters" in the root group using absolute name.
    group_id = H5Gcreate(file_id, "/parameters", 0);

    // Iterate through parameter list 
    h5_write_param_list_recursive(params, group_id);
    
    // Finalize hdf5 file
    status = H5Gclose(group_id);
    status = H5Fclose(file_id);
 }

static herr_t f_operator(hid_t loc_id, const char *name, void *opdata) {
    H5G_stat_t statbuf;
    hid_t dataset_id, space_id, type_id;
    Teuchos::ParameterList* sublist;
    Teuchos::ParameterList* params = (Teuchos::ParameterList*)opdata;
    /*
     * Get type of the object and display its name and type.
     * The name of the object is passed to this function by 
     * the Library. Some magic :-)
     */
    H5Gget_objinfo(loc_id, name, 0, &statbuf);
    switch (statbuf.type) {
    case H5G_GROUP:
        sublist = &params->sublist(name);
        H5Giterate(loc_id, name , NULL, f_operator, sublist);
        break;
    case H5G_DATASET:
        hsize_t len;
        dataset_id = H5Dopen(loc_id, name);
        space_id = H5Dget_space(dataset_id);
        if (H5Sget_simple_extent_ndims(space_id) != 1)
            rError("Dimensionality of parameters must be 1.");
        H5Sget_simple_extent_dims(space_id, &len, NULL);
        type_id = H5Dget_type(dataset_id);
        if (H5Tequal(type_id, H5T_NATIVE_DOUBLE) > 0) {
            double value;
            rAssert(len == 1);
            H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
            params->set(name, value);
        } else if (H5Tequal(type_id, H5T_NATIVE_INT) > 0) {
            int value;
            rAssert(len == 1);
            H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
            params->set(name, value);
        } else if (H5Tequal(type_id, H5T_C_S1) > 0) {
            char* buf = new char[len];
            H5Dread(dataset_id, H5T_C_S1, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
            params->set(name, std::string(buf));
            delete[] buf;
        } else if (H5Tequal(type_id, H5T_NATIVE_USHORT) > 0) {
            unsigned short value;
            rAssert(len == 1);
            H5Dread(dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &value);
            params->set(name, value != 0 ? true : false);
        } else {
            rError("Unsupported datatype for parameter \"%s\".", name);   
        }
        H5Tclose(type_id);
        H5Sclose(space_id);  
        H5Dclose(dataset_id);  
        break;
    default:
        rError("Unsupported type for parameter \"%s\".", name);
    }
    return 0;
}

void h5_read_param_list(const std::string& file_name, Teuchos::ParameterList& params) {
    // Open hdf5 file
    hid_t       file_id, group_id;  /* identifiers */
    herr_t      status;
    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);

    // Create group "parameters" in the root group using absolute name.
    group_id = H5Gopen(file_id, "/parameters");

    // Iterate through parameter list 
    status = H5Giterate(group_id, "/parameters" , NULL, f_operator, &params);
    
    // Finalize hdf5 file
    status = H5Gclose(group_id);
    status = H5Fclose(file_id);
 }

void h5_write_eigenmodes(const std::string& file_name, const Epetra_MultiVector& Q, const Epetra_SerialDenseVector& lambda)
{
        hid_t       file_id, group_id, dset_id;       /* file and dataset identifiers */
        hid_t       filespace, memspace;              /* file and memory dataspace identifiers */
        hid_t       plist_id;                         /* property list identifier */
        herr_t      status;

        // Set up file access property list with parallel I/O access
        plist_id = H5Pcreate(H5P_FILE_ACCESS);
    MPI_Comm mpi_comm  = dynamic_cast<const Epetra_MpiComm&>(Q.Comm()).GetMpiComm();
    MPI_Info mpi_info  = MPI_INFO_NULL;
    H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

    // Open file collectively and release property list identifier.
    // file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
    H5Pclose(plist_id);

    // Create a group named /modes in the file.
    group_id = H5Gcreate(file_id, "/modes", 0);

    //
    // Eigenvectors
    //

    // Redistribute eigenvectors to force linear map
    Epetra_Map target_map(Q.GlobalLength(), 0, Q.Comm());
    Epetra_MultiVector Q0(target_map, Q.NumVectors());
    Q0.Export(Q, Epetra_Export(Q.Map(), target_map), Add);
    assert(Q0.Map().LinearMap());

    // Enforce defined direction for all eigenvectors
    for (int k = 0; k < Q0.NumVectors(); ++ k) {
        Epetra_Vector& q0 = *Q0(k);
        int sign = 1;
        // Proc 0 determines whether vector will be negated
        if (Q0.Comm().MyPID() == 0) {
            int i;
            double val;
            for (i = 0; i < Q0.MyLength() && (val = q0[i]) == 0.0; ++ i);
            if (i < Q0.MyLength() && val < 0.0)
                sign = -1;
        }
        // Brodcast sign to all procs and negate if necessary
        Q0.Comm().Broadcast(&sign, 1, 0);
        if (sign == -1) {
            for (int i = 0; i < Q0.MyLength(); ++ i)
                q0[i] = -q0[i];
        }
    } // for k

#if 0
    // For debugging: fill eigenvectors
    for (int j = 0; j < Q0.NumVectors(); ++ j)
        for (int i = 0; i < Q0.MyLength(); ++ i)
            Q0.Values()[j*Q0.MyLength() + i] = 10000*(j+1) + Q0.Map().GID(i);
    // cout << *(Q0(1)) << "\n";
#endif

    // Create the dataspace for the dataset.
    hsize_t q_dimsf[] = {Q0.NumVectors(), Q0.GlobalLength()};                 /* dataset dimensions */
    filespace = H5Screate_simple(2, q_dimsf, NULL);

    // Create the dataset with default properties and close filespace.
    dset_id = H5Dcreate(group_id, "eigenvectors", H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT);

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

    // Select hyperslab in the file.
    hsize_t offset[] = {0, Q0.Map().GID(0)};
    hsize_t stride[] = {1, 1};
    hsize_t count[] = {Q0.NumVectors(), 1};
    hsize_t block[] = {1, Q0.MyLength()};
    filespace = H5Dget_space(dset_id);
    H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, stride, count, block);

    // Each process defines dataset in memory and writes it to the hyperslab in the file.
    hsize_t dimsm[] = {Q0.NumVectors()*Q0.MyLength()};
    memspace = H5Screate_simple(1, dimsm, NULL);

    // Write hyperslab
    status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, Q0.Values());
    H5Sclose(memspace);
    H5Sclose(filespace);
    H5Dclose(dset_id);

    //
    // Eigenvalues
    //
    hsize_t lambda_dimsf = lambda.Length();
    filespace = H5Screate_simple(1, &lambda_dimsf, NULL);
    dset_id = H5Dcreate(group_id, "eigenvalues", H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT);
    status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, plist_id, lambda.Values());
    H5Dclose(dset_id);
    H5Sclose(filespace);

    // Close/release resources.
    H5Pclose(plist_id);
    H5Gclose(group_id);
    H5Fclose(file_id);
}

void h5_read_eigenmodes(const std::string& file_name, Epetra_MultiVector** Q, Epetra_SerialDenseVector** lambda) {
}

void h5_read_eigenmodes(const std::string& file_name, 
                        colarray::Matrix<double>& Q, 
                        colarray::ColumnVector<double>& lambda) 
{
    hid_t file_id, group_id, dataset_id, space_id;
    hsize_t num_eigenpairs, dims[2];
    herr_t status;
    
    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
        group_id = H5Gopen(file_id, "/modes");
    
        dataset_id = H5Dopen(group_id, "eigenvalues");
        space_id = H5Dget_space(dataset_id);
        if (H5Sget_simple_extent_ndims(space_id) != 1)
        {
                rError("Dimensionality of modes/eigenvalues must be 1.");
        }
        
        status = H5Sget_simple_extent_dims(space_id, &num_eigenpairs, NULL);
        lambda = colarray::ColumnVector<double>(num_eigenpairs);
        status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
                     lambda._v);
    H5Sclose(space_id);  
    H5Dclose(dataset_id);  

    dataset_id = H5Dopen(group_id, "eigenvectors");
    space_id = H5Dget_space(dataset_id);
    if (H5Sget_simple_extent_ndims(space_id) != 2)
        rError("Dimensionality of modes/eigenvectors must be 2.");
    status = H5Sget_simple_extent_dims(space_id, dims, NULL);
    if (dims[0] != num_eigenpairs)
        rError("Number of eigenvalues (%d) does not match number of eigenvectors (%d).", 
               static_cast<int>(num_eigenpairs), 
               static_cast<int>(dims[0]));
    Q = colarray::Matrix<double>(dims[1], dims[0]);
    status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
                     Q._v);
    H5Sclose(space_id);  
    H5Dclose(dataset_id);  

    H5Gclose(group_id);
    H5Fclose(file_id);
}

 void h5_read_eigenmesh(const std::string& file_name)
 {
        ;
 }

void h5_write_eigenmesh(const std::string& file_name,const Epetra_Comm& comm, mesh::TetMesh* tetmesh)
{
        int i,j,k; 
        int res=0; 
        id_t gntet=0; 
        id_t ntet=0; 
        id_t gpointid; 
        int ncol=NCORNERTET+2; 
        mesh::ParallelTetMesh* ptetmesh=NULL;
        ptetmesh=dynamic_cast<mesh::ParallelTetMesh*>(tetmesh);

        hid_t       file_id, group_id, dset_id;       /* file and dataset identifiers */
        hid_t       filespace, memspace;              /* file and memory dataspace identifiers */
        hid_t       plist_id;                         /* property list identifier */
        herr_t      status;

        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
        
        int mpirank=0; 
        int mpisize=0; 
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " reached " << endl;
#endif 
        plist_id = H5Pcreate(H5P_FILE_ACCESS);
        H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
    H5Pclose(plist_id);

    group_id = H5Gcreate(file_id, H5_GROUP_FEMAX_EIGENMESH.c_str(), 0);

        ntet=tetmesh->get_nof_tets();
         
        /* retrieve global number of tetrahedra, required for allocating dimensions of dataset */
        gntet=ptetmesh->get_num_global_tets();
         
        hsize_t dimsf[]={gntet,ncol};
        filespace = H5Screate_simple(DIM2D, dimsf, NULL);

        dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_TETVID.c_str(), H5T_NATIVE_UINT, filespace, H5P_DEFAULT);
        H5Sclose(filespace);

        hsize_t count[2];
        hsize_t offset[2];
        
        count[0]=tetmesh->get_nof_tets();
        count[1]=dimsf[1]; 
        id_t* ntetlist=new id_t[mpisize]; 
        res=MPI_Allgather(&ntet,1,MPI_UNSIGNED, ntetlist, 1, MPI_UNSIGNED, mpi_comm);

        offset[0]=0;
        for(i=0;i<mpirank;i++)
        {
                offset[0]+=ntetlist[i];
        }

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " offset[0] " << offset[0] << endl;
 #endif 
        offset[1]=0; /* no offset w.r.t columns */
        memspace=H5Screate_simple(DIM2D,count,NULL);

        filespace=H5Dget_space(dset_id);
        H5Sselect_hyperslab(filespace, H5S_SELECT_SET,offset,NULL,count,NULL);
        
        unsigned int *tvertexid=new unsigned int [ntet*ncol];
        for(i=0;i<ntet*ncol;i++)
        {
                tvertexid[i]=0;
        }
        
        mesh::TetMesh::tet_iterator tit;
        i=0;
        
        for(tit=tetmesh->tet_begin(); tit != tetmesh->tet_end(); tit++ )
        {
                for(j=0; j<NCORNERTET;j++)
                        {
                                gpointid=ptetmesh->point_gid(tit->get_corner_id(j));
                                
                                tvertexid[i*ncol+j]=gpointid;
                        }
                        
                        tvertexid[i*ncol+NCORNERTET]=mpirank;
                        tvertexid[i*ncol+NCORNERTET+1]=mpisize;
 
//#ifdef __DEBUG__VERBOSE__
//                              cout << "global vertex id: " << tvertexid[i][0] << " : " <<tvertexid[i][1] << " : " << tvertexid[i][2] << " : " << tvertexid[i][3] << endl;
// #endif /**  __DEBUG__VERBOSE__ */
                        
                        i++;
        }
        
        plist_id=H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
        
#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " reached " << endl;
 #endif 
        status=H5Dwrite(dset_id, H5T_NATIVE_UINT, memspace, filespace,plist_id, tvertexid);

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " reached " << endl;
#endif 
        delete []tvertexid;
        delete []ntetlist;

        H5Dclose(dset_id);
        H5Sclose(filespace);
        H5Pclose(plist_id);
        H5Gclose(group_id);
        H5Fclose(file_id);      
 }
 
 
 
void h5_write_eigenpoint(const std::string& file_name,  const Epetra_Comm& comm, mesh::TetMesh* tetmesh)
{
        int i,j,k; 
        int res=0; 
        id_t gntet=0; 
        id_t ntet=0; 
        int ncol=NCORNERTET+2; 
        mesh::ParallelTetMesh* ptetmesh=NULL;
        ptetmesh=dynamic_cast<mesh::ParallelTetMesh*>(tetmesh);

        hid_t           file_id, group_id, dset_id;       /* file and dataset identifiers */
        hid_t           filespace, memspace;              /* file and memory dataspace identifiers */
        hid_t           plist_id;                         /* property list identifier */
        herr_t  status;

        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
        
        int mpirank=0; 
        int mpisize=0; 
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);

        plist_id = H5Pcreate(H5P_FILE_ACCESS);
        H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

  file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
  H5Pclose(plist_id);

  group_id = H5Gopen(file_id, H5_GROUP_FEMAX_EIGENMESH.c_str());
  
        if(group_id < 0) /* there is no such group yet, create one */
        {
                group_id = H5Gcreate(file_id, H5_GROUP_FEMAX_EIGENMESH.c_str(), 0);
        }

        id_t gnpoint=0; 
        gnpoint=ptetmesh->get_num_global_points();
 
        hsize_t dimsf[]={gnpoint,NREALCOORD3D};
        filespace = H5Screate_simple(DIM2D, dimsf, NULL);

        dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_POINTS.c_str(), H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT);
        H5Sclose(filespace);

        hsize_t count[2];
        hsize_t offset[2];
        
        id_t npoint=0;
        id_t npowned=0; 
        mesh::TetMesh::point_iterator pit;
        id_t gpointid=0;
        npoint=tetmesh->get_nof_points();
        bool powned=false;
        double coords3d[NREALCOORD3D];

        npowned=0;                      
        for(pit=tetmesh->point_begin(); pit != tetmesh->point_end(); pit++ )
        {
                gpointid=ptetmesh->point_gid(pit->get_id());
                powned=ptetmesh-> is_owned_point(gpointid);
                if(powned == true)
                {
                        npowned++;
                }
        }
        double* xfercoords=new double[npowned*NREALCOORD3D];
        for(i=0;i<npowned*NREALCOORD3D;i++){xfercoords[i]=0.0; }

        count[0]=npowned;
        count[1]=NREALCOORD3D;
        memspace=H5Screate_simple(DIM2D,count,NULL);
        filespace=H5Dget_space(dset_id);
        
        hsize_t scount[2];
        scount[0]=1;
        scount[1]=NREALCOORD3D;
        
        
        i=0;
        npowned=0;
        
        for(pit=tetmesh->point_begin(); pit != tetmesh->point_end(); pit++ )
        {
                gpointid=ptetmesh->point_gid(pit->get_id());
                powned=ptetmesh-> is_owned_point(gpointid);
                
                if(powned == true)
                {
                        offset[0]=gpointid;
                        offset[1]=0;

                        if(npowned == 0)
                        {
                                H5Sselect_hyperslab(filespace, H5S_SELECT_SET,offset,NULL,scount,NULL);
                        }
                        else
                        {
                                H5Sselect_hyperslab(filespace, H5S_SELECT_OR,offset,NULL,scount,NULL);
                        }
                        for(j=0;j<NREALCOORD3D;j++)
                        {
                                xfercoords[i*NREALCOORD3D+j]=pit->get_coord()[j];
                        }

                        npowned++; 
                        i++;  
                }
        }
        
#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : h5_tools->h5_write_eigenpoints : process : " << mpirank << " with " << npowned << " points " << endl;
#endif 
        plist_id=H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
        status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace,plist_id, xfercoords);

        delete []xfercoords;
        
        H5Pclose(plist_id);
        H5Sclose(filespace);
        H5Dclose(dset_id);
        H5Gclose(group_id);
        H5Fclose(file_id);
}
 
void h5_write_eigenfield(const std::string& file_name,  const Epetra_Comm& comm,
                                                        mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                                                        const Epetra_MultiVector& Q, double* lambda)
{
        int i,j,k,m; /* standard loop indices */
        int res; /* return code */
        int ntet; /* number of tetrahedra local to this process */
        int gntet; /* global number of tetrahedra */
        int nmode; /* store number of calulated eigenmodes */
  int mpirank=0; 
  int mpisize=0; 
  hid_t       file_id, group_id, dset_id, plist_id;       /* file and dataset identifiers */

        mesh::ParallelTetMesh* ptetmesh=dynamic_cast<mesh::ParallelTetMesh*>(tetmesh);

//      NedelecMesh nedelecmesh=femaxmesh->get_nedelec_mesh();
        NedelecElement* nedelecelem=nedelecmesh.get_element(0); 
        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
 
  res=MPI_Comm_rank(mpi_comm, &mpirank);
  res=MPI_Comm_size(mpi_comm, &mpisize);


        int order=nedelecelem->get_order(); 
        int ndofelem=nedelecelem->get_nof_dofs(); 
        /* compute an Epetra_Map */
//      Epetra_Map* targetmap=h5_write_eigenfeield_retrieve_DoF_efficient(file_name,comm, tetmesh, nedelecmesh, Q,lambda); /* the efficient way */
        Epetra_Map* targetmap=h5_write_eigenfield_retrieve_DoF_robust(Q,comm); /* the robust way */
                                                                                                                
        Epetra_MultiVector Q0(*targetmap, Q.NumVectors());
        
        Epetra_Import importer(*targetmap,Q.Map());
//      Epetra_Export exporter(Q.Map(),*targetmap);
        
        Q0.Import(Q,importer,Insert);
//      Q0.Export(Q, exporter,Add);
                 

        /* evaluate the modal fields by looping over all Epetra_Multivector components */
        mesh::TetMesh::tet_iterator tit;
        mesh::Vector3* efield=NULL;
        mesh::Vector3* hfield=NULL;

        const int* tldofid=NULL;
        colarray::Vector<double> dofs(nedelecelem->get_nof_dofs());
        int mapped_id=0;

        mesh::Vector3 samploc(0,0,0);
        mesh::Tet* ttet=NULL;

        ntet=tetmesh->get_nof_tets(); 
        gntet=ptetmesh->get_num_global_tets(); /* retrieve global number of tetrahedra, required for allocating dimensions of dataset */

        /* calculate sampling locations */
        double* sloc=new double[NREALCOORD3D*ntet];
        
        k=0; /* used as index over tetrahedral sampling locations */
        for(tit=tetmesh->tet_begin(); tit != tetmesh->tet_end(); tit++ )
        {
                samploc=tit->get_center_of_gravity();
                        
                sloc[k*NREALCOORD3D]=samploc.x;
                sloc[k*NREALCOORD3D+1]=samploc.y;
                sloc[k*NREALCOORD3D+2]=samploc.z;
                
                k++;
        }       

        efield=new mesh::Vector3[ntet]; 
        hfield=new mesh::Vector3[ntet]; 
        for(i=0;i<ntet;i++)
        {
                /* initialize the electric field */
                efield[i].x=0.0;
                efield[i].y=0.0;
                efield[i].z=0.0;
                
                /* initialize the magnetic field */
                hfield[i].x=0.0;
                hfield[i].y=0.0;
                hfield[i].z=0.0;        
        }

        nmode=Q.NumVectors();
        int nev=NREALCOORD3D*ntet;
        double* ev=new double[nev]; 
        double* hv=new double[nev]; 
        for(i=0;i<nev;i++)
        {
                ev[i]=0.0;
                hv[i]=0.0;              
        }

        int tlid=0;

        for(m=0;m<nmode;m++)
        {
                Epetra_Vector mthmode(Copy, Q0,m);
                
                k=0; 
                for(tit=tetmesh->tet_begin(); tit != tetmesh->tet_end(); tit++ )
                {
                        samploc=tit->get_center_of_gravity();
                        
                        sloc[k*NREALCOORD3D]=samploc.x;
                        sloc[k*NREALCOORD3D+1]=samploc.y;
                        sloc[k*NREALCOORD3D+2]=samploc.z;
                
                        nedelecelem=nedelecmesh.get_element(tit->get_id());
                
                        ttet=&(*tit); 
                        tldofid=nedelecelem->get_dof_ids(tetmesh,ttet);
                
                        for (i=0; i < nedelecelem->get_nof_dofs(); i ++) 
                        {
                                mapped_id=nedelecmesh.map_dof(tldofid[i]);
                        
                                if (mapped_id == mesh::ID_NONE)
                                {
                                        dofs(i) = 0.0;
                                }
                                else
                                {                               
                                        tlid=targetmap->LID(mapped_id);
                                        dofs(i) = mthmode[tlid];
                                }
      }
        
                        efield[k].x=nedelecelem->eval(ttet, samploc, dofs).x;
                        efield[k].y=nedelecelem->eval(ttet, samploc, dofs).y;
                        efield[k].z=nedelecelem->eval(ttet, samploc, dofs).z;

                        hfield[k].x=nedelecelem->eval_curl(ttet, samploc, dofs).x;
                        hfield[k].y=nedelecelem->eval_curl(ttet, samploc, dofs).y;
                        hfield[k].z=nedelecelem->eval_curl(ttet, samploc, dofs).z;
                        
                        k++;
                }
        
                i=0;
                for(tit=tetmesh->tet_begin(); tit != tetmesh->tet_end(); tit++ )
                {
                        for(j=0;j<NREALCOORD3D;j++)
                        {
                                ev[i*NREALCOORD3D+j]=efield[i][j];
                                hv[i*NREALCOORD3D+j]=hfield[i][j];
                        }
                        i++;
                }

                h5_write_eigenfield_write_hdf5e(file_name, comm, ntet, gntet, m, ev);
                h5_write_eigenfield_write_hdf5h(file_name, comm, ntet, gntet, m, hv);
        }

        h5_write_eigenfield_writehdf5sloc(file_name, comm, ntet, gntet, sloc);

        delete []efield;
        delete []hfield;
        delete []ev;
        delete []hv;
        
        /* we are done */;
} 


Epetra_Map* h5_write_eigenfield_retrieve_DoF_robust(const Epetra_MultiVector& Q, const Epetra_Comm& comm)
{
        
        int NumMyElements_target = Q.GlobalLength();
                
//      Epetra_Map* targetmap=new Epetra_Map(-1, NumMyElements_target, 0, comm);;
        Epetra_Map* targetmap=new Epetra_Map(NumMyElements_target, NumMyElements_target, 0, comm);;

        /* we are done */
        return(targetmap);
}


Epetra_Map* h5_write_eigenfeield_retrieve_DoF_efficient(const std::string& file_name,  const Epetra_Comm& comm,
                                                                                                        mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                                                                                                        const Epetra_MultiVector& Q, double* lambda)
{
        int i,j,k,m; /* standard loop indices */
        int res; /* return code */
        int ntet; /* number of tetrahedra local to this process */
        int gntet; /* global number of tetrahedra */
        int nmode; /* store number of calulated eigenmodes */
        
        mesh::ParallelTetMesh* ptetmesh=NULL;
        ptetmesh=dynamic_cast<mesh::ParallelTetMesh*>(tetmesh);

        NedelecElement* nedelecelem=nedelecmesh.get_element(0); 
        int order=nedelecelem->get_order(); 
        int ndofelem=nedelecelem->get_nof_dofs(); 
        id_t nmyedgdofu=0; 
        id_t nmyfacdofu=0; 
        id_t nmydof=0; 
        id_t tmpid=0; 
        id_t* ldofid=NULL; /* pointer to array of local DoF id's */
        id_t edgeoffset=0; 
        id_t geid=0; 
        mesh::TetMesh::edge_iterator eit;
        mesh::TetMesh::face_iterator fit;
        mesh::TetMesh::tet_iterator tit;

        if(order == FIRST_ORDER_FE) 
        {
                nmyedgdofu=tetmesh->get_nof_edges(); 
                nmyfacdofu=0; 
                nmydof=nmyedgdofu+nmyfacdofu;
                tmpid=0;
                ldofid=new id_t[nmydof];

                for(i=0;i<nmydof;i++){ ldofid[i]=0; }

                for(eit=tetmesh->edge_begin(); eit != tetmesh->edge_end();eit++)
                {
                        ldofid[tmpid]=eit->get_id(); 
                        tmpid++;
                }
        }
        else if(order == SECOND_ORDER_FE)
        {
                nmyedgdofu=2*tetmesh->get_nof_edges(); 
                nmyfacdofu=3*tetmesh->get_nof_faces(); 
                nmydof=nmyedgdofu+nmyfacdofu;
                tmpid=0;
                ldofid=new id_t[nmydof];
                
                for(i=0;i<nmydof;i++){ ldofid[i]=0; }

                edgeoffset=tetmesh->get_nof_edges(); 
                for(eit=tetmesh->edge_begin(); eit != tetmesh->edge_end();eit++)
                {
                        ldofid[tmpid]=eit->get_id();
                        ldofid[tmpid+edgeoffset]=ldofid[tmpid]+edgeoffset;
                        tmpid++;
                }
                
        id_t face_id=0;
        int base_face_dof=0;
        
                for(fit=tetmesh->face_begin(); fit != tetmesh->face_end();fit++)
                {
                        face_id=fit->get_id();
                        base_face_dof = 3*face_id + 2*tetmesh->get_nof_edges();
                        
                        for(i=0;i<NCORNERTRIANGLE;i++)
                        {
                                ldofid[tmpid]=base_face_dof + i;
                        }
                }
        }

        int* gdofid=NULL;
        int nmapmydof=0;
        
        for(i=0;i<nmydof;i++)
        {
                if(nedelecmesh.map_dof(ldofid[i]) != mesh::ID_NONE)
                {
                        nmapmydof++;
                }
        }
        
        gdofid=new int[nmapmydof];
        j=0;
        
        for(i=0;i<nmydof;i++)
        {
                if(nedelecmesh.map_dof(ldofid[i]) != mesh::ID_NONE)
                {
                        gdofid[j]=nedelecmesh.gdof(ldofid[i]);
                        j++;
                }
        }

        Epetra_Map* targetmap=new Epetra_Map(-1,nmapmydof,gdofid,0,comm);
        
        /* we are done */
        return(targetmap);
}

void h5_write_eigenfield_write_hdf5e(const std::string& file_name,  const Epetra_Comm& comm,
                                                                        int ntet, int gntet, int mthmode, double* ev)
{
  rInfo("storing electric mode[%d] into HDF5 output file",mthmode);
  
        int i,j,k; /* standard loop indices */
        int res; /* standard return code */
        int ncol=NREALCOORD3D; 
        hid_t       file_id, group_id, dset_id;       /* file and dataset identifiers */
        hid_t       filespace, memspace;              /* file and memory dataspace identifiers */
        hid_t       plist_id;                         /* property list identifier */
        herr_t      status;
        std::string mtheigenfieldbasename(H5_NAME_FEMAX_EIGENFIELD); /* base dataset name

        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
        
        int mpirank=0; 
        int mpisize=0; 
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << " reached line" << __LINE__  << " reached " << endl;
#endif 
        plist_id = H5Pcreate(H5P_FILE_ACCESS);
        H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

  file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
  H5Pclose(plist_id);

  rInfo("opening group in HDF5 file = %s",H5_GROUP_FEMAX_EIGENMODES.c_str());
  group_id = H5Gopen(file_id, H5_GROUP_FEMAX_EIGENMODES.c_str());
  rInfo("opened group in HDF5 file = %s",H5_GROUP_FEMAX_EIGENMODES.c_str());

//      if(group_id < 0) /* there is no such group yet, create one */
//      {
//              group_id = H5Gcreate(file_id, H5_GROUP_FEMAX_EIGENMODES.c_str(), 0);
//      }

        hsize_t dimsf[]={gntet,ncol};
        filespace = H5Screate_simple(DIM2D, dimsf, NULL);

        std::ostringstream out;
        out << mthmode;
        mtheigenfieldbasename.append(out.str());
                
        dset_id = H5Dcreate(group_id, mtheigenfieldbasename.c_str(), H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT);
        H5Sclose(filespace);

        hsize_t count[2];
        hsize_t offset[2];
        
        count[0]=ntet;
        count[1]=dimsf[1]; 
        id_t* ntetlist=new id_t[mpisize]; 
        res=MPI_Allgather(&ntet,1,MPI_UNSIGNED, ntetlist, 1, MPI_UNSIGNED, mpi_comm);

        offset[0]=0;
        for(i=0;i<mpirank;i++)
        {
                offset[0]+=ntetlist[i];
        }

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " offset[0] " << offset[0] << endl;
#endif 
        offset[1]=0; /* no offset w.r.t columns */
        memspace=H5Screate_simple(DIM2D,count,NULL);

        filespace=H5Dget_space(dset_id);
        H5Sselect_hyperslab(filespace, H5S_SELECT_SET,offset,NULL,count,NULL);
        
        plist_id=H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
        status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace,plist_id, ev);

        H5Dclose(dset_id);
        H5Sclose(filespace);
        H5Sclose(memspace);
        H5Gclose(group_id);
        H5Fclose(file_id);

        /* deletey dynamically allocated variables if any */
        delete []ntetlist;
        
        /* we are done */
  rInfo("stored electric mode[%d] into HDF5 output file",mthmode);
}

void h5_write_eigenfield_write_hdf5h(const std::string& file_name,  const Epetra_Comm& comm,
                                                                        int ntet, int gntet, int mthmode, double* hv)
{
  rInfo("storing magnetic field of mode[%d] into HDF5 output file",mthmode);
  
        int i,j,k; /* standard loop indices */
        int res; /* standard return code */
        int ncol=NREALCOORD3D; 
        hid_t       file_id, group_id, dset_id;       /* file and dataset identifiers */
        hid_t       filespace, memspace;              /* file and memory dataspace identifiers */
        hid_t       plist_id;                         /* property list identifier */
        herr_t      status;
        std::string mtheigencurlbasename(H5_NAME_FEMAX_EIGENCURL); /* base dataset name

        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
        
        int mpirank=0; 
        int mpisize=0; 
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << " reached line" << __LINE__  << " reached " << endl;
#endif 
        plist_id = H5Pcreate(H5P_FILE_ACCESS);
        H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

  file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
  H5Pclose(plist_id);

  rInfo("opening group in HDF5 file = %s",H5_GROUP_FEMAX_EIGENMODES.c_str());
  group_id = H5Gopen(file_id, H5_GROUP_FEMAX_EIGENMODES.c_str());
  rInfo("opened group in HDF5 file = %s",H5_GROUP_FEMAX_EIGENMODES.c_str());
  
        hsize_t dimsf[]={gntet,ncol};
        filespace = H5Screate_simple(DIM2D, dimsf, NULL);

        std::ostringstream out;
        out << mthmode;
        mtheigencurlbasename.append(out.str());
                
        dset_id = H5Dcreate(group_id, mtheigencurlbasename.c_str(), H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT);
        H5Sclose(filespace);

        hsize_t count[2];
        hsize_t offset[2];
        
        count[0]=ntet;
        count[1]=dimsf[1]; 
        id_t* ntetlist=new id_t[mpisize]; 
        res=MPI_Allgather(&ntet,1,MPI_UNSIGNED, ntetlist, 1, MPI_UNSIGNED, mpi_comm);

        offset[0]=0;
        for(i=0;i<mpirank;i++)
        {
                offset[0]+=ntetlist[i];
        }

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " offset[0] " << offset[0] << endl;
#endif 
        offset[1]=0; /* no offset w.r.t columns */
        memspace=H5Screate_simple(DIM2D,count,NULL);

        filespace=H5Dget_space(dset_id);
        H5Sselect_hyperslab(filespace, H5S_SELECT_SET,offset,NULL,count,NULL);
        
        plist_id=H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
        status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace,plist_id, hv);

        H5Dclose(dset_id);
        H5Sclose(filespace);
        H5Sclose(memspace);
        H5Gclose(group_id);
        H5Fclose(file_id);

        /* deletey dynamically allocated variables if any */
        delete []ntetlist;
        
        /* we are done */
  rInfo("stored magnetic mode[%d] into HDF5 output file",mthmode);
}


void h5_write_eigenfield_writehdf5sloc(const std::string& file_name,
                                                                const Epetra_Comm& comm,
                                                                int ntet,
                                                                int gntet,
                                                                double* slocv)
{
        int i,j,k; /* standard loop indices */
        int res; /* standard return code */
        int ncol=NREALCOORD3D; 
        hid_t       file_id, group_id, dset_id;       /* file and dataset identifiers */
        hid_t       filespace, memspace;              /* file and memory dataspace identifiers */
        hid_t       plist_id;                         /* property list identifier */
        herr_t      status;
        std::string slocbasename(H5_NAME_FEMAX_SLOCBASENAME); /* base dataset name

        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
        
        int mpirank=0; 
        int mpisize=0; 
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << " reached line" << __LINE__  << " reached " << endl;
#endif 
        plist_id = H5Pcreate(H5P_FILE_ACCESS);
        H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
    H5Pclose(plist_id);

    group_id = H5Gopen(file_id, H5_GROUP_FEMAX_EIGENMODES.c_str());
  
        if(group_id < 0) /* there is no such group yet, create one */
        {
                group_id = H5Gcreate(file_id, H5_GROUP_FEMAX_EIGENMODES.c_str(), 0);
        }

        hsize_t dimsf[]={gntet,ncol};
        filespace = H5Screate_simple(DIM2D, dimsf, NULL);

        dset_id = H5Dcreate(group_id, slocbasename.c_str(), H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT);
        H5Sclose(filespace);

        hsize_t count[2];
        hsize_t offset[2];
        
        count[0]=ntet;
        count[1]=dimsf[1]; 
        id_t* ntetlist=new id_t[mpisize]; 
        res=MPI_Allgather(&ntet,1,MPI_UNSIGNED, ntetlist, 1, MPI_UNSIGNED, mpi_comm);

        offset[0]=0;
        for(i=0;i<mpirank;i++)
        {
                offset[0]+=ntetlist[i];
        }

#ifdef __DEBUG__VERBOSE__
        cout << "[DEBUGOUT]  : " << __FILE__  << " : " << __LINE__  << " offset[0] " << offset[0] << endl;
#endif 
        offset[1]=0; /* no offset w.r.t columns */
        memspace=H5Screate_simple(DIM2D,count,NULL);

        filespace=H5Dget_space(dset_id);
        H5Sselect_hyperslab(filespace, H5S_SELECT_SET,offset,NULL,count,NULL);
        
        plist_id=H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
        status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace,plist_id, slocv);

        H5Dclose(dset_id);
        H5Sclose(filespace);
        H5Sclose(memspace);
        H5Gclose(group_id);
        H5Fclose(file_id);

        /* deletey dynamically allocated variables if any */
        delete []ntetlist;
        
        /* we are done */       
        
        
}


void h5_write_eigenvalue(const std::string& file_name,  const Epetra_Comm& comm,
                                                                                mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                                                                                const Epetra_MultiVector& Q, double* lambda)
{
        int i,j,k; /* standard loop indices */
        int res; /* standard return code */
        int ncol=3; /* number of columns to be written; at present we write the eigenvalue lambda, its associated eigenfrequency f [Hz] and the quality factor */
        hid_t file_id, group_id, dset_id; /* file and dataset identifiers */
        hid_t dataspace, filespace, memspace; /* file and memory dataspace identifiers */
        hid_t plist_id; /* property list identifier */
        herr_t status; /* error signal */
        
        int mpirank=0; 
        int mpisize=0; 
        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 
        
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);
        
        /* go and compute the quality factors of the respective modes */
        int nlambda=Q.NumVectors();
        double* quality=new double[nlambda];
        for(i=0;i<nlambda;i++){quality[i]=0.0; }
        h5_compute_eigenquality(file_name, comm, tetmesh, nedelecmesh, Q, lambda, quality);
        
        if(mpirank==MPIROOTPROCESS) /* if I am the root process then I write the eigenvalue array into the HDF5 file */
        {
    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
  
    group_id = H5Gcreate(file_id, H5_GROUP_FEMAX_EIGENMODES.c_str(), 0);
                
                /* handle HDF5 related stuff for storing the array */
                hsize_t count[]={Q.NumVectors(),ncol};
                hsize_t dimsf[]={Q.NumVectors(),ncol};
        
                /* store eigenvalues and associated eigenfrequencies in 2d array */
                int nlambda=dimsf[0];
                double* lambdafres=new double[nlambda*ncol];
                for(i=0;i<nlambda;i++)
                {
                        /* store eigenvalues into 2d array */
                        lambdafres[i*ncol]=lambda[i];
                        
                        /* calculate associated resonance frequency and store it into the 2d array */
                        lambdafres[i*ncol+1]=sqrt(lambda[i])*(SPEED_OF_LIGHT_VACUUM/(2*PI));
                        
                        /* store the associated quality factir into the 2d array */
                        lambdafres[i*ncol+2]=quality[i];
                }
                
                dataspace = H5Screate_simple(DIM2D, dimsf, NULL);               
                dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_EIGENVALUE.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
                H5Sclose(dataspace);
                
                status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, lambdafres);
                
                /* close down HDF5 operations */
                H5Dclose(dset_id);
                H5Gclose(group_id);
                H5Fclose(file_id);
                
                /* delete dynamically allocate memory if any */
                delete []lambdafres;
                        
        }
        
        /* we are done */
}

void h5_compute_eigenquality(const std::string& file_name,  const Epetra_Comm& comm,
                                                                                mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                                                                                const Epetra_MultiVector& Q, double* lambda, double* &quality)
{
        int i,j,k,m; /* standard loop indices */
        int res; /* return code */
        int ntet; /* number of tetrahedra local to this process */
        int gntet; /* global number of tetrahedra */
        int nmode=Q.NumVectors(); /* store number of calulated eigenmodes */
  
  rInfo("computing quality factor of eigenmodal solutions for %d modes",nmode);
  
        mesh::ParallelTetMesh* ptetmesh=dynamic_cast<mesh::ParallelTetMesh*>(tetmesh);

        NedelecElement* nedelecelem=nedelecmesh.get_element(0); 
        MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
        MPI_Info mpi_info  = MPI_INFO_NULL; 

        int mpirank=0; 
        int mpisize=0; 
        res=MPI_Comm_rank(mpi_comm, &mpirank);
        res=MPI_Comm_size(mpi_comm, &mpisize);

        int order=nedelecelem->get_order(); 
        int ndofelem=nedelecelem->get_nof_dofs(); 
        /* compute an Epetra_Map */
        Epetra_Map* targetmap=h5_write_eigenfield_retrieve_DoF_robust(Q,comm); /* the robust way */
                                                                                                                
        Epetra_MultiVector Q0(*targetmap, Q.NumVectors());
        
        Epetra_Import importer(*targetmap,Q.Map());
        
        Q0.Import(Q,importer,Insert);

        /* nested loop over all modes and all faces for computing the surface integrals */
        mesh::id_t tetid=mesh::ID_NONE;
        mesh::TetMesh::surface_iterator fit;
        mesh::Face* face=NULL;
        const int* tldofid=NULL;
        int tlid=0;
        int mapped_id=0;
        int nface=0;
        
        double storedenergy=0.0;
        
        double* integral=new double[nmode]; /* provide memory for the respective surface integral contributions of the processes */
        double rintegral=0.0; /* used for storing result of reduction operation */
        double tintegral=0.0;
        
        double* omega=new double[nmode]; /* provide memory for the angular resonance frequencies */
        double* fres=new double[nmode]; /* provide memory for the resonance frequencies in Hz */
        double* skindepth=new double[nmode]; /* provide memory for the skin depths at the respective resonance frequencies */
        double* surfaceresistance=new double[nmode]; /* provide memory for surfaceresistance */
        double* avgpowerloss=new double[nmode]; /* provide memory for the average power losses */
        
        id_t t0,t1; /* store id's of tetrahedra adjacent to a face */
        t0=t1=mesh::ID_NONE;
        
        /* loop over modes */
        for(m=0;m<nmode;m++)
        {
//    rInfo("computing mode[%d] quality factor",m);
    
                Epetra_Vector mthmode(Copy, Q0,m);
                tintegral=0.0; /* initialize integral value */
                rintegral=0.0; /* initiaize reduced integral value */
                nface=0; /* initialize number of faces looped over */

                for(fit=tetmesh->surface_begin(); fit != tetmesh->surface_end(); ++fit)
                {
                        face = *fit;

                        /* find tet and elem belonging to boundary face */
                        tetid = face->get_tet_id();
                        nedelecelem = nedelecmesh.get_element(tetid);

                        /* The surface integral evaluates the n x H over all boundary faces;
                         * there are different types of boundary faces, i.e. PEC, PMC, symmetry
                         * boundary faces and interpocessor boundary faces; faces of the type
                         * PEC, PMC and symmetry only have one adjacent tetrahedron while
                         * faces of the interprocessor type have actually two boundary
                         * faces; therefore, we must make sure that the faces for which the
                         * integral is evaluated, has only one single adjacent tetrahedron.
                         * 
                         */
                        face->get_tet_ids(&t0,&t1);
                        if(( (face->on_boundary() == true ) && ( face->on_symmetry_planes() == 0) ) )
                        {
                                /* obtain memory for storing DoF */
                                colarray::Vector<double> dofs(nedelecelem->get_nof_dofs());
                                                                
                                tldofid=nedelecelem->get_dof_ids(tetmesh,tetmesh->get_tet(tetid));
                
                                for (i=0; i < nedelecelem->get_nof_dofs(); i ++) 
                                {
                                        mapped_id=nedelecmesh.map_dof(tldofid[i]);
                        
                                        if (mapped_id == mesh::ID_NONE)
                                        {
                                                dofs(i) = 0.0;
                                        }
                                        else
                                        {                               
                                                tlid=targetmap->LID(mapped_id);
                                                dofs(i) = mthmode[tlid];
                                        }
        }
                        
                                /* calculate the contribution of a face to the integral: int_(face)(curl N_i * curl N_j, ds) */
                                tintegral += nedelecelem->surface_integral_curl(tetmesh->get_tet(tetid), face, dofs);
                                
                                /* increase number of processed faces */
                                nface++;
                        }
                }
        
    /* 1st, we need to reduce the processes' contribution to the surface integral */
//    rInfo("collecting contributions for integral from all processors");
    res=MPI_Allreduce(&tintegral,&rintegral,1, MPI_DOUBLE, MPI_SUM, mpi_comm);
    integral[m]=rintegral; /* save result of modewise integration */
    rInfo("computed mode[%d] quality factor",m);
  
        }
        
        /* compute the quality factor */
        double cintegral=0.0; /*DEBUGING USE ONLY */
        double cquality=0.0;  /* DEBUGGING USE ONLY */
        double cfres=0.0;/*DEBUGING USE ONLY */
        double cskindepth=0.0; /*DEBUGING USE ONLY */
        double csurfaceresistance=0.0; /*DEBUGING USE ONLY */
        double cavgpowerloss=0.0; /*DEBUGING USE ONLY */
        
        /* loop over the modal surface integrals and compute resonance freuencies, skindepths
         * and quality factors */
        for(m=0;m<nmode;m++)
        {
                omega[m]=sqrt(lambda[m])*SPEED_OF_LIGHT_VACUUM;
                fres[m]=sqrt(lambda[m])*SPEED_OF_LIGHT_VACUUM/(2.0*PI);
                skindepth[m]=1.0/(sqrt(PI*fres[m]*MU_ZERO*SIGMA_CU_PSI)); /* Ramo et al., p. 151, eq. 3.16(11) */
                surfaceresistance[m]=(1.0/(SIGMA_CU_PSI*skindepth[m])); /* Ramo et. al., p. 154, eq. 3.17(5) */
                integral[m] = integral[m] / ( MU_ZERO * MU_ZERO *omega[m] * omega[m]);
                avgpowerloss[m]=(surfaceresistance[m] / 2.0) * integral[m]; /* Ramo et al. p. 446, eq. 8.13(19) */
                
//                      const int nofSym = 0;
//                      avgpowerloss[m]*= (1 << nofSym);
                
                storedenergy=EPSILON_ZERO / 2.0;
                quality[m]=omega[m] * storedenergy / avgpowerloss[m]; /* Ramo et al., p. 258. eq. 5.14(4) */
                
                /*DEBUGING USE ONLY */  
                cintegral=integral[m];
                cquality=quality[m];
                cfres=fres[m];
                cskindepth=skindepth[m];
                csurfaceresistance=surfaceresistance[m];
                cavgpowerloss=avgpowerloss[m];
                /*DEBUGING USE ONLY */
                
        }
                
                
//    const int nofSym = 0;           // nofSym is not necessary since it cancels out
//    integral *= (1 << nofSym);
//    double omega = ::sqrt(lambda) * cLight;
//    // double f = omega / 2.0 / pi;
//    double intJpower2 = integral / mu0 / mu0 / omega / omega;
//    double delta = sqrt(2.0 / (mu0 * omega * sigmaCu));
//    double reZs = 1.0 / sigmaCu / delta;
//    double W = reZs / 2 * intJpower2;          // Wall power loss
//    double U = (1 << nofSym) * eps0 / 2;       // stored total energy
//    double Q = omega * U / W;                  // unloaded quality factor
                
        
        /* free dynamically allocated memory if any */
        delete []integral;
        delete []omega;
        delete []fres;
        delete []skindepth;
        delete []surfaceresistance;
        delete []avgpowerloss;
                                
                //    mesh::TetMesh::surface_iterator it;
//    for (it = _mesh->surface_begin(); it != _mesh->surface_end(); ++ it) {
//        Face* face = *it;
//
//        // find tet and elem belonging to boundary face
//        mesh::id_t tet = face->get_tet_id();
//        NedelecElement* elem = get_element(tet);
//        // dofs for elem
//        colarray::Vector<double> dofs(elem->get_nof_dofs());
//        for (int i = 0; i < elem->get_nof_dofs(); ++ i) {
//            const int* dof_ids = elem->get_dof_ids(_mesh, _mesh->get_tet(tet));
//            int mapped_id = map_dof(dof_ids[i]);
//            if (mapped_id == -1)
//                dofs(i) = 0.0;
//            else
//                dofs(i) = q[mapped_id];
//        }
//        // calculate int_(face)(curl N_i * curl N_j, ds)
//        integral += elem->surface_integral_curl(_mesh->get_tet(tet), face, dofs);
//    }
                
                
                
                
        

        
//    double integral = 0.0;
//
//    mesh::TetMesh::surface_iterator it;
//    for (it = _mesh->surface_begin(); it != _mesh->surface_end(); ++ it) {
//        Face* face = *it;
//
//        // find tet and elem belonging to boundary face
//        mesh::id_t tet = face->get_tet_id();
//        NedelecElement* elem = get_element(tet);
//        // dofs for elem
//        colarray::Vector<double> dofs(elem->get_nof_dofs());
//        for (int i = 0; i < elem->get_nof_dofs(); ++ i) {
//            const int* dof_ids = elem->get_dof_ids(_mesh, _mesh->get_tet(tet));
//            int mapped_id = map_dof(dof_ids[i]);
//            if (mapped_id == -1)
//                dofs(i) = 0.0;
//            else
//                dofs(i) = q[mapped_id];
//        }
//        // calculate int_(face)(curl N_i * curl N_j, ds)
//        integral += elem->surface_integral_curl(_mesh->get_tet(tet), face, dofs);
//    }
//  
//    const int nofSym = 0;           // nofSym is not necessary since it cancels out
//    integral *= (1 << nofSym);
//    double omega = ::sqrt(lambda) * cLight;
//    // double f = omega / 2.0 / pi;
//    double intJpower2 = integral / mu0 / mu0 / omega / omega;
//    double delta = sqrt(2.0 / (mu0 * omega * sigmaCu));
//    double reZs = 1.0 / sigmaCu / delta;
//    double W = reZs / 2 * intJpower2;          // Wall power loss
//    double U = (1 << nofSym) * eps0 / 2;       // stored total energy
//    double Q = omega * U / W;                  // unloaded quality factor
//#if 0
//    cout << "lambda=" << lambda << endl;
//    cout << "integral=" << integral << endl;
//    cout << "omega=" << omega << endl;
//    cout << "intJpower2=" << intJpower2 << endl;
//    cout << "delta=" << delta << endl;
//    cout << "reZs=" << reZs << endl;
//    cout << "W=" << W << endl;
//    cout << "U=" << U << endl;
//#endif
//    return Q;
        
        
        
        
  rInfo("computed quality factor of eigenmodal solutions");
        
        
}


void h5_cartesian_sampling(const std::string& file_name,  const Epetra_Comm& comm,
                    mesh::TetMesh* tetmesh, NedelecMesh& nedelecmesh,
                    const Epetra_MultiVector& Q, double* lambda,
                    double xmin, double xmax, const int nx,
                    double ymin, double ymax, const int ny,
                    double zmin, double zmax, const int nz)
{
  rInfo("evaluating eigenmodal solution on cartesian grid for %d modes",Q.NumVectors());
  
  mesh::ParallelTetMesh* ptetmesh=dynamic_cast<mesh::ParallelTetMesh*>(tetmesh);
  
  NedelecElement* nedelecelem=nedelecmesh.get_element(0); 
  int order=nedelecelem->get_order(); 
  int ndofelem=nedelecelem->get_nof_dofs(); 
  /* compute an Epetra_Map */
  Epetra_Map* targetmap=h5_write_eigenfield_retrieve_DoF_robust(Q,comm); /* the robust way */
                            
  Epetra_MultiVector Q0(*targetmap, Q.NumVectors());
  
  Epetra_Import importer(*targetmap,Q.Map());

  Q0.Import(Q,importer,Insert);


  std::vector<std::vector<std::vector<unsigned int> > > positionExists;
  positionExists.resize(nx);
  
  for(std::vector<std::vector<std::vector<unsigned int> > >::iterator it=positionExists.begin(); it != positionExists.end();it++)
  {
    (*it).resize(ny);
    for(std::vector<std::vector<unsigned int> >::iterator it2=(*it).begin(); it2 != (*it).end(); it2++)
    {
      (*it2).resize(nz);
      for(std::vector<unsigned int>::iterator it3=(*it2).begin(); it3 != (*it2).end(); it3++)
      {
        (*it3)=0;
      }
    }
  }


  unsigned int nmode=Q.NumVectors();
  
  std::vector<std::vector<std::vector<std::vector<mesh::Vector3> > > > efield;
  std::vector<std::vector<std::vector<std::vector<mesh::Vector3> > > > hfield;

    efield.resize(nmode);
    hfield.resize(nmode);
    
    for(unsigned int m=0;m<nmode;m++)
    {
      efield[m].resize(nx);
      hfield[m].resize(nx);
      
      for(unsigned int i=0;i<nx;i++)
      {
        efield[m][i].resize(ny);
        hfield[m][i].resize(ny);
  
        for(unsigned int j=0;j<ny;j++)
        {
          efield[m][i][j].resize(nz);
          hfield[m][i][j].resize(nz);
          
          for(unsigned int k=0;k<nz;k++)
          {
            efield[m][i][j][k]=mesh::Vector3(0,0,0); 
            hfield[m][i][j][k]=mesh::Vector3(0,0,0); 
          }        
        }
      }
    }




  rInfo("constructing octree for the terahedral mesh");
  tetmesh->construct_octree();

  rInfo("sampling eigenmodal solution on cartesian grid for %d modes", nmode);
    
  double dx=(xmax-xmin)/((double)(nx-1));
  double dy=(ymax-ymin)/((double)(ny-1));
  double dz=(zmax-zmin)/((double)(nz-1));
  
  rInfo("cartesian sampling interval on x-axis dx = %12.9e", dx);
  rInfo("cartesian sampling interval on y-axis dx = %12.9e", dy);
  rInfo("cartesian sampling interval on z-axis dx = %12.9e", dz);

//  Epetra_Vector* mthmode=NULL;
//  for(unsigned int m=0;m<nmode;m++)
//  {
//    mthmode[m]=new Epetra_Vector;
//    mthmode[m](Copy, Q0,m);
//    
//  }

  std::vector<Epetra_Vector> mthmode;
    for(unsigned int m=0;m<nmode;m++)
    {
      mthmode.push_back(Epetra_Vector(Copy, Q0,0)); 
    }



  for(unsigned int i=0;i<nx;i++)
  {
    for(unsigned int j=0;j<ny;j++)
    {
      for(unsigned int k=0;k<nz;k++)
      {
        double x=xmin+(i*dx);
        double y=ymin+(j*dy);
        double z=zmin+(k*dz);
        
        mesh::Vector3 pos(x,y,z);
        
        std::vector<mesh::Tet*> tet;
        
        tetmesh->find_tets_by_point(pos,tet);
        
        unsigned int ntetfound=tet.size();
 
        if(ntetfound != 0)
        {
          positionExists[i][j][k]=1;

          for(unsigned int m=0; m < nmode; m++)
          {
            
            const int* tldofid=NULL;
            colarray::Vector<double> dofs(nedelecelem->get_nof_dofs());
            int mapped_id=0;
            mesh::Tet* ttet=tet[0];
            int tlid=0;

            nedelecelem=nedelecmesh.get_element(ttet->get_id());
            tldofid=nedelecelem->get_dof_ids(tetmesh,ttet);

            for (unsigned int d=0; d < nedelecelem->get_nof_dofs(); d++) 
            {
              mapped_id=nedelecmesh.map_dof(tldofid[d]);
            
              if (mapped_id == mesh::ID_NONE)
              {
                dofs(d) = 0.0;
              }
              else
              {       
                tlid=targetmap->LID(mapped_id);
                dofs(d) = mthmode[m][tlid];
              }
            }
            
            efield[m][i][j][k].x=nedelecelem->eval(ttet, pos, dofs).x;
            efield[m][i][j][k].y=nedelecelem->eval(ttet, pos, dofs).y;          
            efield[m][i][j][k].z=nedelecelem->eval(ttet, pos, dofs).z;
                        
            hfield[m][i][j][k].x=nedelecelem->eval_curl(ttet, pos, dofs).x;
            hfield[m][i][j][k].y=nedelecelem->eval_curl(ttet, pos, dofs).y; 
            hfield[m][i][j][k].z=nedelecelem->eval_curl(ttet, pos, dofs).z;
          }
        }
        else
        {
          positionExists[i][j][k]=0;
        }
      }   
    }   
  }   
  rInfo("now going to write sampled field into HDF5 file", nmode);
  h5_cartesian_sampling_write_hdf5(comm, file_name, nmode,
                    positionExists,
                    efield,
                    hfield,
                    xmin, xmax, nx,
                    ymin, ymax, ny,
                    zmin, zmax, nz);
}

void h5_cartesian_sampling_write_hdf5(const Epetra_Comm& comm, const std::string& file_name, unsigned int nmode,
                    std::vector<std::vector<std::vector<unsigned int> > >& exists,
                    std::vector<std::vector<std::vector<std::vector<mesh::Vector3> > > >& efield,
                    std::vector<std::vector<std::vector<std::vector<mesh::Vector3> > > >& hfield,
                    double xmin, double xmax, const int nx,
                    double ymin, double ymax, const int ny,
                    double zmin, double zmax, const int nz)
{
  rInfo("storing eigenmodal solution sampled on cartesian grid into hdf5 file");  
  
  MPI_Comm mpi_comm=dynamic_cast<const Epetra_MpiComm&>(comm).GetMpiComm();
  MPI_Info mpi_info  = MPI_INFO_NULL; 
  
  int mpirank=0; /* store the rank of this prcocess within this MPI communicator */
  int mpisize=0; /* store the size of an MPI communicator, i.e. the number of processes engaged */
  int res;       /* standard return code */ 
  
  res=MPI_Comm_rank(mpi_comm, &mpirank);
  res=MPI_Comm_size(mpi_comm, &mpisize);

  hid_t file_id, group_id, dset_id, plist_id;               /* file and dataset identifiers */
  hid_t dataspace;                                          /* file and memory dataspace identifiers */
  herr_t status;                                                                                  /* error signal */


  if(mpirank==MPIROOTPROCESS) /* if I am the root process then I initialize cartesian sampling storage */
  {
    rInfo("initializing cartesian sampling storage of eigenmodal solution");  
      
    file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
  
    group_id = H5Gcreate(file_id, H5_GROUP_NAME_EIGENMODES_CARTESIAN.c_str(), 0);
    rInfo("creating HDF5 group = %s",H5_GROUP_NAME_EIGENMODES_CARTESIAN.c_str());  
    
    hsize_t dimsf[]={3,2};    
    dataspace = H5Screate_simple(DIM2D, dimsf, NULL);
    
    dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_EIGENOMDES_CARTESIAN_LIMITS.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
    H5Sclose(dataspace);
    
    double limit[]={xmin,xmax,ymin,ymax,zmin,zmax};
    status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, limit); 
    
    H5Dclose(dset_id);

    dimsf[0]=3;
    dimsf[1]=1;
    
    dataspace = H5Screate_simple(DIM2D, dimsf, NULL);
    
    dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_EIGENOMDES_CARTESIAN_NUMBER_SAMPLES.c_str(),
                        H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
                        
    H5Sclose(dataspace);
    
    double nsample[]={nx,ny,nz};
    status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nsample); 
    
    H5Dclose(dset_id);

    dimsf[0]=1;
    dimsf[1]=1;
    
    dataspace = H5Screate_simple(DIM2D, dimsf, NULL);
    
    dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_EIGENOMDES_CARTESIAN_NMODE.c_str(),
                        H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
                        
    H5Sclose(dataspace);
    
    double nmode_[]={nmode};
    status=H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, nmode_); 
    
    H5Dclose(dset_id);

    rInfo("now initializing the array which indicates if a position exists");  

    hsize_t dims3Darray[]={nx,ny,nz};
    
    dataspace = H5Screate_simple(DIM3D, dims3Darray, NULL);

    dset_id = H5Dcreate(group_id, H5_NAME_FEMAX_EIGENOMDES_CARTESIAN_EXISTS.c_str(),
                        H5T_NATIVE_UINT, dataspace, H5P_DEFAULT);

    H5Sclose(dataspace);

    unsigned int nbool=nx*ny*nz;
    unsigned int* exists_=new unsigned int[nbool];
    for(unsigned int i=0;i<nbool;i++){exists_[i]=0;}

    status=H5Dwrite(dset_id, H5T_NATIVE_UINT, H5S_ALL, H5S_ALL, H5P_DEFAULT, exists_);

    H5Dclose(dset_id);
    H5Gclose(group_id);
    H5Fclose(file_id);
  }
  
  hid_t dset_id_electric, dset_id_magnetic;                 /* file and dataset identifiers */
//  hid_t filespace_electric, filespace_magnetic;             /* file and memory dataspace identifiers */
  hid_t memspace_electric, memspace_magnetic;               /* file and memory dataspace identifiers */
  hid_t dataspace_electric, dataspace_magnetic;             /* file and memory dataspace identifiers */
  
  
  rInfo("transferring cartesian eigenmodal solution into HDF5 file");  

  unsigned int nlcartesian=0;   
  for(unsigned int i=0;i<nx;i++)
  {
    for(unsigned int j=0;j<ny;j++)
    {
      for(unsigned int k=0;k<nz;k++)
      {
        if(exists[i][j][k] != 0)
        {
          nlcartesian++;
        }     
      }
    }
  }

  rInfo("#[actually used local cartesian samples]=%d",nlcartesian);  

  plist_id = H5Pcreate(H5P_FILE_ACCESS);
  H5Pset_fapl_mpio(plist_id, mpi_comm, mpi_info);

  file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, plist_id);
  H5Pclose(plist_id);

  rInfo("opening HDF5 group = %s",H5_GROUP_NAME_EIGENMODES_CARTESIAN.c_str());  
  group_id = H5Gopen(file_id, H5_GROUP_NAME_EIGENMODES_CARTESIAN.c_str());


  for(unsigned int m=0;m<nmode;m++)
  {
    std::ostringstream out;
    out << m;
    std::string electricMthCartEigenfield(H5_NAME_FEMAX_EIGENOMDES_CARTESIAN_EFIELD_MODE);
    std::string magneticMthCartEigenfield(H5_NAME_FEMAX_EIGENOMDES_CARTESIAN_HFIELD_MODE);  
    
    electricMthCartEigenfield.append(out.str());
    magneticMthCartEigenfield.append(out.str());

//    hsize_t dims4[]={nx,ny,nz,DIM3D};
    unsigned int ncsample=nx*ny*nz;
    hsize_t dim2[]={ncsample,DIM3D};
    
//    rInfo("dim2[0] = %d ; dim2[1] = %d",dim2[0],dim2[1]);  
    
    dataspace_electric = H5Screate_simple(DIM2D, dim2, NULL);
    dataspace_magnetic = H5Screate_simple(DIM2D, dim2, NULL);

    dset_id_electric = H5Dcreate(group_id, electricMthCartEigenfield.c_str(), H5T_NATIVE_DOUBLE, dataspace_electric, H5P_DEFAULT);
    dset_id_magnetic = H5Dcreate(group_id, magneticMthCartEigenfield.c_str(), H5T_NATIVE_DOUBLE, dataspace_magnetic, H5P_DEFAULT);

    H5Sclose(dataspace_electric);
    H5Sclose(dataspace_magnetic);
    
    double* xfer3dEfield=new double[nlcartesian*NREALCOORD3D];
    double* xfer3dHfield=new double[nlcartesian*NREALCOORD3D];
    
    for(unsigned int i=0;i<(nlcartesian*NREALCOORD3D);i++)
    {
      xfer3dEfield[i]=0.0;
      xfer3dHfield[i]=0.0;
    }

    hsize_t count[]={nlcartesian,DIM3D};
    
//    rInfo("count[0] = %d ; count[1] = %d",count[0],count[1]);  
    
    
    hsize_t offset[DIM2D];
    
    
    
    memspace_electric=H5Screate_simple(DIM2D,count,NULL);
    memspace_magnetic=H5Screate_simple(DIM2D,count,NULL);
        
    dataspace_electric=H5Dget_space(dset_id_electric);
    dataspace_magnetic=H5Dget_space(dset_id_magnetic);
  
    hsize_t scount[2];
    scount[0]=1;
    scount[1]=NREALCOORD3D;
  
    unsigned int nlcartesianprocessed=0;  /* count number of locally processed cartesian sampling locations */
    
    for(unsigned int i=0;i<nx;i++)
    {
      for(unsigned int j=0;j<ny;j++)
      {
        for(unsigned int k=0;k<nz;k++)
        {
          
          if(exists[i][j][k] != 0)
          {
            offset[0]=(i * ny * nz) + (j * nz ) + k;
            offset[1]=0;

            if(nlcartesianprocessed == 0) /* the very first hyperslab uses a different logic for construction */
            {
              H5Sselect_hyperslab(dataspace_electric, H5S_SELECT_SET,offset,NULL,scount,NULL);
              H5Sselect_hyperslab(dataspace_magnetic, H5S_SELECT_SET,offset,NULL,scount,NULL);         
            }
            else
            {
              H5Sselect_hyperslab(dataspace_electric, H5S_SELECT_OR,offset,NULL,scount,NULL);
              H5Sselect_hyperslab(dataspace_magnetic, H5S_SELECT_OR,offset,NULL,scount,NULL);
            }
            
            for(unsigned int n=0;n<NREALCOORD3D;n++)
            {
              xfer3dEfield[nlcartesianprocessed*NREALCOORD3D+n]=efield[m][i][j][k][n];
              xfer3dHfield[nlcartesianprocessed*NREALCOORD3D+n]=hfield[m][i][j][k][n];
            }
            
            nlcartesianprocessed++;
          }
           
        }
      }
    }
    
    plist_id=H5Pcreate(H5P_DATASET_XFER);
    H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
    
    status=H5Dwrite(dset_id_electric, H5T_NATIVE_DOUBLE, memspace_electric, dataspace_electric, plist_id, xfer3dEfield);
//    status=H5Dwrite(dset_id_magnetic, H5T_NATIVE_DOUBLE, memspace_magnetic, filespace_magnetic, plist_id, xfer3dHfield);
    
    delete []xfer3dEfield;
    delete []xfer3dHfield;
    

    H5Pclose(plist_id);

    H5Sclose(memspace_electric);
    H5Sclose(memspace_magnetic);
    
    H5Dclose(dset_id_electric);
    H5Dclose(dset_id_magnetic);
    
  }


  H5Gclose(group_id);
  H5Fclose(file_id);  
}

---