MPIIO.cc

// Import necessary stuff
#include "MPIIO.h"
#include <cstdlib> // To get the exit function
#include <iostream>

/* -----------------------------------------------------------------------------
 Authors: Niels Aage, Erik Andreassen, Boyan Lazarov, August 2013
 Updated: June 2019, Niels Aage
 Copyright (C) 2013-2019,

 This MPIIO implementation is licensed under Version 2.1 of the GNU
 Lesser General Public License.

 This MMA implementation is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public
 License as published by the Free Software Foundation; either
 version 2.1 of the License, or (at your option) any later version.

 This Module is distributed in the hope that it will be useful,implementation
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public
 License along with this Module; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 -------------------------------------------------------------------------- */

/*
 * Modified by Zhidong Brian Zhang in May 2020, University of Waterloo
 */

// Constructor
#define _NO_SUCH_FILE 35

MPIIO::MPIIO (DM da_nodes, int nPf, std::string pnames, int nCf, std::string cnames) {

  // User defined string
  std::string infoString = "TopOpt result version 1.1";
  // Maximum number of points per element
#if DIM == 2  // # new
  int nPEl = 4; // 2D includes 4 nodes per element
#elif DIM == 3
  int nPEl = 8;
#endif
  // Number of domains (where domain refers to different regions in the
  // optimization)
  const int nDom = 1;

  // Number of point fields per domain:
  int nPFields[nDom] = { nPf };
  // The names of the point fields
  std::string pFieldNames = pnames;
  // Number of cell (element) fields per domain:
  int nCFields[nDom] = { nCf };
  // The names of the cell fields
  std::string cFieldNames = cnames;

  // Points/cells in each of the domains:
  nPointsMyrank = new unsigned long int[nDom];
  nCellsMyrank = new unsigned long int[nDom];

  // ----------------- Get the mesh: -------------------
  // Coordinates and a pointer
  Vec coordinates;
  PetscScalar *coordinatesPointer;

  // Get coordinates in local node numbering (including ghosts)
  DMGetCoordinatesLocal (da_nodes, &coordinates);
  VecGetArray (coordinates, &coordinatesPointer);

  // Get the global dof number (and ghosts)
  PetscInt nn;
  VecGetSize (coordinates, &nn);

  // Get the FE mesh structure (from the nodal mesh)
  PetscInt nel, nen;
  const PetscInt *necon;
#if DIM == 2 // # new
  DMDAGetElements_2D (da_nodes, &nel, &nen, &necon);
#elif DIM == 3
  DMDAGetElements_3D (da_nodes, &nel, &nen, &necon);
#endif

  // Number of points/cells in each domain for this rank
#if DIM == 2  // # new
  PetscInt numnodaldof = 2;
#elif DIM == 3
  PetscInt numnodaldof = 3;
#endif
  nPointsMyrank[0] = nn / numnodaldof;
  nCellsMyrank[0] = nel; // We have this number from when we called DMDAGetElements_2D/DMDAGetElements_3D

  // --------- Allocate the output object: -------
  Allocate (infoString, nDom, nPFields, nCFields, nPointsMyrank, nCellsMyrank, nPEl, pFieldNames, cFieldNames);

  // Write the points (or coordinates of the points)
  float *pointsDomain0 = new float[3 * nPointsMyrank[0]]; // always use "3" because point data in VTK All point data must use 3 coordinates
#if DIM == 2  // # new
  for (unsigned long int i = 0; i < nPointsMyrank[0]; i++) { // 2D
    // Convert to single precission to save space
    pointsDomain0[3 * i] = float (coordinatesPointer[DIM * i]);
    pointsDomain0[3 * i + 1] = float (coordinatesPointer[DIM * i + 1]);
    pointsDomain0[3 * i + 2] = float (0.0);
  }
#elif DIM == 3
  for (unsigned long int i = 0; i < 3 * nPointsMyrank[0]; i++) {
    // Convert to single precission to save space
    pointsDomain0[i] = float (coordinatesPointer[i]);
  }
#endif

  writePoints (0, pointsDomain0);
  // Restore coordinates array
  VecRestoreArray (coordinates, &coordinatesPointer);

  // Run through the elements of the domain (we have already called
  // DMDAGetElements)
  unsigned long int *cellsDomain0 = new unsigned long int[nPEl * nCellsMyrank[0]];
  unsigned long int *cellsOffset0 = new unsigned long int[nCellsMyrank[0]];
  unsigned long int *cellsTypes0 = new unsigned long int[nCellsMyrank[0]];
  unsigned long int CellOffset = 0;

  for (unsigned long int i = 0; i < nCellsMyrank[0]; i++) {
#if DIM == 2  // # new
    // Element type is the first number outputted:
    if (nen == 4) { // QUAD element
      cellsTypes0[i] = 9; // (in vtk QUAD element type is 9), http://victorsndvg.github.io/FEconv/formats/vtk.xhtml
    }
    // Then run through the nodes
    for (int j = 0; j < nen; j++) {
      cellsDomain0[i * nPEl + j] = necon[i * nen + j];
    }
    // Create the offset
    if (nen == 4) { // Quad element
      CellOffset += nen;
      cellsOffset0[i] = CellOffset; // (in vtk Quad element type is 9)
    }
    // Finally, in case we have varying elements size, make sure the extra nodes
    // are put to zero
    // REMARK: This is only an example, and is never used in this implementation
    for (int j = 4; j < nPEl; j++) {
      cellsDomain0[i * (nPEl + 1) + j + 1] = 0;
    }
#elif DIM == 3
    // Element type is the first number outputted:
    if (nen == 8) { // Hex element
      cellsTypes0[i] = 12; // (in vtk hex element type is 12)
    }
    // Then run through the nodes
    for (int j = 0; j < nen; j++) {
      cellsDomain0[i * nPEl + j] = necon[i * nen + j];
    }
    // Create the offset
    if (nen == 8) { // Hex element
      CellOffset += nen;
      cellsOffset0[i] = CellOffset; // (in vtk hex element type is 12)
    }
    // Finally, in case we have varying elements size, make sure the extra nodes
    // are put to zero
    // REMARK: This is only an example, and is never used in this implementation
    for (int j = 8; j < nPEl; j++) {
      cellsDomain0[i * (nPEl + 1) + j + 1] = 0;
    }
#endif
  }
  writeCells (0, cellsDomain0, cellsOffset0, cellsTypes0); // First domain

  // Allocate working arrays for outputting fields from timesteps:
  workPointField = new float[nPointsMyrank[0] * nPFields[0]]; // For first domain
  workCellField = new float[nCellsMyrank[0] * nCFields[0]]; // For first domain

  delete[] pointsDomain0;
  delete[] cellsDomain0;
  delete[] cellsOffset0;
  delete[] cellsTypes0;
}

