MPIIO.h

#ifndef MPIIO_H
#define MPIIO_H

// Include necessary libraries
#include <petsc.h>
//#include <petsc-private/dmdaimpl.h>
#include <mpi.h>
#include <petsc/private/dmdaimpl.h>
#include <petscdmda.h>
#include <string>

/* -----------------------------------------------------------------------------
Authors: Niels Aage, Erik Andreassen, Boyan Lazarov, August 2013
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
-------------------------------------------------------------------------- */

class MPIIO {
  public:
    // ------------- METHODS ------------------------------------------

    MPIIO(DM da_nodes, int nPfields, std::string pnames, int nCfields, std::string cnames,std::string filename);
    ~MPIIO();

    // NOT CLEAN INTERFACE: REPLACE BY STD::PAIR OR SUCH !!!!!!
    PetscErrorCode WriteVTK(DM da_nodes, Vec x, Vec xTilde, Vec xPhys,   PetscInt itr);
    
	 PetscErrorCode WriteVTK(DM da_nodes, Vec U, Vec x,      Vec xTilde, Vec xPhys, PetscInt itr);
	 PetscErrorCode WriteVTK(DM da_nodes, Vec U, Vec xPhys, PetscInt itr);
    PetscErrorCode WriteVTK(DM da_nodes, Vec U, Vec DU,     Vec DDU,    Vec x, Vec xtilde, Vec xPhys, PetscInt itr);

  private:
    // -------------- METHODS -----------------------------------------
    int lenWorkPointField;
    void abort(std::string errorMsg, std::string position);

    unsigned long int sum(unsigned long int* startPos, unsigned long int nel);
    // --------------- MEMBERS ----------------------------------------
    int         MPI_IS;               //!< The size of an MPI unsigned long integer in bytes
    int         MPI_FS;               //!< The size of an MPI float in bytes
    int         MPI_CS;               //!< The size of an MPI char in bytes
    int         nDom;                 //!< Number of domains
    int*        nPFields;             //!< Number of point fields in each domain
    int*        nCFields;             //!< Number of cell/element fields in each domain
    MPI_Offset  offset;               //!< The offset of each thread in the file
    int         rank;                 //!< The processor rank
    int         ncpu;                 //!< The number of cpus
    int         nodesPerElement;      //!< Number of nodes per element
    bool        firstFieldOutputDone; //!< Will be set to true after first field output
    std::string filename;             //!< Output filename
    MPI_File    fh;                   //!< Filehandle

    void 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);
    // std::string filename = "/home/naage/PETSc/output2.dat");

    void writePoints(int domain, float coordinates[]);

    void writeCells(int domain, unsigned long int elements[], unsigned long int cellsOffset0[],
                    unsigned long int cellsTypes0[]);

    void writePointFields(unsigned long int timeStep, int domain, float fields[],
                          std::string newFilename = "notDefined");

    void writeCellFields(int domain, float fields[]);
    // ------------ MEMBERS  - maybe they can be private too -----------
    unsigned long int* nPoints;  //!< The number of points in each domain in each thread
    unsigned long int* nCells;   //!< The number of elements/cells in each domain in each thread
    unsigned long int* nPointsT; //!< The total number of points in each domain
    unsigned long int* nCellsT;  //!< The total number of cells in each domain

    // Converters needed for PETSc adaptation
    unsigned long int *nPointsMyrank, *nCellsMyrank;
    float *            workPointField, *workCellField;
    PetscErrorCode     DMDAGetElements_3D(DM dm, PetscInt* nel, PetscInt* nen, const PetscInt* e[]);
};
/** @example
        An illustrative example explaining how to use the class
        @code
  string userDefined = "something"
  int nDom 		 = 2;
  int nPFields[nDom]= {2, 1};
  string pFieldNames = "density, pdensity, density";
  int nCFields[nDom]= {1, 1};
  string cFieldNames = "density, density";
  int nPointsMyrank[nDom] = {10, 20}; // Arbitrary numbers
  int nCellsMyrank[nDom] = {9, 19};

  // Initialize
  MPIIO outputObject = DFEMMMPIIO(userDefined, nDom, nPFields, pFieldNames,
                                                                nCFields,
  cFieldNames, rank, nPoints, nCells);
  // Both the total number of points/cells and the corresponding
  // processor dependent numbers are determined in the constructor!
  // This does of course require MPI communication, but is only done once.


  // Output coordinates/points/vertices:
  // First domain:
  float points1[3*nPoints[0]];
  // ... put coordinates into the array points1
  outputObject.writePoints(0, points1);
  // Second domain:
  float points1[3*nPoints[1]];
  // ... put coordinates into the array points2
  outputObject.writePoints(1, points2);

  // Output elements:
  // First domain:
  int elements1[USER_SPECIFIED_SIZE];
  // ... put elements into the array elements1
  outputObject.writeCells(0, elements1);
  // Second domain
  int elements2[USER_SPECIFIED_SIZE];
  // ... put elements into the array
  outputObject.writeCells(1,  elements2);

  // The thought is that elements should be a long list with integers, where one
  integer
  // specifies the element type and the following integers specify point
  connectivity.

  // ...  Perform your computations

  // Output fields at a time point (timeStep)
  // First point fields are outputted, then cell fields

  // Output point field:
  // First domain:
  // Put all field variables into an array pFieldsD1, should of course
  // be in the same order as your point list, and the fields should come
  // after each other.
  outputObject.writePointField(timeStep, 0, pFieldsD1);
  // Second domain:
  // Put all field variables into an array pFieldsD2.
  outputObject.writePointField(1, pFieldsD2);

  // Output cell/element fields:
  // First domain:
  // Put all field variables into an array cFieldsD1, should of course
  // be in the same order as your element list, and the fields should come after
  each other. outputObject.writePointField(0, cFieldsD1);
  // Second domain:
  // Put all field variables into an array cFieldsD2.
  outputObject.writePointField(1, cFieldsD2);

  // Repeat as many times as you want.
        @endcode
*/

#endif