SF_mesh_io.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#ifndef _SF_MESH_IO
#define _SF_MESH_IO
 
#include <cstddef>
#include <list>
#include <string>
#include <cstring>
#include <unistd.h>
 
#include <mpi.h>
 
#include "SF_container.h"
#include "SF_io_base.h"
#include "SF_network.h"
#include "SF_globals.h"
#include "SF_vector.h"
 
#define HDR_SIZE 1024
 
namespace SF {
 
inline size_t read_num_pts(std::string basename)
{
  size_t numpts;
 
  bool read_binary = fileExists(basename + ".bpts");
  std::string pts_file = read_binary ? basename + ".bpts" : basename + ".pts";
 
  FILE* pts_fd = fopen(pts_file.c_str(), "r");
  if(!pts_fd) {
    fprintf(stderr, "Error: could not open file: %s. Aborting!\n", pts_file.c_str());
    exit(1);
  }
 
  // read number of points
  char buffer[2048];
  if(read_binary) fread(buffer, sizeof(char), HDR_SIZE, pts_fd);
  else fgets( buffer, 2048, pts_fd);
  sscanf(buffer, "%lu", &numpts);
  fclose(pts_fd);
 
  return numpts;
}
 
inline void read_headers(FILE* ele_fd, FILE* fib_fd, bool read_binary, size_t & numelem, int & nFib)
{
  char buffer[HDR_SIZE];
  // first read number of elems
  if(read_binary) {
    fread(buffer, sizeof(char), HDR_SIZE, ele_fd);
  }
  else {
    fgets(buffer, HDR_SIZE, ele_fd);
  }
  sscanf(buffer, "%lu", &numelem);
 
  // now read number of fiber axes
  if(fib_fd) {
    if(read_binary) {
      fread(buffer, sizeof(char), HDR_SIZE, fib_fd);
      sscanf(buffer, "%d", &nFib);
    }
    else {
      // fiber files can have no header, thus we can potentially read 1, 3 or 6 values
      fgets(buffer, HDR_SIZE, fib_fd);
      float fbuff[6];
      short nread = sscanf(buffer, "%f %f %f %f %f %f", fbuff, fbuff+1, fbuff+2, fbuff+3, fbuff+4, fbuff+5);
      switch(nread) {
        case 1: nFib = fbuff[0]; break;
        case 3: nFib = 1;        break;
        case 6: nFib = 2;        break;
        default:
          fprintf(stderr, "%s error: read %d values in line 1 of fiber file."
                          "This is unexpected! Setting number of fibers to 1.\n",
                          __func__, nread);
          nFib = 1;
          break;
      }
 
      // if the fiber file had no header, we have to reset the FD
      if(nread != 1)
        fseek(fib_fd, 0, SEEK_SET);
    }
  }
  else
    nFib = 0;
}
inline void write_elem_headers(FILE* & ele_fd, FILE* & fib_fd, bool binary, size_t numelem, int nFib)
{
 
  // we need to write the CARP header
  char header[HDR_SIZE];
  int checksum = 666;
  
  // element header
  if(binary) {
    // first elements, second endianness, third checksum
    snprintf(header, sizeof header, "%lu %d %d", numelem, MT_ENDIANNESS, checksum);
    fwrite(header, HDR_SIZE, sizeof(char), ele_fd);
  }
  else
    fprintf(ele_fd, "%lu\n", numelem);
 
  // fiber header
  if(nFib >= 1) {
    if(binary) {
      // first numaxes, second numelems, third endianness, last checksum
      snprintf(header, sizeof header, "%d %lu %d %d", nFib, (unsigned long int)numelem, MT_ENDIANNESS, checksum);
      fwrite(header, sizeof(char), HDR_SIZE, fib_fd);
    }
    else
      fprintf(fib_fd, "%d\n", nFib);
  }
}
inline void write_pts_header(FILE* & pts_fd, bool binary, size_t numpts)
{
 
  // we need to write the CARP header
  char header[HDR_SIZE];
  int checksum = 666;
 
  // header
  if(binary) {
    snprintf(header, sizeof header, "%lu %d %d", numpts, MT_ENDIANNESS, checksum);
    fwrite(header, HDR_SIZE, sizeof(char), pts_fd);
  }
  else
    fprintf(pts_fd, "%lu\n", numpts);
}
 
template<class T, class S>
inline void read_elem_block(FILE* & fd, bool read_binary, size_t bstart, size_t bsize, meshdata<T, S> & mesh)
{
  const int bufsize = 2048;
  const int max_elem_size = 9;
 
  char buffer[bufsize];
  char etype_str[8];
  int n[max_elem_size];
 
  vector<T> & ref_eidx = mesh.register_numbering(NBR_ELEM_REF);
  ref_eidx.resize(bsize);
 
  mesh.dsp.resize(bsize+1);
  mesh.tag.resize(bsize);
  mesh.type.resize(bsize);
  mesh.con.resize(bsize*max_elem_size);
 
  size_t ele_read = 0, con_size = 0;
  mesh.dsp[0] = 0;
  T* elem = mesh.con.data();
 
  for(size_t i=0; i<bsize; i++)
  {
    if(read_binary) {
      int etbuff;
      size_t r = fread(&etbuff, sizeof(int), 1, fd);
      if(r != 1) break;
      mesh.type[i] = (elem_t)etbuff;
    }
    else {
      char* ptr = fgets( buffer, bufsize, fd);
      if(ptr == NULL) break;
      sscanf(buffer, "%s %d %d %d %d %d %d %d %d %d", etype_str, n, n+1, n+2, n+3, n+4, n+5, n+6, n+7, n+8);
      mesh.type[i] = getElemTypeID(etype_str);
    }
 
    T nodes;
    switch(mesh.type[i])
    {
      case Line:
        nodes = 2;
        break;
 
      case Tri:
        nodes = 3;
        break;
 
      case Quad:
        nodes = 4;
        break;
 
      case Tetra:
        nodes = 4;
        break;
 
      case Pyramid:
        nodes = 5;
        break;
 
      case Prism:
        nodes = 6;
        break;
 
      case Hexa:
        nodes = 8;
        break;
 
      default:
        fprintf(stderr, "Error: Unsupported element type!\n");
        exit(1);
    }
    mesh.dsp[i+1] = mesh.dsp[i] + nodes;
    con_size += nodes;
 
    // copy the element connectivity
    if(read_binary) {
      fread(n, sizeof(int), nodes+1, fd);
    }
    for(int j=0; j<nodes; j++) elem[j] = n[j];
 
    mesh.tag[i] = n[nodes];
    ref_eidx[i] = bstart + i;
    ele_read++;
 
    elem += nodes;
  }
  mesh.dsp.resize(ele_read+1);
  mesh.tag.resize(ele_read);
  mesh.type.resize(ele_read);
  ref_eidx.resize(ele_read);
  mesh.con.resize(con_size);
 
  mesh.l_numelem = ele_read;
}
 
 
template<class T, class S>
inline void write_elem_block(FILE* fd, bool write_binary, const meshdata<T, S> & mesh)
{
  const T* con = mesh.con.data();
 
  if(write_binary)
  {
    int wbuff[9];
 
    for(size_t eidx=0; eidx < mesh.l_numelem; eidx++)
    {
      fwrite(&mesh.type[eidx], 1, sizeof(elem_t), fd);
      T esize = mesh.dsp[eidx+1] - mesh.dsp[eidx];
      vec_assign(wbuff, con, esize); // copy-convert to int
      wbuff[esize] = mesh.tag[eidx]; // convert to int
 
      fwrite(wbuff, esize+1, sizeof(int), fd);
      con += esize;
    }
  }
  else
  {
    for(size_t eidx=0; eidx < mesh.l_numelem; eidx++)
    {
      switch(mesh.type[eidx])
      {
        case Line:
          fprintf(fd, "Ln %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 2;
          break;
 
        case Tri:
          fprintf(fd, "Tr %lld %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(con[2]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 3;
          break;
 
        case Quad:
          fprintf(fd, "Qd %lld %lld %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(con[2]),
                  static_cast<long long>(con[3]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 4;
          break;
 
        case Tetra:
          fprintf(fd, "Tt %lld %lld %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(con[2]),
                  static_cast<long long>(con[3]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 4;
          break;
 
        case Pyramid:
          fprintf(fd, "Py %lld %lld %lld %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(con[2]),
                  static_cast<long long>(con[3]),
                  static_cast<long long>(con[4]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 5;
          break;
 
        case Prism:
          fprintf(fd, "Pr %lld %lld %lld %lld %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(con[2]),
                  static_cast<long long>(con[3]),
                  static_cast<long long>(con[4]),
                  static_cast<long long>(con[5]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 6;
          break;
 
        case Hexa:
          fprintf(fd, "Hx %lld %lld %lld %lld %lld %lld %lld %lld %lld\n",
                  static_cast<long long>(con[0]),
                  static_cast<long long>(con[1]),
                  static_cast<long long>(con[2]),
                  static_cast<long long>(con[3]),
                  static_cast<long long>(con[4]),
                  static_cast<long long>(con[5]),
                  static_cast<long long>(con[6]),
                  static_cast<long long>(con[7]),
                  static_cast<long long>(mesh.tag[eidx]));
          con += 8;
          break;
 
        default:
          fprintf(stderr, "Error: Unsupported element type!\n");
          exit(1);
      }
    }
  }
}
 
 
template<class T, class S>
inline void read_fib_block(FILE* & fd, bool read_binary, int nFib, size_t bsize, meshdata<T, S> & mesh)
{
  const int bufsize = 2048;
  char buffer[bufsize];
 
  if(nFib == 2) {
    float fib[6];
    size_t nr = 0;
 
    mesh.fib.resize(bsize*3);
    mesh.she.resize(bsize*3);
 
    for(size_t i=0; i<bsize; i++) {
      if(read_binary) {
        size_t r = fread(fib, sizeof(float), 6, fd);
        if(r != 6) break;
      }
      else {
        char* ptr = fgets( buffer, bufsize, fd);
        if(ptr == NULL) break;
        sscanf(buffer, "%f %f %f %f %f %f", fib, fib+1, fib+2, fib+3, fib+4, fib+5);
      }
      mesh.fib[nr*3+0] = fib[0];
      mesh.fib[nr*3+1] = fib[1];
      mesh.fib[nr*3+2] = fib[2];
      mesh.she[nr*3+0] = fib[3];
      mesh.she[nr*3+1] = fib[4];
      mesh.she[nr*3+2] = fib[5];
 
      nr++;
    }
    mesh.fib.resize(nr*3);
    mesh.she.resize(nr*3);
  }
  else if (nFib == 1) {
    float fib[3];
    size_t nr = 0;
 
    mesh.fib.resize(bsize*3);
 
    for(size_t i=0; i<bsize; i++) {
      if(read_binary) {
        size_t r = fread(fib, sizeof(float), 3, fd);
        if(r != 3) break;
      }
      else {
        char* ptr = fgets( buffer, bufsize, fd);
        if(ptr == NULL) break;
        sscanf(buffer, "%f %f %f", fib, fib+1, fib+2);
      }      
      
      mesh.fib[nr*3+0] = fib[0];
      mesh.fib[nr*3+1] = fib[1];
      mesh.fib[nr*3+2] = fib[2];
 
      nr++;
    }
    mesh.fib.resize(nr*3);
  }
}
 
 
template<class T, class S>
inline void write_fib_block(FILE* & fd, bool write_binary, const meshdata<T, S> & mesh)
{
  float fib[6];
  int nFib = 0;
  if(mesh.fib.size() > 0) {
    nFib = 1;
    if(mesh.she.size() == mesh.fib.size())
      nFib = 2;
  }
 
  if(nFib == 2)
  {
    for(size_t i=0; i<mesh.l_numelem; i++)
    {
      fib[0] = mesh.fib[i*3+0];
      fib[1] = mesh.fib[i*3+1];
      fib[2] = mesh.fib[i*3+2];
      fib[3] = mesh.she[i*3+0];
      fib[4] = mesh.she[i*3+1];
      fib[5] = mesh.she[i*3+2];
 
      if(write_binary)
        fwrite(fib, 6, sizeof(float), fd);
      else
        fprintf(fd, "%f %f %f %f %f %f\n", fib[0], fib[1], fib[2], fib[3], fib[4], fib[5]);
    }
  }
  else if(nFib == 1)
  {
    for(size_t i=0; i<mesh.l_numelem; i++)
    {
      fib[0] = mesh.fib[i*3+0];
      fib[1] = mesh.fib[i*3+1];
      fib[2] = mesh.fib[i*3+2];
 
      if(write_binary)
        fwrite(fib, 3, sizeof(float), fd);
      else
        fprintf(fd, "%f %f %f\n", fib[0], fib[1], fib[2]);
    }
  }
}
 
template<class T, class S>
inline void read_element_tags(meshdata<T, S> & mesh, std::string filename)
{
  MPI_Comm comm = mesh.comm;
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  FILE* fd = NULL;
  size_t gnumelems_file = 0;
  int err = 0;
 
  if (rank == 0) {
    fd = fopen(filename.c_str(), "r");
    if (!fd) {
      fprintf(stderr, "Error: could not open tag file: %s. Aborting!\n", filename.c_str());
      err++;
    }
  }
  MPI_Allreduce(MPI_IN_PLACE, &err, 1, MPI_INT, MPI_SUM, comm);
  if (err) exit(EXIT_FAILURE);
 
  if (rank == 0) {
    char buffer[HDR_SIZE];
    if (!fgets(buffer, HDR_SIZE, fd) || sscanf(buffer, "%lu", &gnumelems_file) != 1) {
      fprintf(stderr, "Error: could not read element count from tag file: %s. Aborting!\n", filename.c_str());
      err++;
    }
  }
  MPI_Allreduce(MPI_IN_PLACE, &err, 1, MPI_INT, MPI_SUM, comm);
  if (err) {
    if (rank == 0) fclose(fd);
    exit(EXIT_FAILURE);
  }
 
  MPI_Bcast(&gnumelems_file, sizeof(size_t), MPI_BYTE, 0, comm);
 
  if (gnumelems_file != mesh.g_numelem) {
    if (rank == 0)
      fprintf(stderr, "Error: tag file %s declares %lu elements, mesh has %lu. Aborting!\n",
              filename.c_str(), gnumelems_file, mesh.g_numelem);
    if (rank == 0) fclose(fd);
    exit(EXIT_FAILURE);
  }
 
  size_t blocksize = (mesh.g_numelem + size - 1) / size;
 
  if (rank == 0) {
    for (size_t i = 0; i < mesh.l_numelem; i++) {
      int t;
      if (fscanf(fd, "%d", &t) != 1) {
        fprintf(stderr, "Error: unexpected end of tag file %s at element %lu. Aborting!\n",
                filename.c_str(), (unsigned long)i);
        fclose(fd);
        exit(EXIT_FAILURE);
      }
      mesh.tag[i] = static_cast<T>(t);
    }
 
    vector<T> buf;
    for (int pid = 1; pid < size; pid++) {
      size_t bsize = std::min(blocksize, mesh.g_numelem - (size_t)pid * blocksize);
      buf.resize(bsize);
      for (size_t i = 0; i < bsize; i++) {
        int t;
        if (fscanf(fd, "%d", &t) != 1) {
          fprintf(stderr, "Error: unexpected end of tag file %s at element %lu. Aborting!\n",
                  filename.c_str(), (unsigned long)((size_t)pid * blocksize + i));
          fclose(fd);
          exit(EXIT_FAILURE);
        }
        buf[i] = static_cast<T>(t);
      }
      MPI_Send(buf.data(), bsize * sizeof(T), MPI_BYTE, pid, SF_MPITAG, comm);
    }
    fclose(fd);
  } else {
    MPI_Status stat;
    MPI_Recv(mesh.tag.data(), mesh.l_numelem * sizeof(T), MPI_BYTE, 0, SF_MPITAG, comm, &stat);
  }
}
 
 
template<class T, class S>
inline void read_elements(meshdata<T, S> & mesh, std::string basename, bool require_fibers = true)
{
  // we use the communicator provided by the mesh
  MPI_Comm comm = mesh.comm;
 
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  FILE* ele_fd = NULL, *fib_fd = NULL;
  size_t gnumelems = 0;
  int nFib;
  bool read_binary = false;
  int err = 0;
 
  // open the files on root and read number of elements
  if(rank == 0) {
    read_binary = fileExists(basename + ".belem");
    std::string ele_file = read_binary ? basename + ".belem" : basename + ".elem";
    std::string fib_file = read_binary ? basename + ".blon"  : basename + ".lon";
 
    ele_fd = fopen(ele_file.c_str(), "r");
    if(!ele_fd) {
      fprintf(stderr, "Error: could not open file: %s. Aborting!\n", ele_file.c_str());
      err++;
    }
 
    fib_fd = fopen(fib_file.c_str(), "r");
    if((!fib_fd) && require_fibers) {
      fprintf(stderr, "Error: could not open file: %s. Aborting!\n", fib_file.c_str());
      err++;      
    }
  }
 
  MPI_Allreduce(MPI_IN_PLACE, &err, 1, MPI_INT, MPI_SUM, comm);
  if(err) {
    exit(EXIT_FAILURE);
  }
 
  if(rank == 0)
    read_headers(ele_fd, fib_fd, read_binary, gnumelems, nFib);
 
  MPI_Bcast(&gnumelems, sizeof(size_t), MPI_BYTE, 0, comm);
  MPI_Bcast(&nFib, sizeof(int), MPI_BYTE, 0, comm);
  mesh.g_numelem = gnumelems;
 
  // compute the size of the block stored on each process
  size_t blocksize = (gnumelems + size - 1) / size;
 
  // rank 0 reads mesh and distributes
  if(rank == 0) {
    // read own block
    read_elem_block(ele_fd, read_binary, 0, blocksize, mesh);
    if(nFib>=1)
      read_fib_block(fib_fd, read_binary, nFib, blocksize, mesh);
 
    // read blocks of other ranks and communicate
    meshdata<T, S> meshbuff;
    for(int pid=1; pid<size; pid++) {
      // read block in buffer
      read_elem_block(ele_fd, read_binary, pid*blocksize, blocksize, meshbuff);
      if(nFib>=1)
        read_fib_block(fib_fd, read_binary, nFib, blocksize, meshbuff);
      vector<T> & ref_eidx = meshbuff.get_numbering(NBR_ELEM_REF);
 
      // communicate
      MPI_Send(&meshbuff.l_numelem,     sizeof(size_t),                     MPI_BYTE, pid, SF_MPITAG, comm);
      MPI_Send(meshbuff.dsp.data(),     meshbuff.dsp.size()*sizeof(T),      MPI_BYTE, pid, SF_MPITAG, comm);
      MPI_Send(meshbuff.tag.data(),     meshbuff.tag.size()*sizeof(T),      MPI_BYTE, pid, SF_MPITAG, comm);
      MPI_Send(ref_eidx.data(),         ref_eidx.size()*sizeof(T),          MPI_BYTE, pid, SF_MPITAG, comm);
      MPI_Send(meshbuff.type.data(),    meshbuff.type.size()*sizeof(elem_t),MPI_BYTE, pid, SF_MPITAG, comm);
      if(nFib>=1)
        MPI_Send(meshbuff.fib.data(),     meshbuff.fib.size()*sizeof(S),      MPI_BYTE, pid, SF_MPITAG, comm);
      if(nFib==2)
        MPI_Send(meshbuff.she.data(),   meshbuff.she.size()*sizeof(S),      MPI_BYTE, pid, SF_MPITAG, comm);
 
      size_t con_size = meshbuff.con.size();
      MPI_Send(&con_size,             sizeof(size_t),     MPI_BYTE, pid, SF_MPITAG, comm);
      MPI_Send(meshbuff.con.data(), con_size*sizeof(T), MPI_BYTE, pid, SF_MPITAG, comm);
    }
 
    fclose(ele_fd);
    if(nFib>=1) fclose(fib_fd);
  }
  else {
    // the rest receives
    MPI_Status stat;
    MPI_Recv(&mesh.l_numelem, sizeof(size_t), MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    
    vector<T> & ref_eidx = mesh.register_numbering(NBR_ELEM_REF);
    ref_eidx.resize(mesh.l_numelem);
    
    mesh.dsp.resize(mesh.l_numelem+1);
    mesh.tag.resize(mesh.l_numelem);
    mesh.type.resize(mesh.l_numelem);
 
    if(nFib>=1) mesh.fib.resize(mesh.l_numelem*3);
    if(nFib==2) mesh.she.resize(mesh.l_numelem*3);
 
    MPI_Recv(mesh.dsp.data(),     mesh.dsp.size()*sizeof(T),      MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    MPI_Recv(mesh.tag.data(),     mesh.tag.size()*sizeof(T),      MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    MPI_Recv(ref_eidx.data(),     ref_eidx.size()*sizeof(T),      MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    MPI_Recv(mesh.type.data(),    mesh.type.size()*sizeof(elem_t),MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    if(nFib>=1)
      MPI_Recv(mesh.fib.data(),     mesh.fib.size()*sizeof(S),      MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    if(nFib==2)
      MPI_Recv(mesh.she.data(),   mesh.she.size()*sizeof(S),      MPI_BYTE, 0, SF_MPITAG, comm, &stat);
 
    size_t con_size;
    MPI_Recv(&con_size, sizeof(size_t), MPI_BYTE, 0, SF_MPITAG, comm, &stat);
    mesh.con.resize(con_size);
    MPI_Recv(mesh.con.data(), con_size*sizeof(T), MPI_BYTE, 0, SF_MPITAG, comm, &stat);
  }
 
  // localize connectivity (w.r.t. reference numbering)
  mesh.localize(NBR_REF);
}
 
 
template<class T, class S>
inline void write_elements(const meshdata<T, S> & mesh, bool binary, std::string basename)
{
  const MPI_Comm comm = mesh.comm;
 
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  FILE* ele_fd = NULL, *fib_fd = NULL;
  int nFib = 0;
  if(mesh.fib.size() > 0) {
    nFib = 1;
    if(mesh.she.size() == mesh.fib.size())
      nFib = 2;
  }  
  std::string ele_file = binary ? basename + ".belem" : basename + ".elem";
  std::string fib_file = binary ? basename + ".blon"  : basename + ".lon";
 
  // open the files on root and write headers
  if(rank == 0) {
    ele_fd = fopen(ele_file.c_str(), "w");
    if(!ele_fd) {
      fprintf(stderr, "Error: could not open file: %s. Aborting!\n", ele_file.c_str());
      exit(1);
    }
    if(nFib>=1) {
      fib_fd = fopen(fib_file.c_str(), "w");
      if(!fib_fd) {
        fprintf(stderr, "Error: could not open file: %s. Aborting!\n", fib_file.c_str());
        fclose(ele_fd);
        exit(1);
      }
    }
    write_elem_headers(ele_fd, fib_fd, binary, mesh.g_numelem, nFib);
    fclose(ele_fd);
    if (nFib>=1) fclose(fib_fd);
  }
 
  // write mesh sequentially rank after rank
  for(int pid=0; pid < size; pid++)
  {
    if(pid == rank)
    {
      ele_fd = fopen(ele_file.c_str(), "a");
      if(!ele_fd) {
        fprintf(stderr, "Error: could not open file: %s. Aborting!\n", ele_file.c_str());
        exit(1);
      }
      if(nFib >= 1) {
        fib_fd = fopen(fib_file.c_str(), "a");
        if(!fib_fd) {
          fprintf(stderr, "Error: could not open file: %s. Aborting!\n", fib_file.c_str());
          fclose(ele_fd);
          exit(1);
        }
      }
      write_elem_block(ele_fd, binary, mesh);
      if(nFib>=1) write_fib_block(fib_fd, binary, mesh);
 
      fclose(ele_fd);
      if(nFib>=1) fclose(fib_fd);
    }
    MPI_Barrier(comm);
  }
}
 
template<class T, class S>
inline void write_surface(const meshdata<T, S> & surfmesh, std::string surffile)
{
  const MPI_Comm comm = surfmesh.comm;
  const bool binary = false;
 
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  FILE* ele_fd = NULL;
 
  // open the files on root and write headers
  if(rank == 0) {
    ele_fd = fopen(surffile.c_str(), "w");
    if(!ele_fd) {
      fprintf(stderr, "Error: could not open file: %s. Aborting!\n", surffile.c_str());
      exit(1);
    }
    fprintf(ele_fd, "%lu\n", surfmesh.g_numelem);
    fclose(ele_fd);
  }
 
  // write surface sequentially rank after rank
  for(int pid=0; pid < size; pid++)
  {
    if(pid == rank) {
      ele_fd = fopen(surffile.c_str(), "a");
      if(!ele_fd) {
        fprintf(stderr, "Error: could not open file: %s. Aborting!\n", surffile.c_str());
        exit(1);
      }
 
      write_elem_block(ele_fd, binary, surfmesh);
      fclose(ele_fd);
    }
    MPI_Barrier(comm);
  }
}
 
 
template<class S>
inline void read_pts_block(FILE* & fd, bool read_binary, size_t bsize, vector<S> & xyz)
{
  const int bufsize = 2048;
  char buffer[bufsize];
  float pts[3];
  xyz.resize(bsize*3);
 
  size_t nr = 0;
 
  for(size_t i=0; i<bsize; i++) {
    if(read_binary) {
      size_t r = fread(pts, sizeof(float), 3, fd);
      if(r != 3) break;
    }
    else {
      char* ptr = fgets( buffer, bufsize, fd);
      if(ptr == NULL) break;
      sscanf(buffer, "%f %f %f", pts, pts+1, pts+2);
    }
    xyz[nr*3+0] = pts[0];
    xyz[nr*3+1] = pts[1];
    xyz[nr*3+2] = pts[2];
 
    nr++;
  }
  xyz.resize(nr*3);
}
 
 
template<class S>
inline void write_pts_block(FILE* & fd, bool write_binary, const vector<S> & xyz)
{
  size_t nnodes = xyz.size() / 3;
  float pt[3];
 
  for(size_t i=0; i<nnodes; i++)
  {
    pt[0] = xyz[i*3+0];
    pt[1] = xyz[i*3+1];
    pt[2] = xyz[i*3+2];
 
    if(write_binary)
      fwrite(pt, 3, sizeof(float), fd);
    else
      fprintf(fd, "%f %f %f\n", pt[0], pt[1], pt[2]);
  }
}
 
 
 
template<class T, class S>
inline void read_points(const std::string basename, const MPI_Comm comm, vector<S> & pts, vector<T> & ptsidx)
{
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  FILE* pts_fd = NULL;
  size_t gnumpts = 0;
  bool read_binary = false;
 
  // open the file on root and read number of points
  if(rank == 0) {
    read_binary = fileExists(basename + ".bpts");
    std::string pts_file = read_binary ? basename + ".bpts" : basename + ".pts";
 
    pts_fd = fopen(pts_file.c_str(), "r");
    if(!pts_fd) {
      fprintf(stderr, "Error: could not open file: %s. Aborting!\n", pts_file.c_str());
      exit(1);
    }
 
    // read number of points
    char buffer[2048];
    if(read_binary) fread(buffer, sizeof(char), HDR_SIZE, pts_fd);
    else fgets( buffer, 2048, pts_fd);
    sscanf(buffer, "%lu", &gnumpts);
  }
  MPI_Bcast(&gnumpts, sizeof(size_t), MPI_BYTE, 0, comm);
 
  // compute the size of the block to read .. this could also be set to a constant
  size_t blocksize = (gnumpts + size - 1) / size;
 
  if(rank == 0) {
    // first read own block
    read_pts_block(pts_fd, read_binary, blocksize, pts);
 
    // now read remaining blocks
    vector<S> buff;
    for(int pid = 1; pid < size; pid++) {
      read_pts_block(pts_fd, read_binary, blocksize, buff);
      long int numsend = buff.size();
 
      MPI_Send(&numsend, 1, MPI_LONG, pid, SF_MPITAG, comm);
      MPI_Send(buff.data(), numsend*sizeof(S), MPI_BYTE, pid, SF_MPITAG, comm);
    }
  }
  else {
    MPI_Status stat;
    long int numrecv = 0;
    MPI_Recv(&numrecv, 1, MPI_LONG, 0, SF_MPITAG, comm, &stat);
 
    pts.resize(numrecv);
    MPI_Recv(pts.data(), numrecv*sizeof(S), MPI_BYTE, 0, SF_MPITAG, comm, &stat);
  }
 
  long int mysize = pts.size() / 3;
  vector<long int> layout;
  layout_from_count(mysize, layout, comm);
 
  interval(ptsidx, layout[rank], layout[rank+1]);
 
  if(rank == 0) fclose(pts_fd);
}
 
template<class T, class S>
inline void insert_points(const vector<S> & pts, const vector<T> & ptsidx,
                          std::list<meshdata<T, S>*> & meshlist)
{
  // list must not be empty
  assert(meshlist.size() > 0);
  assert(pts.size() == (ptsidx.size() * 3));
 
  // we use the communicator provided by the first mesh in the list
  // the code does not support a meshlist with different communicators
  MPI_Comm comm = (*meshlist.begin())->comm;
 
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  vector<S> ptsbuff;
  vector<T> idxbuff;
 
  for(int pid = 0; pid < size; pid++) {
    size_t numsend = pts.size();
    MPI_Bcast(&numsend, sizeof(size_t), MPI_BYTE, pid, comm);
 
    if(rank == pid) {
      ptsbuff = pts;
      idxbuff = ptsidx;
    }
    else {
      ptsbuff.resize(numsend);
      idxbuff.resize(numsend / 3);
    }
 
    MPI_Bcast(ptsbuff.data(), ptsbuff.size()*sizeof(S), MPI_BYTE, pid, comm);
    MPI_Bcast(idxbuff.data(), idxbuff.size()*sizeof(T), MPI_BYTE, pid, comm);
 
    // iterate over all meshes in meshlist. if the interval of read indices contains indices of the mesh,
    // we copy the respective nodes from the point buffer into the mesh
    for(auto it = meshlist.begin(); it != meshlist.end(); ++it)
    {
      meshdata<T, S> & mesh = *(*it);  // reference to current mesh in meshlist
      const vector<T> & rnod = mesh.get_numbering(NBR_REF);
 
      // this function is used at the earliest stages of setting up a mesh,
      // thus we cannot yet use the parallel layout functionalities
      hashmap::unordered_map<T,T> g2l;
      for(size_t i=0; i<rnod.size(); i++)
        g2l[rnod[i]] = i;
 
      mesh.xyz.resize(mesh.l_numpts*3);
 
      vector<T> ridx, widx;
      ridx.reserve(rnod.size()), widx.reserve(rnod.size());
 
      for(size_t j=0; j<idxbuff.size(); j++) {
        if(g2l.count(idxbuff[j])) {
          ridx.push_back(j);
          widx.push_back(g2l[idxbuff[j]]);
        }
      }
 
      // insert those nodes which are in the local range
      for(size_t j=0; j<ridx.size(); j++) {
        T w = widx[j], r = ridx[j];
        mesh.xyz[w*3+0] = ptsbuff[r*3+0];
        mesh.xyz[w*3+1] = ptsbuff[r*3+1];
        mesh.xyz[w*3+2] = ptsbuff[r*3+2];
      }
    }
  }
}
 
template<class T, class S>
inline void writeVTKmesh_binary(const meshdata<T, S> & mesh, std::string file)
{
  FILE* vtk_file = fopen(file.c_str(), "w");
  if(vtk_file == NULL) return;
 
 
  fprintf (vtk_file, "# vtk DataFile Version 3.0\n");
  fprintf (vtk_file, "vtk output\n");
  fprintf (vtk_file, "binary\n");
  fprintf (vtk_file, "DATASET UNSTRUCTURED_GRID\n\n");
  fprintf (vtk_file, "POINTS %lu float\n", mesh.l_numpts);
 
  float pts[3];
  const S* p = mesh.xyz.data();
 
  for (unsigned long int i=0; i<mesh.l_numpts; i++ ) {
    pts[0] = htobe(p[0]);
    pts[1] = htobe(p[1]);
    pts[2] = htobe(p[2]);
    fwrite(pts, sizeof(float), 3, vtk_file);
    p += 3;
  }
  fprintf(vtk_file, "\n");
 
  fprintf (vtk_file, "CELL_TYPES %lu\n", mesh.l_numelem);
  unsigned long int valcount = 0;
  int vtk_type;
  for(unsigned long int i=0; i< mesh.l_numelem; i++) {
    elem_t etype = mesh.type[i];
    switch(etype) {
      // Line for Purkinje
      case Line:
        vtk_type =  3; // Lines are encoded as index 3
        break;
      case Tri:
        vtk_type = 5; // Triangles are encoded as index 5
        break;
      case Tetra:
        vtk_type = 10; // Tetras are encoded as index 10
        break;
      case Quad:
        vtk_type = 9; // Quads are encoded as index 9
        break;
      case Pyramid:
        vtk_type = 14; // Pyramids are encoded as index 14
        break;
      case Prism:
        vtk_type = 13; // Prisms are encoded as index 13
        break;
      case Hexa:
        vtk_type = 12; // Hexahedras are encoded as index 12
        break;
      default: break;
    }
    vtk_type = htobe(vtk_type);
    fwrite(&vtk_type, sizeof(int), 1, vtk_file);
    valcount += mesh.dsp[i+1] - mesh.dsp[i] + 1;
  }
  fprintf(vtk_file, "\n");
 
  fprintf(vtk_file, "CELLS %lu %lu\n", mesh.l_numelem, valcount);
  const T* elem = mesh.con.data();
  for(unsigned long int i=0; i<mesh.l_numelem; i++)
  {
    int nodes = mesh.dsp[i+1] - mesh.dsp[i], n;
    int be_nodes = htobe(nodes);
    fwrite(&be_nodes, sizeof(int), 1, vtk_file);
 
    for(int j=0; j<nodes; j++) {
      n = elem[j]; n = htobe(n);
      fwrite(&n, sizeof(int), 1, vtk_file);
    }
    elem += nodes;
  }
  fprintf(vtk_file, "\n");
 
  fprintf (vtk_file, "CELL_DATA %lu \n", mesh.l_numelem);
  fprintf (vtk_file, "SCALARS elemTag int 1\n");
  fprintf (vtk_file, "LOOKUP_TABLE default\n");
  for (unsigned long int i=0; i<mesh.l_numelem; i++ ) {
    int t = htobe(int(mesh.tag[i]));
    fwrite(&t, sizeof(int), 1, vtk_file);
  }
 
  // write fiber data
  int nFib = 0;
  if(mesh.fib.size() > 0) {
    nFib = 1;
    if(mesh.she.size() == mesh.fib.size())
      nFib = 2;
  }
  
  fprintf (vtk_file, "VECTORS fiber float\n");
  for (unsigned long int i=0; i<mesh.l_numelem; i++ ) {
    pts[0] = htobe(float(mesh.fib[i*3+0])),
    pts[1] = htobe(float(mesh.fib[i*3+1])),
    pts[2] = htobe(float(mesh.fib[i*3+2]));
    fwrite(pts, sizeof(float), 3, vtk_file);
  }
  fprintf(vtk_file, "\n");
 
  if(nFib == 2) {
    fprintf (vtk_file, "VECTORS sheet float\n");
    for (unsigned long int i=0; i<mesh.l_numelem; i++ ) {
      pts[0] = htobe(float(mesh.she[i*3+0])),
      pts[1] = htobe(float(mesh.she[i*3+1])),
      pts[2] = htobe(float(mesh.she[i*3+2]));
      fwrite(pts, sizeof(float), 3, vtk_file);
    }
    fprintf(vtk_file, "\n");
  }
 
  fclose(vtk_file);
}
 
}
#endif