// Destructor
MPIIO::~MPIIO () {
  // Delete the allocated arrays
  delete[] workPointField;
  delete[] workCellField;
  delete[] nPointsMyrank;
  delete[] nCellsMyrank;

  delete[] nPoints;
  delete[] nCells;
  delete[] nPointsT;
  delete[] nCellsT;
  delete[] nPFields;
  delete[] nCFields;
}

PetscErrorCode MPIIO::WriteVTK (DM da_nodes, Vec U, Vec nodeDensity, Vec x, Vec xTilde, Vec xPhys, Vec xPassive0, Vec xPassive1,Vec xPassive2, Vec xPassive3, PetscInt itr) {

  // Here we only have one "timestep" (no optimization)
  unsigned long int timestep = itr;

  PetscErrorCode ierr;

  // POINT FIELD(S)
  // Displacement
  Vec Ulocal;
  DMCreateLocalVector (da_nodes, &Ulocal);
  ierr = VecSet (Ulocal, 0.0);
  CHKERRQ(ierr);
  // Update the local vector from global solution
  ierr = DMGlobalToLocalBegin (da_nodes, U, INSERT_VALUES, Ulocal);
  CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd (da_nodes, U, INSERT_VALUES, Ulocal);
  CHKERRQ(ierr);
  // We need a pointer to the local vector
  PetscScalar *UlocalPointer;
  ierr = VecGetArray (Ulocal, &UlocalPointer);
  CHKERRQ(ierr);

  // Node density
  Vec NDlocal;
  DMCreateLocalVector (da_nodes, &NDlocal);
  ierr = VecSet (NDlocal, 0.0);
  CHKERRQ(ierr);
  // Update the local vector from global solution
  ierr = DMGlobalToLocalBegin (da_nodes, nodeDensity, INSERT_VALUES, NDlocal);
  CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd (da_nodes, nodeDensity, INSERT_VALUES, NDlocal);
  CHKERRQ(ierr);
  // We need a pointer to the local vector
  PetscScalar *NDlocalPointer;
  ierr = VecGetArray (NDlocal, &NDlocalPointer);
  CHKERRQ(ierr);

#if PHYSICS == 0 || PHYSICS == 1
  for (unsigned long int i = 0; i < nPointsMyrank[0]; i++) {
    // Ux
    workPointField[i] = float (UlocalPointer[DIM * i]);
    // Uy
    workPointField[i + nPointsMyrank[0]] = float (UlocalPointer[DIM * i + 1]);
#if DIM == 2   // # new
    // Uz fake, because all point data must use 3 coordinates in VTK files
    workPointField[i + 2 * nPointsMyrank[0]] = float (0.0);
#elif DIM == 3
    // Uz
    workPointField[i + 2 * nPointsMyrank[0]] = float (UlocalPointer[DIM * i + 2]);
#endif
    // Node density
    workPointField[i + 3 * nPointsMyrank[0]] = float (NDlocalPointer[DIM * i]);
  }
#elif PHYSICS == 2   // # new
  for (unsigned long int i = 0; i < nPointsMyrank[0]; i++) {
    // Ux
    workPointField[i] = float (UlocalPointer[i]);
    // Uy
    workPointField[i + nPointsMyrank[0]] = float (0.0);
#if DIM == 2
    // Uz fake, because all point data must use 3 coordinates in VTK files
    workPointField[i + 2 * nPointsMyrank[0]] = float (0.0);
#elif DIM == 3
    // Uz
    workPointField[i + 2 * nPointsMyrank[0]] = float (0.0);
#endif
    // Node density
    workPointField[i + 3 * nPointsMyrank[0]] = float (NDlocalPointer[i]);
  }
#endif
  writePointFields (timestep, 0, workPointField);
  // Restore Ulocal array
  ierr = VecRestoreArray (Ulocal, &UlocalPointer);
  CHKERRQ(ierr);
  ierr = VecRestoreArray (NDlocal, &NDlocalPointer);
  CHKERRQ(ierr);

  // CELL FIELD(S)
  PetscScalar *xpp, *xp, *xt;
  VecGetArray (x, &xp);
  VecGetArray (xTilde, &xt);
  VecGetArray (xPhys, &xpp);
  PetscScalar *xPassive0p, *xPassive1p, *xPassive2p, *xPassive3p;  // # new
  VecGetArray (xPassive0, &xPassive0p);  // # new
  VecGetArray (xPassive1, &xPassive1p);  // # new
  VecGetArray (xPassive2, &xPassive2p);  // # new
  VecGetArray (xPassive3, &xPassive3p);  // # new

  for (unsigned long int i = 0; i < nCellsMyrank[0]; i++) { // 2D/3D use the same code for cell field
    // Density
    workCellField[i + 0 * nCellsMyrank[0]] = float (xp[i]);
    workCellField[i + 1 * nCellsMyrank[0]] = float (xt[i]);
    workCellField[i + 2 * nCellsMyrank[0]] = float (xpp[i]);
    workCellField[i + 3 * nCellsMyrank[0]] = float (xPassive0p[i]);  // # new
    workCellField[i + 4 * nCellsMyrank[0]] = float (xPassive1p[i]);  // # new
    workCellField[i + 5 * nCellsMyrank[0]] = float (xPassive2p[i]);  // # new
    workCellField[i + 6 * nCellsMyrank[0]] = float (xPassive3p[i]);  // # new
  }
  writeCellFields (0, workCellField);

  // Restore arrays
  VecRestoreArray (x, &xp);
  VecRestoreArray (xTilde, &xt);
  VecRestoreArray (xPhys, &xpp);
  VecRestoreArray (xPassive0, &xPassive0p); // # new
  VecRestoreArray (xPassive1, &xPassive1p); // # new
  VecRestoreArray (xPassive2, &xPassive2p); // # new
  VecRestoreArray (xPassive3, &xPassive3p); // # new

  // clean up
  ierr = VecDestroy (&Ulocal);
  CHKERRQ(ierr);
  ierr = VecDestroy (&NDlocal);  // # new
  CHKERRQ(ierr);  // # new

  return ierr;
}

