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
-------------------------------------------------------------------------- */

// 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
    int nPEl = 8;
    // 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;
    DMDAGetElements_3D(da_nodes, &nel, &nen, &necon);

    // Number of points/cells in each domain for this rank
    PetscInt numnodaldof = 3;
    nPointsMyrank[0]     = nn / numnodaldof;
    nCellsMyrank[0]      = nel; // We have this number from when we called DMDAGetElements_3D

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

    // Write the points (or coordinates of the points)
    float* pointsDomain0 = new float[3 * nPointsMyrank[0]];
    for (unsigned long int i = 0; i < 3 * nPointsMyrank[0]; i++) {
        // Convert to single precission to save space
        pointsDomain0[i] = float(coordinatesPointer[i]);
    }
    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++) {
        // 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;
        }
    }
    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 x, Vec xTilde, Vec xPhys, 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);
    for (unsigned long int i = 0; i < nPointsMyrank[0]; i++) {
        // Ux
        workPointField[i] = float(UlocalPointer[3 * i]);
        // Uy
        workPointField[i + nPointsMyrank[0]] = float(UlocalPointer[3 * i + 1]);
        // Uz
        workPointField[i + 2 * nPointsMyrank[0]] = float(UlocalPointer[3 * i + 2]);
    }
    writePointFields(timestep, 0, workPointField);
    // Restore Ulocal array
    ierr = VecRestoreArray(Ulocal, &UlocalPointer);
    CHKERRQ(ierr);

    // CELL FIELD(S)
    PetscScalar *xpp, *xp, *xt;
    VecGetArray(x, &xp);
    VecGetArray(xTilde, &xt);
    VecGetArray(xPhys, &xpp);

    for (unsigned long int i = 0; i < nCellsMyrank[0]; i++) {
        // 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]);
    }
    writeCellFields(0, workCellField);
    // Restore arrays
    VecRestoreArray(x, &xp);
    VecRestoreArray(xTilde, &xt);
    VecRestoreArray(xPhys, &xpp);

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

    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;
}

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);
}