void MPIIO::Allocate (std::string info, const int nDom, const int nPFields[],
    const int nCFields[], unsigned long int nPointsMyrank[],
    unsigned long int nCellsMyrank[], unsigned long int nodesPerElement,
    std::string pFNames, std::string cFNames) //, std::string filename)
    /*  info = string with user defined info
     nDom = number of domains
     nPFields = array with number of point fields in each domain (the number
     we want to write) pFNames = string with point field names nCFields = array
     with number of cell fields in each domain (the number we want to write)
     cFNames = string with cell field names
     nPointsMyrank = array with number of points in each domain in
     thread=Myrank nCellsMyrank = array with number of cells in each domain in
     thread=Myrank nodesPerElement = (max) number of nodes per element. filename =
     name of output file (default = "output.dat")
     */
    {
  // default name
  std::string filename = "output_00000.dat";

  // Check PETSc input for a work directory
  char filenameChar[PETSC_MAX_PATH_LEN];
  PetscBool flg = PETSC_FALSE;
  PetscOptionsGetString (NULL, NULL, "-workdir", filenameChar,
      sizeof(filenameChar), &flg);

  // If input, change path of the file in filename
  if (flg) {
    filename = "";
    filename.append (filenameChar);
    filename.append ("/output.dat");
  }

  PetscPrintf (PETSC_COMM_WORLD,
      "##############################################"
          "##########################\n");
  PetscPrintf (MPI_COMM_WORLD, "Outputfile is written to: %s \n",
      filename.c_str ());
  PetscPrintf (MPI_COMM_WORLD,
      "To change the working directory, specify '-workdir' at runtime\n");

  // Continue to allocate
  firstFieldOutputDone = false; // Initialize
  this->filename = filename;
  int ierror;
  int headerLen;
  MPI_File fh; // File handle
  // Find how many bytes are used to store an unsigned long integer
  MPI_Type_size (MPI_UNSIGNED_LONG, &MPI_IS);
  // Bytes used to store a float
  MPI_Type_size (MPI_FLOAT, &MPI_FS);
  // Bytes used to store a char
  MPI_Type_size (MPI_CHAR, &MPI_CS);
  // Find out how many cpus we have and their ranks
  ierror = MPI_Comm_rank (MPI_COMM_WORLD, &rank);
  if (ierror) {
    abort ("Problems getting rank", "MPIIO:MPIIO");
  }
  ierror = MPI_Comm_size (MPI_COMM_WORLD, &ncpu);
  if (ierror) {
    abort ("Problems getting number of cpus", "MPIIO:MPIIO");
  }

  // Communicate number of points and cells to all processors
  // First allocate space for the arrays
  // Check how many domains your trying to allocate (to avoid crash)
  if (nDom > 1000) {
    abort ("ERROR: More than 1000 domains!", "MPIIO:MPIIO");
  }
  nPoints = new unsigned long int[nDom * ncpu];
  nCells = new unsigned long int[nDom * ncpu];
  for (int i = 0; i < nDom; i++) {
    ierror = MPI_Allgather(&nPointsMyrank[i], 1, MPI_UNSIGNED_LONG,
        &nPoints[i * ncpu], 1, MPI_UNSIGNED_LONG, MPI_COMM_WORLD);
    if (ierror) {
      abort ("Problems exchanging number of points", "MPIIO:MPIIO");
    }
    ierror = MPI_Allgather(&nCellsMyrank[i], 1, MPI_UNSIGNED_LONG,
        &nCells[i * ncpu], 1, MPI_UNSIGNED_LONG, MPI_COMM_WORLD);
    if (ierror) {
      abort ("Problems exchanging number of cells", "MPIIO:MPIIO");
    }
  }
  // Also allocate and save the other data provided
  this->nDom = nDom;
  this->nodesPerElement = nodesPerElement;
  this->nPFields = new int[nDom];
  this->nCFields = new int[nDom];
  this->nPointsT = new unsigned long int[nDom]; // Total number of points
  this->nCellsT = new unsigned long int[nDom]; // Total number of cells
  // All processors position in the file is moved below the outputted data
  headerLen = 2 + 4 * nDom;
  unsigned long int *header = new unsigned long int[headerLen];
  // Put the number of domains into the buffer
  header[0] = nDom;
  for (int i = 0; i < nDom; i++) {
    // Sum up total number of points and cells in the domain
    header[1 + i] = 0;
    header[1 + nDom + i] = 0;
    for (int j = 0; j < ncpu; j++) {
      header[1 + i] += nPoints[i * ncpu + j];
      header[1 + nDom + i] += nCells[i * ncpu + j];
    }
    nPointsT[i] = header[1 + i];
    nCellsT[i] = header[1 + nDom + i];
    // And then the number of point/cell fields in each domain
    header[1 + 2 * nDom + i] = nPFields[i];
    header[1 + 3 * nDom + i] = nCFields[i];
    this->nPFields[i] = nPFields[i];
    this->nCFields[i] = nCFields[i];
  }
  // Last entry is the nodes per element
  header[headerLen - 1] = nodesPerElement;
  // Save the number of characters to output
  int numberOfCharacters = info.size () + pFNames.size () + cFNames.size () + 4;
  // The first processor outputs total number of points, cells, and fields
  if (rank == 0) { // Can this part be done as standard C++ binary output?
    // If there is an old file, delete this
    ierror = MPI_File_delete (&filename[0], MPI_INFO_NULL);
    // The below line is commented since it caused errors with some MPI
    // implementations
    // if (ierror != _NO_SUCH_FILE && ierror) {abort("Problems deleting old
    // file", "MPIIO:MPIIO");}
    // Then open file
    ierror = MPI_File_open (MPI_COMM_SELF, &filename[0],
    MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
    if (ierror) {
      abort ("Problems opening file", "MPIIO::MPIIO");
    }
    // No need to create filetype here, its just a buffer with integers
    offset = 0; // Start at the beginning of the file
    // Set view
    ierror = MPI_File_set_view (fh, offset, MPI_CHAR, MPI_CHAR,
        (char*) "native", MPI_INFO_NULL);
    if (ierror) {
      abort ("Problems setting view", "MPIIO::MPIIO");
    }
    // Write to the file
    info.append ("\n\x01"); // Make sure the string ends with an endline
    ierror = MPI_File_write (fh, (char*) info.c_str (), info.size (), MPI_CHAR,
    MPI_STATUS_IGNORE);
    if (ierror) {
      abort ("Problems writing to file", "MPIIO::MPIIO");
    }
    // Set view
    offset += MPI_CS * info.size (); // Adjust offset
    ierror = MPI_File_set_view (fh, offset, MPI_UNSIGNED_LONG,
    MPI_UNSIGNED_LONG, (char*) "native", MPI_INFO_NULL);
    if (ierror) {
      abort ("Problems setting view", "MPIIO::MPIIO");
    }
    // Write to the file
    ierror = MPI_File_write (fh, header, headerLen, MPI_UNSIGNED_LONG,
    MPI_STATUS_IGNORE);
    if (ierror) {
      abort ("Problems writing to file", "MPIIO::MPIIO");
    }
    // Set view
    offset += MPI_IS * headerLen; // Adjust offset
    ierror = MPI_File_set_view (fh, offset, MPI_CHAR, MPI_CHAR,
        (char*) "native", MPI_INFO_NULL);
    if (ierror) {
      abort ("Problems setting view", "MPIIO::MPIIO");
    }
    // Write to the file
    pFNames.append ("\x01"); // Make sure the string ends with an endline
    cFNames.append ("\x01"); // Make sure the string ends with an endline
    pFNames.append (cFNames); // Output both strings at once
    // Write to the file
    ierror = MPI_File_write (fh, (char*) pFNames.c_str (), pFNames.size (),
    MPI_CHAR, MPI_STATUS_IGNORE);
    // Close the file (I don't think we need a barrier here)
    ierror = MPI_File_close (&fh);
    if (ierror) {
      abort ("Problems closing file", "MPIIO::MPIIO");
    }
  }
  MPI_Barrier (MPI_COMM_WORLD);
  // All processors position in the file is moved below the outputted data
  offset = MPI_IS * headerLen + MPI_CS * numberOfCharacters;
  // ALWAYS remember to deallocate:
  delete[] header;
}

// Output coordinates - only done once
void MPIIO::writePoints (int domain, float coordinates[])
    /*
     domain = The domain number (not the name)
     coordinates = An array with the coordinates
     */
    {
  int ierror;
  // Open file
  ierror = MPI_File_open (MPI_COMM_WORLD, &filename[0],
  MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
  if (ierror) {
    abort ("Problems opening file", "MPIIO::writePoints");
  }
  // Compute offset for the different processors
  if (domain == 0) { // For the first domain sum up to the current position
    offset += 3 * sum (nPoints, rank) * MPI_FS;
  } else { // For the rest of the domains add the points written since last time
    offset += 3 * sum (&nPoints[ncpu * (domain - 1) + rank], ncpu) * MPI_FS;
  }
  // Set view
  ierror = MPI_File_set_view (fh, offset, MPI_FLOAT, MPI_FLOAT,
      (char*) "native", MPI_INFO_NULL);
  if (ierror) {
    abort ("Problems setting view", "MPIIO::writePoints");
  }
  // Write to the file
  int len = 3 * nPoints[domain * ncpu + rank]; // Number of floats to write
  ierror = MPI_File_write_all (fh, coordinates, len, MPI_FLOAT,
  MPI_STATUS_IGNORE);
  if (ierror) {
    abort ("Problems writing to file", "MPIIO::writePoints");
  }
  // Close the file (I don't think we need a barrier here)
  ierror = MPI_File_close (&fh);
  if (ierror) {
    abort ("Problems closing file", "MPIIO::writePoints");
  }
}

// Output cells - only done once
// You provide the elements with "local numbering" -
// the method will convert to "global numbering".
void MPIIO::writeCells (int domain, unsigned long int elements[],
    unsigned long int cellsOffset0[], unsigned long int cellsTypes0[])
    /*
     domain = The domain number
     elements = Array with node connectivities. It is assumed that
     the node numbering is local in both domain and thread.
     And that we want a an overall global numbering for all
     domains.
     Furthermore, every nodesPerElement number should be
     an element type number.
     */
    {
  int ierror;
  unsigned long int shift;
  // Compute the shift number (from local to global node number)
  // This is done by summing up all points outputted before the points in this
  // domain from this rank
  shift = sum (nPoints, ncpu * domain + rank);
  // Run through all "elements" and make them global by adding "shift":
  for (unsigned long int i = 0;
      i < (nodesPerElement) * nCells[ncpu * domain + rank]; i++) {
    // but shift all the nodes to global numbering
    elements[i] += shift;
  }
  // Open file
  ierror = MPI_File_open (MPI_COMM_WORLD, &filename[0],
  MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
  if (ierror) {
    abort ("Problems opening file", "MPIIO::writeCells");
  }

  // Compute offset FOR ELEMENT CONN for the different processors
  if (domain == 0) { // For the first domain sum up to the current position
    // First, add for the remaining points written
    offset += 3 * sum (&nPoints[ncpu * (nDom - 1) + rank], ncpu - rank)
              * MPI_FS;
    // Then, add the cells
    offset += (nodesPerElement) * sum (nCells, rank) * MPI_IS;
  } else { // For the rest of the domains add the elements written since last
           // time
    offset += (nodesPerElement)
              * sum (&nCells[ncpu * (domain - 1) + rank], ncpu)
              * MPI_IS;
  }
  // Set view
  ierror = MPI_File_set_view (fh, offset, MPI_UNSIGNED_LONG, MPI_UNSIGNED_LONG,
      (char*) "native", MPI_INFO_NULL);
  if (ierror) {
    abort ("Problems setting view", "MPIIO::writeCells");
  }
  // Length of data stream to write
  unsigned long int len = (nodesPerElement) * nCells[domain * ncpu + rank]; // Number of integers
  // Write ELEMENT Conn to file
  ierror = MPI_File_write_all (fh, elements, len, MPI_UNSIGNED_LONG,
  MPI_STATUS_IGNORE);
  if (ierror) {
    abort ("Problems writing ELEMENTS to file", "MPIIO::writeCells");
  }

  // Write the VTK OFFSET
  // Update the write offset
  // First jump past ALL the connectivity
  offset += (nodesPerElement)
            * sum (&nCells[ncpu * (nDom - 1) + rank], ncpu - rank)
            * MPI_IS;
  // Next jump past the previous ranks offset list
  offset += sum (nCells, rank) * MPI_IS;
  unsigned long int addToOffsetList = nodesPerElement * sum (nCells, rank);
  for (int i = 0; i < (int) nCells[ncpu * domain + rank]; i++) {
    cellsOffset0[i] += addToOffsetList;
  }
  // Length of the offset to write
  len = nCells[domain * ncpu + rank]; // Number of integers
  // Move the offset in the file
  ierror = MPI_File_set_view (fh, offset, MPI_UNSIGNED_LONG, MPI_UNSIGNED_LONG,
      (char*) "native", MPI_INFO_NULL);
  if (ierror) {
    abort ("Problems setting view OFFSET", "MPIIO::writeCells");
  }
  // write the offset list
  ierror = MPI_File_write_all (fh, cellsOffset0, len, MPI_UNSIGNED_LONG,
  MPI_STATUS_IGNORE);

  // Write the VTK ELEMENT TYPE
  // First jump past ALL the offsets
  offset += sum (&nCells[ncpu * (nDom - 1) + rank], ncpu - rank) * MPI_IS;
  // Nextjump past the previous ranks type list
  offset += sum (nCells, rank) * MPI_IS;
  // Length of type list for this rank
  len = nCells[domain * ncpu + rank]; // Number of integers
                                      // Move the offset in the file
  ierror = MPI_File_set_view (fh, offset, MPI_UNSIGNED_LONG, MPI_UNSIGNED_LONG,
      (char*) "native", MPI_INFO_NULL);
  // write the type list to file
  ierror = MPI_File_write_all (fh, cellsTypes0, len, MPI_UNSIGNED_LONG,
  MPI_STATUS_IGNORE);

  // Close the file (I don't think we need a barrier here)
  ierror = MPI_File_close (&fh);
  if (ierror) {
    abort ("Problems closing file", "MPIIO::writeCells");
  }
}

// Output point fields
void MPIIO::writePointFields (unsigned long int timeStep, int domain,
    float fields[], std::string newFilename)
    // timeStep = integer specifying time step
// domain = domain number (integer, should always start at zero and increase by
// one) fields = array containing field values (for all fields, the fields come
// after each other),
//          which means the user has to put the field values in an array
// before passing it to             this method. This is done because, the entries have to be
// converted from double            to single precision anyway, and then it does not matter
// that the user            has to allocate an extra array and store the single precision
// values in this.          Furthermore, it simplifies the MPI_IO commands.
// newFilename = optional argument with filename to the file you want to write
// NB: It is assumed that all domains are written to the same file
//     and that they are always written in the same chronological order (this
//     could be changed)
    {
  int ierror;
  if (newFilename != "notDefined" && newFilename != filename) {
    if (domain != 0) {
      abort ("Only new filename when first domain!", "MPIIO::writePointFields");
    }
    filename = newFilename;
    // Reset positions
    offset = 0;
  }
  // Compute the offset
  else if (domain == 0) {
    if (!firstFieldOutputDone) {
      // Add for the remaining elements written
      offset += sum (&nCells[ncpu * (nDom - 1) + rank], ncpu - rank) * MPI_IS;
    } else { // We have outputted fields earlier
      // Add for the remaining cell field values written
      offset += sum (&nCells[ncpu * (nDom - 1) + rank], ncpu - rank) * MPI_FS;
    }
  }
  if (domain == 0) { // For the first domain sum up to the current position
    // Add the point field values written by other ranks
    offset += sum (nPoints, rank) * MPI_FS;
  } else { // For the rest of the domains add the point field values written
           // since last time
    // From the previous domain
    offset += sum (&nPoints[ncpu * (domain - 1) + rank], ncpu - rank) * MPI_FS;
    // From the current domain (if rank != 0)
    offset += sum (&nPoints[ncpu * domain], rank) * MPI_FS;
  }
  // If it is the first domain (domain=0), the timeStep should be written by
  // rank 0 and all offsets should be increased by one
  if (domain == 0) {
    if (rank == 0) {
      // abort(" --------- HERE---------------", "MPIIO::writePointFields");
      ierror = MPI_File_open (MPI_COMM_SELF, &filename[0],
      MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
      if (ierror) {
        abort ("Problems opening file", "MPIIO::writePointFields");
      }
      // Set view
      ierror = MPI_File_set_view (fh, offset, MPI_UNSIGNED_LONG,
      MPI_UNSIGNED_LONG, (char*) "native", MPI_INFO_NULL);
      if (ierror) {
        abort ("Problems setting view", "MPIIO::writePointFields");
      }
      // Write to the file
      ierror = MPI_File_write (fh, &timeStep, 1, MPI_UNSIGNED_LONG,
      MPI_STATUS_IGNORE);
      if (ierror) {
        abort ("Problems writing to file", "MPIIO::writePointFields");
      }
      // Close the file
      ierror = MPI_File_close (&fh);
      if (ierror) {
        abort ("Problems closing file", "MPIIO::writePointFields");
      }
    }
    MPI_Barrier (MPI_COMM_WORLD);
    // Increase all offsets by one integer:
    offset += 1 * MPI_IS;
  }
  // MPI_Barrier(MPI_COMM_WORLD);
  // Open the file
  ierror = MPI_File_open (MPI_COMM_WORLD, &filename[0],
  MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
  if (ierror) {
    abort ("Problems opening file", "MPIIO::writePointFields");
  }
  // Set filetype
  MPI_Datatype filetype;
  // The number of points for the specified rank in this domain
  int blocklength = nPoints[ncpu * domain + rank];
  // The total number of points in the domain
  int stride = nPointsT[domain];
  // Number of blocks to write
  int count = nPFields[domain];
  ierror = MPI_Type_vector (count, blocklength, stride, MPI_FLOAT, &filetype);
  if (ierror) {
    abort ("Problems creating MPI vector", "MPIIO::writePointFields");
  }
  ierror = MPI_Type_commit (&filetype);
  if (ierror) {
    abort ("Problems creating filetype", "MPIIO::writePointFields");
  }
  // Set view
  ierror = MPI_File_set_view (fh, offset, MPI_FLOAT, filetype, (char*) "native",
  MPI_INFO_NULL);
  if (ierror) {
    abort ("Problems setting view", "MPIIO::writePointFields");
  }

  // Since the array fields contain all fields, we can output them all at once
  // Remember: datatype specifies the layout in memory, while filetype specifies
  // the layout in the file. They are both of the MPI_Datatype kind.
  // Set datatype such that the access in memory is correct. In this case the
  // field values are in the right order already, so no need to make a datatype.
  // Write to the file (one datatype is written)
  ierror = MPI_File_write_all (fh, fields, count * blocklength, MPI_FLOAT,
  MPI_STATUS_IGNORE);
  if (ierror) {
    abort ("Problems writing to file", "MPIIO::writePointFields");
  }
  // Close the file
  ierror = MPI_File_close (&fh);
  if (ierror) {
    abort ("Problems closing file", "MPIIO::writePointFields");
  }
  // Check if it was the first time this function has been called
  if (!firstFieldOutputDone) {
    firstFieldOutputDone = true;
  }
  // Free the memory used for filetype
  ierror = MPI_Type_free (&filetype);
  if (ierror) {
    abort ("Problems freeing datatype", "MPIIO::writePointFields");
  }
  // Finally, update the offset to the beginning of the last field we wrote
  offset += stride * (count - 1) * MPI_FS;
}

// Output cell fields
void MPIIO::writeCellFields (int domain, float fields[])
    // timeStep = integer specifying time step
// domain = domain number
// fields = field values at cell points. Should contain all fields!
// NB: It is always assumed that cell fields are written to the same file as
// point fields
//     and that all domains are written to the same file, and that point fields
//     are called first
    {
  int ierror;
  if (domain == 0) {
    // First, add for the remaining points written
    offset += sum (&nPoints[ncpu * (nDom - 1) + rank], ncpu - rank) * MPI_FS;
    // Add the cell field values written by other ranks
    offset += sum (nCells, rank) * MPI_FS;
  } else { // For the rest of the domains add the cell field values written
           // since last time
    // From the previous domain
    offset += sum (&nCells[ncpu * (domain - 1) + rank], ncpu - rank) * MPI_FS;
    // From the current domain (if rank != 0)
    offset += sum (&nCells[ncpu * domain], rank) * MPI_FS;
  }
  // Open the file
  ierror = MPI_File_open (MPI_COMM_WORLD, &filename[0],
  MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &fh);
  if (ierror) {
    abort ("Problems opening file", "MPIIO::writeCellFields");
  }

  // Set filetype
  MPI_Datatype filetype;
  // The number of points for the specified rank in this domain
  int blocklength = nCells[ncpu * domain + rank];
  // The total number of points in the domain
  unsigned long int stride = nCellsT[domain];
  // Number of blocks to write
  int count = nCFields[domain];
  ierror = MPI_Type_vector (count, blocklength, stride, MPI_FLOAT, &filetype);
  if (ierror) {
    abort ("Problems creating MPI vector", "MPIIO::writeCellFields");
  }
  ierror = MPI_Type_commit (&filetype);
  if (ierror) {
    abort ("Problems creating filetype", "MPIIO::writeCellFields");
  }
  // Set view
  ierror = MPI_File_set_view (fh, offset, MPI_FLOAT, filetype, (char*) "native",
  MPI_INFO_NULL);
  if (ierror) {
    abort ("Problems setting view", "MPIIO::writeCellFields");
  }

  // Since the array fields contain all fields, we can output them all at once
  // Remember: datatype specifies the layout in memory, while filetype specifies
  // the layout in the file. They are both of the MPI_Datatype kind.
  // Set datatype such that the access in memory is correct. In this case the
  // field values are in the right order already, so no need to make a datatype.
  // Write to the file (one datatype is written)
  ierror = MPI_File_write_all (fh, fields, count * blocklength, MPI_FLOAT,
  MPI_STATUS_IGNORE);
  if (ierror) {
    abort ("Problems writing to file", "MPIIO::writeCellFields");
  }
  // Close the file
  ierror = MPI_File_close (&fh);
  if (ierror) {
    abort ("Problems closing file", "MPIIO::writeCellFields");
  }
  // Free the memory used for filetype
  ierror = MPI_Type_free (&filetype);
  if (ierror) {
    abort ("Problems freeing datatype", "MPIIO::writeCellFields");
  }
  // Finally, update the offset to the beginning of the last field we wrote
  offset += stride * (count - 1) * MPI_FS;
}

// Method to do MPI errors
void MPIIO::abort (std::string errorMsg, std::string position)
    // errorMsg = the error message the programmer has written
// position = the methode in which the error occured
    {
  std::cerr << errorMsg << " in " << position << std::endl;
  // Stop the execution of the program
  MPI_Barrier (MPI_COMM_WORLD);
  std::cerr << "rank = " << rank << std::endl;
  MPI_Abort (MPI_COMM_WORLD, -1);
  exit (0);
}

unsigned long int MPIIO::sum (unsigned long int *startPos,
    unsigned long int nel) {
  unsigned long int total = 0; // The number of points
  for (unsigned long int i = 0; i < nel; i++) {
    total += startPos[i];
  }
  return total;
}

#if DIM == 2   // # new
PetscErrorCode MPIIO::DMDAGetElements_2D (DM dm, PetscInt *nel, PetscInt *nen,
    const PetscInt *e[]) {
  PetscErrorCode ierr;
  DM_DA *da = (DM_DA*) dm->data;
  PetscInt i, xs, xe, Xs, Xe;
  PetscInt j, ys, ye, Ys, Ye;
  PetscInt cnt = 0, cell[4], ns = 1, nn = 4;
  PetscInt c;
  if (!da->e) {
    if (da->elementtype == DMDA_ELEMENT_Q1) {
      ns = 1;
      nn = 4;
    }
    ierr = DMDAGetCorners (dm, &xs, &ys, NULL, &xe, &ye, NULL);
    CHKERRQ(ierr);
    ierr = DMDAGetGhostCorners (dm, &Xs, &Ys, NULL, &Xe, &Ye, NULL);
    CHKERRQ(ierr);
    xe += xs;
    Xe += Xs;
    if (xs != Xs) xs -= 1;
    ye += ys;
    Ye += Ys;
    if (ys != Ys) ys -= 1;
    da->ne = ns * (xe - xs - 1) * (ye - ys - 1);
    PetscMalloc((1 + nn * da->ne) * sizeof(PetscInt), &da->e);
    for (j = ys; j < ye - 1; j++) {
      for (i = xs; i < xe - 1; i++) {
        cell[0] = (i - Xs) + (j - Ys) * (Xe - Xs);
        cell[1] = (i - Xs + 1) + (j - Ys) * (Xe - Xs);
        cell[2] = (i - Xs + 1) + (j - Ys + 1) * (Xe - Xs);
        cell[3] = (i - Xs) + (j - Ys + 1) * (Xe - Xs);
        if (da->elementtype == DMDA_ELEMENT_Q1) {
          for (c = 0; c < ns * nn; c++)
            da->e[cnt++] = cell[c];
        }
      }
    }
  }
  *nel = da->ne;
  *nen = nn;
  *e = da->e;
  return (0);
}

#elif DIM == 3
PetscErrorCode MPIIO::DMDAGetElements_3D (DM dm, PetscInt *nel, PetscInt *nen,
    const PetscInt *e[]) {
  PetscErrorCode ierr;
  DM_DA *da = (DM_DA*) dm->data;
  PetscInt i, xs, xe, Xs, Xe;
  PetscInt j, ys, ye, Ys, Ye;
  PetscInt k, zs, ze, Zs, Ze;
  PetscInt cnt = 0, cell[8], ns = 1, nn = 8;
  PetscInt c;
  if (!da->e) {
    if (da->elementtype == DMDA_ELEMENT_Q1) {
      ns = 1;
      nn = 8;
    }
    ierr = DMDAGetCorners (dm, &xs, &ys, &zs, &xe, &ye, &ze);
    CHKERRQ(ierr);
    ierr = DMDAGetGhostCorners (dm, &Xs, &Ys, &Zs, &Xe, &Ye, &Ze);
    CHKERRQ(ierr);
    xe += xs;
    Xe += Xs;
    if (xs != Xs) xs -= 1;
    ye += ys;
    Ye += Ys;
    if (ys != Ys) ys -= 1;
    ze += zs;
    Ze += Zs;
    if (zs != Zs) zs -= 1;
    da->ne = ns * (xe - xs - 1) * (ye - ys - 1) * (ze - zs - 1);
    PetscMalloc((1 + nn * da->ne) * sizeof(PetscInt), &da->e);
    for (k = zs; k < ze - 1; k++) {
      for (j = ys; j < ye - 1; j++) {
        for (i = xs; i < xe - 1; i++) {
          cell[0] = (i - Xs) + (j - Ys) * (Xe - Xs)
                    + (k - Zs) * (Xe - Xs) * (Ye - Ys);
          cell[1] = (i - Xs + 1) + (j - Ys) * (Xe - Xs)
                    + (k - Zs) * (Xe - Xs) * (Ye - Ys);
          cell[2] = (i - Xs + 1) + (j - Ys + 1) * (Xe - Xs)
                    + (k - Zs) * (Xe - Xs) * (Ye - Ys);
          cell[3] = (i - Xs) + (j - Ys + 1) * (Xe - Xs)
                    + (k - Zs) * (Xe - Xs) * (Ye - Ys);
          cell[4] = (i - Xs) + (j - Ys) * (Xe - Xs)
                    + (k - Zs + 1) * (Xe - Xs) * (Ye - Ys);
          cell[5] = (i - Xs + 1) + (j - Ys) * (Xe - Xs)
                    + (k - Zs + 1) * (Xe - Xs) * (Ye - Ys);
          cell[6] = (i - Xs + 1) + (j - Ys + 1) * (Xe - Xs)
                    + (k - Zs + 1) * (Xe - Xs) * (Ye - Ys);
          cell[7] = (i - Xs) + (j - Ys + 1) * (Xe - Xs)
                    + (k - Zs + 1) * (Xe - Xs) * (Ye - Ys);
          if (da->elementtype == DMDA_ELEMENT_Q1) {
            for (c = 0; c < ns * nn; c++)
              da->e[cnt++] = cell[c];
          }
        }
      }
    }
  }
  *nel = da->ne;
  *nen = nn;
  *e = da->e;
  return (0);
}
#endif