samplingdata.h

Go to the documentation of this file.

// Copyright (C) 2024 The m-AIA AUTHORS
//
// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
//
// SPDX-License-Identifier: LGPL-3.0-only
 
#ifndef SAMPLINGDATA_H_
#define SAMPLINGDATA_H_
 
#include <numeric>
#include "COMM/mpioverride.h"
#include "GEOM/geometry.h"
#include "GEOM/geometry2d.h"
#include "GEOM/geometry3d.h"
#include "GEOM/geometryelement.h"
#include "IO/parallelio.h"
#include "UTIL/functions.h"
#include "UTIL/tensor.h"
#include "globals.h"
#include "typetraits.h"
 
// TODO labels:PP,COMM Put all ranks which record data with this feature in a seperate MPI
// communicator and use that during I/O for optimization
 
// Predeclaration of auxiliary claas
class SamplingDataSeries;
template <MInt nDim, class SolverType>
class SamplingData {
 public:
  SamplingData(SolverType& solver) : m_solver(solver) {}
 
  void init();
  void save(MBool finalTimestep);
  void setInputOutputProperties();
  MBool enabled() const { return m_enabled; }
 
  virtual MInt loadInputFile(const MString NotUsed(inputFileName),
                             const MInt NotUsed(fileNo),
                             std::vector<MFloat>& NotUsed(coordinates)) {
    TERMM(1, "ERROR: not implemented in base class");
    return -1;
  }
  virtual MString getBaseName() const {
    TERMM(1, "ERROR: not implemented in base class");
    return "";
  }
 
 protected:
  // Methods
  void init(SamplingDataSeries& timeSeries);
  void save(SamplingDataSeries& timeSeries, MBool finalTimeStep);
 
  virtual MBool hasInputDataFile(const MInt fileNo) const {
    const MString propName = getBaseName() + "DataFileName_" + std::to_string(fileNo);
    if(Context::propertyExists(propName, solverId())) {
      const MString fileName = Context::getSolverProperty<MString>(propName, solverId(), AT_);
      return fileName != "";
    } else {
      return false;
    }
  }
 
  virtual MBool generatePoints() const { return false; }
 
  virtual MString getInputDataFile(const MInt fileNo) const {
    const MString propName = getBaseName() + "DataFileName_" + std::to_string(fileNo);
    return Context::getSolverProperty<MString>(propName, solverId(), AT_);
  }
 
  virtual MInt noSamplingPoints(const MInt NotUsed(fileNo)) {
    TERMM(1, "Not implemented in derived class!");
    return -1;
  }
 
  virtual void getSamplingPoint(const MInt NotUsed(fileNo), const MInt NotUsed(pointId),
                                MFloat* const NotUsed(coordinates)) {
    TERMM(1, "Not implemented in derived class!");
  }
 
  void indexXD(const MInt index, const MInt* const nPoints, MInt* indexXD);
 
  void saveSamplingPointCoordinates(SamplingDataSeries& timeSeries);
 
  virtual void readAdditionalProperties(const MInt NotUsed(fileNo)){};
 
  // Convenience wrappers to reduce direct solver access
  MInt solverId() const { return m_solver.solverId(); }
  MBool isMpiRoot() const { return m_solver.domainId() == 0; }
  MPI_Comm mpiComm() const { return m_solver.mpiComm(); }
 
  // Other convenience methods
  MInt noVars() const { return std::accumulate(m_noSolverSamplingVars.begin(), m_noSolverSamplingVars.end(), 0); }
  MString getFileBaseName() const { return getBaseName() + "_data_"; }
 
  SolverType& solver() { return m_solver; }
  const SolverType& solver() const { return m_solver; }
 
  // Members
  // Used to disable this entire feature
  MBool m_enabled = false;
  MBool m_interpolatePointData;
 
  std::vector<SamplingDataSeries> m_timeSeries{};
 
  std::vector<MInt> m_solverSamplingVarIds{};
  std::vector<MInt> m_noSolverSamplingVars{};
  std::vector<std::vector<MString>> m_solverSamplingVarNames{};
 
 private:
  // Members
  // Store reference to solver
  SolverType& m_solver;
};
 
 
template <MInt nDim, class SolverType>
class PointData : public SamplingData<nDim, SolverType> {
 public:
  using Base = SamplingData<nDim, SolverType>;
  // Allow the use of base class functions without the need to use this->function()
  using Base::isMpiRoot;
  using Base::m_enabled;
  using Base::m_timeSeries;
  using Base::mpiComm;
  using Base::solver;
  using Base::solverId;
 
  // Constructor that must be called; sets reference to solver
  PointData(SolverType& solver_) : SamplingData<nDim, SolverType>(solver_) {}
  virtual ~PointData(){};
  // Override virtual methods of base class
  MInt loadInputFile(const MString inputFileName, const MInt NotUsed(fileNo),
                     std::vector<MFloat>& coordinates) override;
 
  MString getBaseName() const override { return "point"; }
};
 
 
template <MInt nDim, class SolverType>
class SurfaceData : public SamplingData<nDim, SolverType> {
 public:
  using Base = SamplingData<nDim, SolverType>;
  // Allow the use of base class functions without the need to use this->function()
  using Base::isMpiRoot;
  using Base::m_enabled;
  using Base::m_timeSeries;
  using Base::mpiComm;
  using Base::solver;
  using Base::solverId;
 
  // Constructor that must be called; sets reference to solver
  SurfaceData(SolverType& solver_) : SamplingData<nDim, SolverType>(solver_) {}
  virtual ~SurfaceData();
 
  // Store geometric information for surface
  void saveGeomInfo(const MString inputFileName, const MInt fileNo, MFloat* coordinates);
 
  // Override virtual methods of base class
  MInt loadInputFile(const MString inputFileName, const MInt fileNo, std::vector<MFloat>& coordinates) override;
  MString getBaseName() const override { return "surface"; }
 
 private:
  // Geometry object with surface elements
  typename GeometryXD<nDim>::type* m_geometry = nullptr;
};
 
 
template <MInt nDim, class SolverType>
class VolumeData : public SamplingData<nDim, SolverType> {
 public:
  using Base = SamplingData<nDim, SolverType>;
  // Allow the use of base class functions without the need to use this->function()
  using Base::indexXD;
  using Base::isMpiRoot;
  using Base::m_enabled;
  using Base::m_timeSeries;
  using Base::mpiComm;
  using Base::solver;
  using Base::solverId;
 
  // Constructor that must be called; sets reference to solver
  VolumeData(SolverType& solver_) : SamplingData<nDim, SolverType>(solver_) {}
  virtual ~VolumeData(){};
  // Override virtual methods of base class
  MString getBaseName() const override { return "volume"; }
 
 protected:
  MBool generatePoints() const override { return true; }
 
  MBool hasInputDataFile(const MInt fileNo) const override {
    const MString propName = getBaseName() + "DataShape_" + std::to_string(fileNo);
    if(Context::propertyExists(propName, solverId())) {
      const MString shape = Context::getSolverProperty<MString>(propName, solverId(), AT_);
      return shape != "";
    } else {
      return false;
    }
  }
 
  MString getInputDataFile(const MInt NotUsed(fileNo)) const {
    return ""; // No input file, points are generated
  }
 
  void readAdditionalProperties(const MInt fileNo) override;
  // Sampling point generation
  MInt noSamplingPoints(const MInt fileNo) override;
  void getSamplingPoint(const MInt fileNo, const MInt pointId, MFloat* const coordinates) override;
 
 private:
  std::vector<MString> m_volumeShape{};
  std::vector<std::vector<MFloat>> m_volumeParameterFloat{};
  std::vector<std::vector<MInt>> m_volumeParameterInt{};
};
 
 
class SamplingDataSeries {
 public:
  // C'tor
  SamplingDataSeries(const MString inputFileName, const MBool generatePoints, const MInt sampleInterval,
                     const MInt writeInterval, const MInt startTimeStep, const MInt endTimeStep, const MInt noFile)
    : m_inputFileName(inputFileName),
      m_generatePoints(generatePoints),
      m_sampleInterval(sampleInterval),
      m_writeInterval(writeInterval),
      m_startTimeStep(startTimeStep),
      m_endTimeStep(endTimeStep),
      m_fileNo(noFile) {}
 
  // Member
  // File name of file containing points
  const MString m_inputFileName = "";
  // Enable generation of points instead of reading from file
  const MBool m_generatePoints = false;
  // Interval in which samplingData should be saved
  const MInt m_sampleInterval = -1;
  // Interval in which samplingData should be written
  const MInt m_writeInterval = -1;
  // Time step after which recording of conservative variables starts
  const MInt m_startTimeStep = -1;
  // Time step after which recording of conservative variables ends
  const MInt m_endTimeStep = -1;
 
  // The following three vectors contain only the data relevant to the rank
  // Vector containing the point coordinates at which states should be written
  // out
  // Contains nDim tuples which represent one coordinate
  std::vector<MFloat> m_coordinates{};
  // Vector containing the element/cell ids of the points
  std::vector<MInt> m_elementIds{};
  // Sort index: each entry points to the index in the input file
  std::vector<MInt> m_sortIndex{};
  // Number of total points at which states should be written out
  MInt m_noPoints = 0;
  // Numver of points relevant to this rank
  MInt m_noLocalPoints = 0;
  // Buffer which holds timesteps intbetween writes of point data
  std::vector<MInt> m_timeStepBuffer{};
  // Buffer which holds times intbetween writes of point data
  std::vector<MFloat> m_timeBuffer{};
  // Buffer which holds sampling data states of variables inbetween writes
  // Should be accessed as 3D array with the following dimensions:
  // 1st dim: point index
  // 2nd dim: sample index
  // 3rd dim: var index
  std::vector<MFloat> m_stateBuffer{};
  // How many times point data was saved to the buffer
  MInt m_sampleIndex = 0;
  // Maximum times point data will be buffered
  MInt m_maxNoSample = 0;
  // MPI Communicator containing point data ranks
  MPI_Comm m_mpiComm = MPI_COMM_NULL;
  // Identifies the filenumber of current object
  const MInt m_fileNo = -1;
  // Identifies new root for point data ranks
  MBool m_isMpiRoot = false;
};
 
 
/* \brief Sets I/O properties concerning this feature
 *
 * \author Fabian Klemp <f.klemp@aia.rwth-aachen.de>
 * \date 2016-09-06
 *
 */
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::setInputOutputProperties() {
  TRACE();
 
  // Get base name of derived class, all properties are defined as "[baseName][PropertyName]"
  const MString type = getBaseName();
 
  if(Context::propertyExists(type + "DataEnabled", solverId())) {
    m_enabled = Context::getSolverProperty<MBool>(type + "DataEnabled", solverId(), AT_);
    if(!m_enabled) {
      return;
    }
  }
 
  // Check for all input files
  if(hasInputDataFile(0)) {
    m_enabled = true;
 
    // Check for new default variables
    MInt sampleIntervalGlobal = 1;
    sampleIntervalGlobal =
        Context::getSolverProperty<MInt>(type + "DataSampleInterval", solverId(), AT_, &sampleIntervalGlobal);
 
    MInt writeIntervalGlobal = -1;
    const MInt restartInterval = m_solver.restartInterval();
    if(restartInterval > 0) {
      writeIntervalGlobal =
          Context::getSolverProperty<MInt>(type + "DataWriteInterval", solverId(), AT_, &restartInterval);
    } else {
      writeIntervalGlobal =
          Context::getSolverProperty<MInt>(type + "DataWriteInterval", solverId(), AT_, &writeIntervalGlobal);
    }
 
    MInt startTimeStepGlobal = 0;
    startTimeStepGlobal =
        Context::getSolverProperty<MInt>(type + "DataStartTimeStep", solverId(), AT_, &startTimeStepGlobal);
 
    MInt endTimeStepGlobal = std::numeric_limits<MInt>::max();
    endTimeStepGlobal = Context::getSolverProperty<MInt>(type + "DataEndTimeStep", solverId(), AT_, &endTimeStepGlobal);
 
    m_solverSamplingVarIds.clear();
    m_noSolverSamplingVars.clear();
    m_solverSamplingVarNames.clear();
    // Get sampling variables from solver
    solver().getSolverSamplingProperties(m_solverSamplingVarIds, m_noSolverSamplingVars, m_solverSamplingVarNames);
 
    // Allocate additional storage in the solver for the sampling variables
    solver().initSolverSamplingVariables(m_solverSamplingVarIds, m_noSolverSamplingVars);
 
    // TODO labels:PP add output about sampling variables
 
    // Read all Input files, properties and create point data objects
    MInt noFiles = 0;
    while(hasInputDataFile(noFiles)) {
      const MString samplingDataFileName = getInputDataFile(noFiles);
 
      // Check if points should be generated if no input file name is specified
      MBool genPoints = false;
      if(samplingDataFileName == "") {
        genPoints = generatePoints(); // Default behaviour
        genPoints = Context::getSolverProperty<MBool>(type + "DataGeneratePoints_" + std::to_string(noFiles),
                                                      solverId(), AT_, &genPoints);
      }
 
      // Check if other property is set by file number
      const MInt sampleInterval = Context::getSolverProperty<MInt>(
          type + "DataSampleInterval_" + std::to_string(noFiles), solverId(), AT_, &sampleIntervalGlobal);
 
      const MInt writeInterval = Context::getSolverProperty<MInt>(type + "DataWriteInterval_" + std::to_string(noFiles),
                                                                  solverId(), AT_, &writeIntervalGlobal);
 
      if(sampleInterval < 1) {
        TERMM(1, type + "DataSampleInterval must be a positive integer.");
      }
 
      if(writeInterval < 1) {
        TERMM(1, type
                     + "DataWriteInterval musst be a positive integer. Check if "
                       "all write intervals are set.");
      }
 
      if(writeInterval % sampleInterval != 0) {
        TERMM(1, "samplingDataSampleInterval must be a divisor of "
                 "samplingDataWriteInterval (default value: restartInterval) for "
                 "technical reasons. Also it makes more sense this way.");
      }
      if(m_solver.restartInterval() > 0 && m_solver.restartInterval() % writeInterval != 0) {
        TERMM(1, type
                     + "DataWriteInterval must be a divisor of restartInterval because otherwise "
                       "restarting a simulation without data loss is not possible.");
      }
 
      // For property description see samplingDataStartTimeStep
      const MInt startTimeStep = Context::getSolverProperty<MInt>(type + "DataStartTimeStep_" + std::to_string(noFiles),
                                                                  solverId(), AT_, &startTimeStepGlobal);
      if(startTimeStep < 0) {
        TERMM(1,
              "Point data recording start time step is negative. It must have "
              "a positive value.");
      }
 
      const MInt endTimeStep = Context::getSolverProperty<MInt>(type + "DataEndTimeStep_" + std::to_string(noFiles),
                                                                solverId(), AT_, &endTimeStepGlobal);
 
      if(endTimeStep < 0) {
        TERMM(1, "Point data recording end time step is negative. It must have "
                 "a positive value.");
      }
      if(endTimeStep < startTimeStep) {
        TERMM(1,
              "End time step of point data recording is smaller than start time"
              " step of point data recording. Check your property file.");
      }
 
      m_interpolatePointData = true;
      if(Context::propertyExists("interpolatePointData")) {
        m_interpolatePointData =
            Context::getSolverProperty<MBool>("interpolatePointData", solverId(), AT_, &m_interpolatePointData);
      }
 
      // Read additional properties if required by derived class
      readAdditionalProperties(noFiles);
 
      // Create auxiliary object
      SamplingDataSeries samplingData(samplingDataFileName, genPoints, sampleInterval, writeInterval, startTimeStep,
                                      endTimeStep, noFiles);
      // Store all auxiliary objects to edit them in a loop
      m_timeSeries.push_back(samplingData);
      noFiles++;
    }
  } else {
    m_enabled = false;
  }
}
 
 
/* \brief Initializes all auxiliary sampling data objects
 *
 * \author Marcus Wiens <m.wiens@aia.rwth-aachen.de>
 * \date 2016-11-02
 */
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::init() {
  if(!m_enabled || !solver().isActive()) {
    return;
  }
 
  for(MUint i = 0; i < m_timeSeries.size(); i++) {
    init(m_timeSeries[i]);
  }
}
 
 
/* \brief Initializes all sampling data related members
 * Load all points relevant to this rank at which states should be
 * written out.
 *
 * \author Fabian Klemp <f.klemp@aia.rwth-aachen.de>
 * \date 2016-06-19
 */
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::init(SamplingDataSeries& timeSeries) {
  using namespace maia;
  using namespace std;
 
 
  // save points and corresponding data in this container for sorting
  // tuple contains following data in that order:
  // coordinates of point, element id of point, line index in input file
  std::vector<std::tuple<std::array<MFloat, nDim>, MInt, MInt>> localTmpPoints;
 
  m_log << "Solver #" << solverId() << ": ";
  // Check for an input file
  if(timeSeries.m_inputFileName != "") {
    m_log << "Loading points for " << getBaseName() << " data class...";
 
    std::vector<MFloat> tmpCoordinates{};
 
    // read in file on root using the I/O method of the derived class
    if(isMpiRoot()) {
      timeSeries.m_noPoints = loadInputFile(timeSeries.m_inputFileName, timeSeries.m_fileNo, tmpCoordinates);
    }
 
    MPI_Bcast(&timeSeries.m_noPoints, 1, type_traits<MInt>::mpiType(), 0, mpiComm(), AT_, "timeSeries.m_noPoints");
 
    if(timeSeries.m_noPoints <= 0) {
      TERMM(1, "ERROR: no sampling points found for input file '" + timeSeries.m_inputFileName + "'");
      m_log << "ERROR: no sampling points found for input file '" << timeSeries.m_inputFileName << "'" << std::endl;
      return;
    }
 
    ScratchSpace<MFloat> coordinatesScratch(timeSeries.m_noPoints * nDim, FUN_, "coordinatesScratch");
    if(isMpiRoot()) {
      copy(tmpCoordinates.begin(), tmpCoordinates.end(), coordinatesScratch.begin());
    }
    MPI_Bcast(&coordinatesScratch[0], timeSeries.m_noPoints * nDim, type_traits<MFloat>::mpiType(), 0, mpiComm(), AT_,
              "coordinatesScratch[0]");
 
    localTmpPoints.resize(timeSeries.m_noPoints);
 
    MFloatTensor tmpCoordinatesTensor(&coordinatesScratch[0], timeSeries.m_noPoints, nDim);
    std::array<MFloat, nDim> curCoordinates{};
    MInt curElementId = -1;
    for(MInt i = 0; i < timeSeries.m_noPoints; i++) {
      for(MInt j = 0; j < nDim; j++) {
        curCoordinates[j] = tmpCoordinatesTensor(i, j);
      }
      // check if the point is relevant to this rank, each point should be globally unique, i.e.,
      // appear only once
      curElementId = solver().getIdAtPoint(&curCoordinates[0], true);
 
      if(curElementId != -1) {
        localTmpPoints[timeSeries.m_noLocalPoints] = make_tuple(curCoordinates, curElementId, i);
        timeSeries.m_noLocalPoints++;
 
        // Precompute interpolation information if possible/required
        solver().initInterpolationForCell(curElementId);
      }
    }
 
    // sanity check if each point ended up on exactly one rank:
    ScratchSpace<MInt> loadedPointsVector(timeSeries.m_noPoints, FUN_, "loadedPointsVector");
    loadedPointsVector.fill(0);
    for(MInt i = 0; i < timeSeries.m_noLocalPoints; i++) {
      loadedPointsVector[get<2>(localTmpPoints[i])]++;
    }
 
    if(isMpiRoot()) {
      MPI_Reduce(MPI_IN_PLACE, &loadedPointsVector[0], loadedPointsVector.size(), type_traits<MInt>::mpiType(), MPI_SUM,
                 0, mpiComm(), AT_, "MPI_IN_PLACE", "loadedPointsVector[0]");
      for(MUint i = 0; i < loadedPointsVector.size(); i++) {
        if(loadedPointsVector[i] != 1) {
          std::ostringstream err;
          err << "Point at line " << i + 1 << " of file " // start line count at one
              << timeSeries.m_inputFileName << " was loaded " << loadedPointsVector[i] << " times. Aborting."
              << std::endl;
          TERMM(1, err.str());
        }
      }
    } else {
      MPI_Reduce(&loadedPointsVector[0], nullptr, loadedPointsVector.size(), type_traits<MInt>::mpiType(), MPI_SUM, 0,
                 mpiComm(), AT_, "loadedPointsVector[0]", "nullptr");
    }
  } else if(timeSeries.m_generatePoints) {
    m_log << "Generating points for " << getBaseName() << " data class...";
    // Generate points based on given properties
    // Get number of points to generate
    const MInt noPoints = noSamplingPoints(timeSeries.m_fileNo);
    timeSeries.m_noPoints = noPoints;
 
    std::array<MFloat, nDim> coordinates{};
    MLong checkSum = 0;
 
    for(MInt i = 0; i < noPoints; i++) {
      // Get the coordinates for this sampling point id
      getSamplingPoint(timeSeries.m_fileNo, i, &coordinates[0]);
 
      // Check if point is relevant for local rank
      const MInt curElementId = solver().getIdAtPoint(&coordinates[0], true);
      if(curElementId != -1) {
        localTmpPoints.push_back(make_tuple(coordinates, curElementId, i));
        timeSeries.m_noLocalPoints++;
        checkSum += (i + 1);
 
        // Precompute interpolation information if possible/required
        solver().initInterpolationForCell(curElementId);
      }
    }
 
    MInt globalNoPoints = timeSeries.m_noLocalPoints;
    MPI_Allreduce(MPI_IN_PLACE, &globalNoPoints, 1, type_traits<MInt>::mpiType(), MPI_SUM, mpiComm(), AT_,
                  "MPI_IN_PLACE", "globalNoPoints");
    if(globalNoPoints != noPoints) {
      TERMM(1, "Global number of points does not match " + std::to_string(globalNoPoints)
                   + " != " + std::to_string(noPoints));
    }
 
    // Use Gauss sum equation to check if all points ended up on exactly one rank
    // Note: this is a heuristic check but it should be very unlikely to get a false positive; this
    // will work up to 2 billion sampling points
    MPI_Allreduce(MPI_IN_PLACE, &checkSum, 1, type_traits<MLong>::mpiType(), MPI_SUM, mpiComm(), AT_, "MPI_IN_PLACE",
                  "checkSum");
    const MLong noPointsTmp = noPoints;
    if(2 * checkSum != (noPointsTmp * noPointsTmp + noPointsTmp)) {
      TERMM(1, "Incorrect check sum for " + std::to_string(noPoints) + " points: " + std::to_string(checkSum));
    }
  } else {
    m_log << "Warning: no sampling data input file specified and generation of points disabled." << std::endl;
    return;
  }
 
  // Create a new MPI Communicator
  ScratchSpace<MInt> pointDataScratch(globalNoDomains(), FUN_, "pointDataScratch");
  std::fill_n(&pointDataScratch[0], globalNoDomains(), -1);
 
  // Get Rank if points are relevant
  MInt rank = -1;
  if(timeSeries.m_noLocalPoints > 0) {
    MPI_Comm_rank(mpiComm(), &rank);
  }
 
  // Combine all ranks to get relevant ranks
  MPI_Allgather(&rank, 1, type_traits<MInt>::mpiType(), &pointDataScratch[0], 1, type_traits<MInt>::mpiType(),
                mpiComm(), AT_, "rank", "pointDataScratch[0]");
 
  // Check for relevant ranks and save them to create the new communicator
  const MInt noPdDomains =
      count_if(pointDataScratch.begin(), pointDataScratch.end(), [](const MInt a) { return a != -1; });
  MIntScratchSpace pdDomains(noPdDomains, AT_, "pdDomains");
  MInt position = 0;
  for(MUint i = 0; i < pointDataScratch.size(); i++) {
    if(pointDataScratch[i] != -1) {
      pdDomains[position] = pointDataScratch[i];
      position++;
    }
  }
 
  // Create new point data mpi group
  MPI_Group globalGroup, localGroup;
  MPI_Comm_group(mpiComm(), &globalGroup, AT_, "globalGroup");
  MPI_Group_incl(globalGroup, noPdDomains, &pdDomains[0], &localGroup, AT_);
 
  // Create new communicator and clean up
  MPI_Comm_create(mpiComm(), localGroup, &timeSeries.m_mpiComm, AT_, "timeSeries.m_mpiComm");
  MPI_Group_free(&globalGroup, AT_);
  MPI_Group_free(&localGroup, AT_);
 
  // Set root for new communicator
  if(timeSeries.m_noLocalPoints > 0) {
    MInt subRank = -1;
    MPI_Comm_rank(timeSeries.m_mpiComm, &subRank);
    if(subRank == 0) {
      timeSeries.m_isMpiRoot = true;
    }
  }
 
 
  // sort points by global id
  stable_sort(localTmpPoints.begin(),
              localTmpPoints.begin() + timeSeries.m_noLocalPoints,
              // Read all input files and save all points into one array
              [this](const std::tuple<std::array<MFloat, nDim>, MInt, MInt>& a,
                     const std::tuple<std::array<MFloat, nDim>, MInt, MInt>& b) {
                return solver().grid().tree().globalId(solver().getCellIdByIndex(get<1>(a)))
                       < solver().grid().tree().globalId(solver().getCellIdByIndex(get<1>(b)));
              });
 
  // split localTmpPoints for it to be easier to handle later on
  if(timeSeries.m_noLocalPoints > 0) {
    timeSeries.m_coordinates.resize(timeSeries.m_noLocalPoints * nDim);
    MFloatTensor coordinatesTensor(&timeSeries.m_coordinates[0], timeSeries.m_noLocalPoints, nDim);
    timeSeries.m_elementIds.resize(timeSeries.m_noLocalPoints);
    timeSeries.m_sortIndex.resize(timeSeries.m_noLocalPoints);
    for(MInt i = 0; i < timeSeries.m_noLocalPoints; i++) {
      for(MInt j = 0; j < nDim; j++) {
        coordinatesTensor(i, j) = get<0>(localTmpPoints[i])[j];
      }
      timeSeries.m_elementIds[i] = get<1>(localTmpPoints[i]);
      timeSeries.m_sortIndex[i] = get<2>(localTmpPoints[i]);
    }
  }
 
  // Save sampling point coordinates if they were not loaded from a file
  if(timeSeries.m_inputFileName == "" && timeSeries.m_noLocalPoints > 0) {
    saveSamplingPointCoordinates(timeSeries);
  }
 
  // Allocate memory to save/buffer point data inbetween writes
  timeSeries.m_maxNoSample = std::max(1, timeSeries.m_writeInterval / timeSeries.m_sampleInterval);
  if(timeSeries.m_noLocalPoints > 0) {
    // Note: +1 for initial data file containing also timestep 0
    timeSeries.m_stateBuffer.resize(timeSeries.m_noLocalPoints * (timeSeries.m_maxNoSample + 1) * noVars());
  }
 
  // Allocate storage for the time step count and the simulation time
  if(timeSeries.m_isMpiRoot) {
    timeSeries.m_timeStepBuffer.resize(timeSeries.m_maxNoSample + 1);
    timeSeries.m_timeBuffer.resize(timeSeries.m_maxNoSample + 1);
  } else {
    timeSeries.m_timeStepBuffer.resize(1);
    timeSeries.m_timeBuffer.resize(1);
  }
 
  m_log << " done" << std::endl << "#" << timeSeries.m_fileNo;
  if(timeSeries.m_inputFileName != "") {
    m_log << " Loaded " << timeSeries.m_noPoints << " points from " << timeSeries.m_inputFileName << " for "
          << getBaseName() << " sampling feature successfully." << std::endl;
  } else if(timeSeries.m_generatePoints) {
    m_log << " Generated " << timeSeries.m_noPoints << " points for " << getBaseName()
          << " sampling feature successfully." << std::endl;
  }
  m_log << "sampleInterval = " << timeSeries.m_sampleInterval << "; startTimeStep = " << timeSeries.m_startTimeStep
        << "; endTimeStep = " << timeSeries.m_endTimeStep << "; writeInterval = " << timeSeries.m_writeInterval
        << std::endl;
}
 
 
/* \brief Buffers/Save all auxiliary sampling data objects
 *
 * \author Marcus Wiens <m.wiens@aia.rwth-aachen.de>
 * \date 2016-11-02
 */
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::save(MBool finalTimeStep) {
  if(!m_enabled) {
    return;
  }
 
  for(MUint i = 0; i < m_timeSeries.size(); i++) {
    save(m_timeSeries[i], finalTimeStep);
  }
}
 
 
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::save(SamplingDataSeries& timeSeries, MBool finalTimeStep) {
  using namespace std;
 
  // Check if time step inverval is reached
  if(solver().getCurrentTimeStep() < timeSeries.m_startTimeStep
     || solver().getCurrentTimeStep() > timeSeries.m_endTimeStep) {
    return;
  }
 
  // Check for correct timestep interval
  if(!(solver().getCurrentTimeStep() % timeSeries.m_sampleInterval == 0 || finalTimeStep)) {
    return;
  }
 
  // Compute sampling variables and exchange data for interpolation (if not done already in this time step)
  solver().calcSamplingVariables(m_solverSamplingVarIds, true);
 
  // Check if rank belongs to an sampling data MPI communicator
  // Note: sampling variables are computed by all ranks as the interpolation might require halo cell
  // values
  if(timeSeries.m_noLocalPoints < 1 || timeSeries.m_noPoints == 0) {
    return;
  }
 
  // buffer time and timestep
  if(timeSeries.m_isMpiRoot) {
    timeSeries.m_timeStepBuffer[timeSeries.m_sampleIndex] = solver().getCurrentTimeStep();
    timeSeries.m_timeBuffer[timeSeries.m_sampleIndex] = solver().time();
  }
 
  // buffer point data states
  MFloatTensor coordinatesTensorBuf(&timeSeries.m_coordinates[0], timeSeries.m_noLocalPoints, nDim);
  MFloatTensor stateTensorBuf(&timeSeries.m_stateBuffer[0], timeSeries.m_noLocalPoints, timeSeries.m_maxNoSample + 1,
                              noVars());
 
  MBool interpolate = m_interpolatePointData;
  for(MInt i = 0; i < timeSeries.m_noLocalPoints; i++) {
    MInt offset = 0;
    for(MUint var = 0; var < m_solverSamplingVarIds.size(); var++) {
      solver().calcSamplingVarAtPoint(&coordinatesTensorBuf(i, 0), timeSeries.m_elementIds[i],
                                      m_solverSamplingVarIds[var], &stateTensorBuf(i, timeSeries.m_sampleIndex, offset),
                                      interpolate);
      offset += m_noSolverSamplingVars[var];
    }
  }
 
  // increment sample index
  timeSeries.m_sampleIndex++;
 
  MBool isFirstDataFile = false;
  // Note: first data file for simulation starting at timestep 0 contains one value more than
  // all other data files; at a restart the restart time step is not considered here even if it is
  // equal to m_startTimeStep
  // if (solver().getCurrentTimeStep() >= timeSeries.m_startTimeStep
  //   && solver().getCurrentTimeStep() <= timeSeries.m_startTimeStep + timeSeries.m_writeInterval) {
  if(solver().getCurrentTimeStep() >= 0 && solver().getCurrentTimeStep() <= timeSeries.m_writeInterval) {
    isFirstDataFile = true;
  }
 
  const MBool isRestartStep =
      (solver().restartInterval() > 0 && solver().getCurrentTimeStep() % solver().restartInterval() == 0);
  const MBool isFirstSample = (isFirstDataFile && timeSeries.m_sampleIndex == 1);
  const MBool writeAtRestartStep = (isRestartStep && !isFirstSample);
  const MBool writeFirstDataFile = (timeSeries.m_sampleIndex > timeSeries.m_maxNoSample && isFirstDataFile);
  const MBool writeDataFile = ((timeSeries.m_sampleIndex + 1 > timeSeries.m_maxNoSample) && !isFirstDataFile);
 
  // write file if it's time to do so otherwise return
  if(!(writeDataFile || finalTimeStep || writeAtRestartStep || writeFirstDataFile)) {
    return;
  }
 
  using namespace maia::parallel_io;
 
  // calculate how many samples have been recorded in total
  MInt sampleCount = (solver().getCurrentTimeStep() / timeSeries.m_sampleInterval) + 1; // also saving at timeStep 0
 
  // Correct number of samples if finalTimeStep is not in the writeInterval
  if(finalTimeStep && !(solver().getCurrentTimeStep() % timeSeries.m_writeInterval == 0)) {
    sampleCount++;
  }
 
  MInt fileNo = 0;
  if(!finalTimeStep || (solver().getCurrentTimeStep() % timeSeries.m_writeInterval == 0)) {
    fileNo = (sampleCount - 1) / (timeSeries.m_writeInterval / timeSeries.m_sampleInterval);
  } else {
    fileNo = ceil((static_cast<MFloat>(sampleCount - 1)
                   / static_cast<MFloat>(timeSeries.m_writeInterval / timeSeries.m_sampleInterval)));
  }
 
  // Output file name
  std::ostringstream fileName;
  MString baseFileName = getFileBaseName();
  fileName << solver().outputDir() << baseFileName << std::to_string(timeSeries.m_fileNo) << "_" << setw(8)
           << setfill('0') << fileNo << ParallelIo::fileExt();
 
  ParallelIo file(fileName.str(), PIO_REPLACE, timeSeries.m_mpiComm);
 
  file.defineArray(PIO_INT, "timeStep", timeSeries.m_sampleIndex);
  file.defineArray(PIO_FLOAT, "time", timeSeries.m_sampleIndex);
  file.defineArray(PIO_INT, "sortIndex", timeSeries.m_noPoints);
 
  file.setAttribute("si[i] = index in input", "description", "sortIndex");
 
  ParallelIo::size_type dimSizes[] = {timeSeries.m_noPoints, timeSeries.m_sampleIndex, noVars()};
  file.defineArray(PIO_FLOAT, "pointStates", 3, &dimSizes[0]);
 
  std::ostringstream descPs;
  // Please be aware that these attributes get parsed in the corresponding post
  // processing script. So if you change them please check if it breaks the
  // script
  file.setAttribute(descPs.str(), "description", "pointStates");
  file.setAttribute("point index", "dim_0", "pointStates");
  file.setAttribute("sample index", "dim_1", "pointStates");
  file.setAttribute("var index", "dim_2", "pointStates");
  file.setAttribute(timeSeries.m_inputFileName, "inputFileName", "pointStates");
 
  MInt varOffset = 0;
  for(MUint var = 0; var < m_solverSamplingVarIds.size(); var++) {
    for(MInt i = 0; i < m_noSolverSamplingVars[var]; i++) {
      file.setAttribute(m_solverSamplingVarNames[var][i], "var_" + std::to_string(varOffset + i), "pointStates");
    }
    varOffset += m_noSolverSamplingVars[var];
  }
 
  // ensure data is consecutive in memory
  if(timeSeries.m_sampleIndex < timeSeries.m_maxNoSample + 1) {
    MFloatTensor stateTensor(&timeSeries.m_stateBuffer[0], timeSeries.m_noLocalPoints, timeSeries.m_maxNoSample + 1,
                             noVars());
    MFloatMatrix finalStateTensor(&timeSeries.m_stateBuffer[0], timeSeries.m_noLocalPoints, timeSeries.m_sampleIndex,
                                  noVars());
    for(MInt i = 0; i < timeSeries.m_noLocalPoints; i++) {
      for(MInt j = 0; j < timeSeries.m_sampleIndex; j++) {
        for(MInt k = 0; k < noVars(); k++) {
          finalStateTensor(i, j, k) = stateTensor(i, j, k);
        }
      }
    }
  }
 
  // root writes timestep and time data
  if(timeSeries.m_isMpiRoot) {
    file.setOffset(timeSeries.m_sampleIndex, 0);
  } else {
    file.setOffset(0, 0);
  }
  file.writeArray(&timeSeries.m_timeStepBuffer[0], "timeStep");
  file.writeArray(&timeSeries.m_timeBuffer[0], "time");
 
  // calculate offset for reverse index and pointStates
  ParallelIo::size_type offset, totalCount;
  file.calcOffset(timeSeries.m_noLocalPoints, &offset, &totalCount, timeSeries.m_mpiComm);
 
  // write sort index
  file.setOffset(timeSeries.m_noLocalPoints, offset);
  file.writeArray(&timeSeries.m_sortIndex[0], "sortIndex");
 
  // write state data to file
  file.setOffset(timeSeries.m_noLocalPoints, offset, 3);
  file.writeArray(&timeSeries.m_stateBuffer[0], "pointStates");
 
  // reset sample index
  timeSeries.m_sampleIndex = 0;
}
 
 
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::indexXD(const MInt index, const MInt* const nPoints, MInt* indexXD) {
  const MInt noPoints2D = nPoints[0] * nPoints[1];
  const MInt maxNoPoints = (nDim == 2) ? noPoints2D : noPoints2D * nPoints[2];
  if(index < 0 || index >= maxNoPoints) {
    TERMM(1, "invalid index: " + std::to_string(index));
  }
 
  // Last index (varies fastest)
  indexXD[nDim - 1] = index % nPoints[nDim - 1];
  // First/second index in 2D/3D
  indexXD[nDim - 2] = ((index - indexXD[nDim - 1]) / nPoints[nDim - 1]) % nPoints[nDim - 2];
  IF_CONSTEXPR(nDim == 3) {
    // First index in 3D
    indexXD[0] = (index - indexXD[1] * nPoints[1] - indexXD[2]) / (nPoints[2] * nPoints[1]);
  }
}
 
 
template <MInt nDim, class SolverType>
void SamplingData<nDim, SolverType>::saveSamplingPointCoordinates(SamplingDataSeries& timeSeries) {
  TRACE();
  using namespace maia::parallel_io;
 
  std::ostringstream fileName;
  fileName << solver().outputDir() << getFileBaseName() << "coordinates_" << timeSeries.m_fileNo
           << ParallelIo::fileExt();
  ParallelIo file(fileName.str(), PIO_REPLACE, timeSeries.m_mpiComm);
 
  // Add spatial dimensions
  file.setAttribute(nDim, "nDim");
 
  file.defineArray(PIO_INT, "sortIndex", timeSeries.m_noPoints);
  file.setAttribute("si[i] = index in input", "description", "sortIndex");
 
  ParallelIo::size_type dimSizes[] = {timeSeries.m_noPoints, nDim};
  file.defineArray(PIO_FLOAT, "coordinates", 2, &dimSizes[0]);
 
  ParallelIo::size_type offset, totalCount;
  file.calcOffset(timeSeries.m_noLocalPoints, &offset, &totalCount, timeSeries.m_mpiComm);
 
  // Write sort index
  file.setOffset(timeSeries.m_noLocalPoints, offset);
  file.writeArray(&timeSeries.m_sortIndex[0], "sortIndex");
 
  // Write coordinates
  file.setOffset(timeSeries.m_noLocalPoints, offset, 2);
  file.writeArray(&timeSeries.m_coordinates[0], "coordinates");
}
 
 
template <MInt nDim, class SolverType>
MInt PointData<nDim, SolverType>::loadInputFile(const MString inputFileName,
                                                const MInt NotUsed(fileNo),
                                                std::vector<MFloat>& coordinates) {
  TRACE();
  if(!isMpiRoot()) {
    TERMM(1, "PointData::loadInputFile() should only be called by mpi root process!");
  }
 
  return loadPointCoordinatesFromFile(inputFileName, nDim, coordinates);
}
 
 
// Destructor
template <MInt nDim, class SolverType>
SurfaceData<nDim, SolverType>::~SurfaceData() {
  // Free geometry object if existing
  if(m_geometry != nullptr) {
    // TODO labels:PP,IO fix deallocation
    // delete m_geometry;
  }
}
 
 
template <MInt nDim, class SolverType>
MInt SurfaceData<nDim, SolverType>::loadInputFile(const MString inputFileName,
                                                  const MInt fileNo,
                                                  std::vector<MFloat>& coordinates) {
  TRACE();
  if(!isMpiRoot()) {
    TERMM(1, "SurfaceData::loadInputFile() should only be called by mpi root process!");
  }
 
  // Create geometry object
  auto geometry = new typename GeometryXD<nDim>::type(0, inputFileName, MPI_COMM_SELF);
  // Store pointer in class member variable
  m_geometry = geometry;
 
  // Determine number of geometry/surface elements
  MInt noElements = geometry->GetNoElements();
  coordinates.resize(nDim * noElements);
 
  std::array<MFloat, nDim * nDim> vertex;
 
  // Compute element centroids and store in coordinates
  for(MInt i = 0; i < noElements; i++) {
    geometry->elements[i].getVertices(&vertex[0]);
    geometry->elements[i].calcCentroid(&vertex[0], &coordinates[i * nDim]);
  }
 
  // Write specific data for all elements into an output file
  saveGeomInfo(inputFileName, fileNo, &coordinates[0]);
 
  return noElements;
}
 
 
template <MInt nDim, class SysEqn>
void SurfaceData<nDim, SysEqn>::saveGeomInfo(const MString inputFileName, const MInt fileNo, MFloat* coordinates) {
  TRACE();
  using namespace maia::parallel_io;
 
  const MInt noElements = m_geometry->GetNoElements();
  // Define output file
  std::ostringstream fileName;
  fileName << solver().outputDir() << "surface_geometryInfo_" << fileNo << ParallelIo::fileExt();
 
  ParallelIo file(fileName.str(), PIO_REPLACE, MPI_COMM_SELF);
  file.setAttribute(inputFileName, "inputFileName");
 
  // Get vertices of an element and calculate length/surface of element
  std::vector<MFloat> vertices(noElements * nDim * nDim);
  std::vector<MFloat> geomMeasures(noElements);
  std::vector<MFloat> normalVecs(noElements * nDim);
 
  MFloat geomCenter[3]{};
 
  // Check if geometry center is defined by user
  if(Context::propertyExists("surfaceDataGeometryCenter_" + std::to_string(fileNo), solverId())) {
    for(MInt i = 0; i < nDim; i++) {
      geomCenter[i] =
          Context::getSolverProperty<MFloat>("surfaceDataGeometryCenter_" + std::to_string(fileNo), solverId(), AT_, i);
    }
  } else {
    // Determine geometry center as average of all element centroids
    for(MInt i = 0; i < noElements; i++) {
      for(MInt dim = 0; dim < nDim; dim++) {
        geomCenter[dim] += coordinates[i * nDim + dim];
      }
    }
    for(MInt dim = 0; dim < nDim; dim++) {
      geomCenter[dim] /= static_cast<MFloat>(noElements);
    }
  }
 
  m_log << "Surface geometry info: geomCenter =";
  for(MInt dim = 0; dim < nDim; dim++) {
    m_log << " " << geomCenter[dim];
  }
  m_log << std::endl;
 
  MInt counter = 0;
  for(MInt i = 0; i < noElements; i++) {
    for(MInt j = 0; j < nDim; j++) {
      for(MInt k = 0; k < nDim; k++) {
        vertices[counter] = m_geometry->elements[i].m_vertices[j][k];
        counter++;
      }
    }
 
    auto& element = m_geometry->elements[i];
    // Calc surface area/length of every element
    IF_CONSTEXPR(nDim == 3) {
      // Define two vectors for triangular surface in 3D
      const MFloat v12_x = element.m_vertices[0][0] - element.m_vertices[1][0];
      const MFloat v12_y = element.m_vertices[0][1] - element.m_vertices[1][1];
      const MFloat v12_z = element.m_vertices[0][2] - element.m_vertices[1][2];
      const MFloat v13_x = element.m_vertices[0][0] - element.m_vertices[2][0];
      const MFloat v13_y = element.m_vertices[0][1] - element.m_vertices[2][1];
      const MFloat v13_z = element.m_vertices[0][2] - element.m_vertices[2][2];
 
      // Cross product
      const MFloat v_x = v12_y * v13_z - v12_z * v13_y;
      const MFloat v_y = v12_z * v13_x - v12_x * v13_z;
      const MFloat v_z = v12_x * v13_y - v12_y * v13_x;
      // Surface area is 1/2 the length of the cross product
      geomMeasures[i] = 0.5 * sqrt(pow(v_x, 2) + pow(v_y, 2) + pow(v_z, 2));
    }
    else {
      const MFloat p1_x = element.m_vertices[0][0];
      const MFloat p2_x = element.m_vertices[1][0];
      const MFloat p1_y = element.m_vertices[0][1];
      const MFloat p2_y = element.m_vertices[1][1];
      // Element length
      geomMeasures[i] = sqrt(pow(p1_x - p2_x, 2) + pow(p1_y - p2_y, 2));
    }
 
    std::array<MFloat, nDim * nDim> vertex;
    m_geometry->elements[i].getVertices(&vertex[0]);
 
    // Calculate surface normal vector
    m_geometry->elements[i].calcNormal(&vertex[0], &normalVecs[i * nDim]);
 
    // Compute dot product of surface normal vector and the vector from the geometric center to the
    // surface centroid
    MFloat dotProduct = 0.0;
    for(MInt k = 0; k < nDim; k++) {
      const MFloat vecCenter = coordinates[i * nDim + k] - geomCenter[k];
      dotProduct += normalVecs[i * nDim + k] * vecCenter;
    }
    // If the dot product is positive the two vectors are pointing in the same 'direction', i.e.,
    // the angle between the vectors is less than 90deg; thus if not >0 change the sign of the
    // normal vector to point outside the geometry
    if(!(dotProduct > 0)) {
      for(MInt k = 0; k < nDim; k++) {
        normalVecs[i * nDim + k] = -normalVecs[i * nDim + k];
      }
    }
  }
 
  const MFloat totalGeomMeasure = std::accumulate(geomMeasures.begin(), geomMeasures.end(), 0.0);
  const MString measure = (nDim == 3) ? "area" : "length";
  m_log << "Total segment " << measure << ": " << totalGeomMeasure << std::endl;
 
  // TODO labels:PP for DEBUG
  std::stringstream dumpFileName;
  dumpFileName << solver().outputDir() << "surface_geometryInfo_" << fileNo << ".dump";
  // DEBUG output for plotting
  FILE* datei;
  datei = fopen(dumpFileName.str().c_str(), "w");
  /* datei = fopen("geometryInfo.dump", "a+"); */
  for(MInt i = 0; i < noElements; i++) {
    fprintf(datei, "%d", i);
 
    for(MInt j = 0; j < nDim; j++) {
      fprintf(datei, " %.8f", coordinates[i * nDim + j]);
    }
 
    for(MInt j = 0; j < nDim; j++) {
      fprintf(datei, " %.8f", normalVecs[i * nDim + j]);
    }
    fprintf(datei, "\n");
  }
  fclose(datei);
 
  const MInt sizeArray1 = noElements * nDim;
  const MInt sizeArray2 = noElements * nDim * nDim;
 
  // Write out geomInfo in output file
  file.defineArray(PIO_FLOAT, "vertices", sizeArray2);
  file.defineArray(PIO_FLOAT, "centroid", sizeArray1);
  file.defineArray(PIO_FLOAT, "normalVector", sizeArray1);
 
  file.defineArray(PIO_FLOAT, "geomMeasure", noElements);
  IF_CONSTEXPR(nDim == 3) { file.setAttribute("area", "description", "geomMeasure"); }
  else {
    file.setAttribute("length", "description", "geomMeasure");
  }
 
  file.setOffset(sizeArray2, 0);
  file.writeArray(&vertices[0], "vertices");
 
  file.setOffset(sizeArray1, 0);
  file.writeArray(&coordinates[0], "centroid");
 
  file.setOffset(sizeArray1, 0);
  file.writeArray(&normalVecs[0], "normalVector");
 
  file.setOffset(noElements, 0);
  file.writeArray(&geomMeasures[0], "geomMeasure");
}
 
 
template <MInt nDim, class SolverType>
void VolumeData<nDim, SolverType>::readAdditionalProperties(const MInt fileNo) {
  TRACE();
 
  const MString volumeShape =
      Context::getSolverProperty<MString>(getBaseName() + "DataShape_" + std::to_string(fileNo), solverId(), AT_);
 
  MInt noInputValuesFloat = -1;
  MInt noInputValuesInt = -1;
 
  if(volumeShape == "box") {
    // Box defined by two points and number of points in each spatial direction
    noInputValuesFloat = 2 * nDim;
    noInputValuesInt = nDim;
  } else if(volumeShape == "cylinder") {
    IF_CONSTEXPR(nDim == 2) { TERMM(1, "Cylinder shape not defined in 2D."); }
    // Cylinder defined by first center point, second center axial position and radius
    noInputValuesFloat = nDim + 2;
    // TODO labels:PP allow to generate non-equidistant points in radial direction to avoid that lots of
    // points are clustered in the center of the cylinder
    // Number of points in each direction: axial, radial, azimuthal; Cylinder axial direction
    // Note: given number of points in radial direction does not include center point
    noInputValuesInt = nDim + 1;
  } else {
    TERMM(1, "Unknown volume data shape: " + volumeShape);
  }
  m_volumeShape.push_back(volumeShape);
 
  if(noInputValuesFloat < 1 || noInputValuesInt < 1) {
    TERMM(1, "Invalid number of input values for volume data shape " + volumeShape + ": "
                 + std::to_string(noInputValuesFloat) + " " + std::to_string(noInputValuesInt));
  }
 
  m_log << "Volume data #" << fileNo << ": " << volumeShape << " with parameters";
 
  // Read float parameters
  std::vector<MFloat> paramFloat(noInputValuesFloat);
  for(MInt i = 0; i < noInputValuesFloat; i++) {
    paramFloat[i] = Context::getSolverProperty<MFloat>(getBaseName() + "DataParameterFloat_" + std::to_string(fileNo),
                                                       solverId(), AT_, i);
    m_log << " " << paramFloat[i];
  }
  m_volumeParameterFloat.push_back(paramFloat);
 
  // Read int parameters
  std::vector<MInt> paramInt(noInputValuesInt);
  for(MInt i = 0; i < noInputValuesInt; i++) {
    paramInt[i] = Context::getSolverProperty<MInt>(getBaseName() + "DataParameterInt_" + std::to_string(fileNo),
                                                   solverId(), AT_, i);
    m_log << " " << paramInt[i];
  }
  m_volumeParameterInt.push_back(paramInt);
 
  m_log << std::endl;
}
 
 
template <MInt nDim, class SolverType>
MInt VolumeData<nDim, SolverType>::noSamplingPoints(const MInt fileNo) {
  TRACE();
 
  MInt noPoints = 0;
  if(m_volumeShape[fileNo] == "box") {
    noPoints = 1;
    for(MInt i = 0; i < nDim; i++) {
      noPoints *= m_volumeParameterInt[fileNo][i];
    }
  } else if(m_volumeShape[fileNo] == "cylinder") {
    // nZ * (nPhi * nR + 1) // +1 for center point
    noPoints =
        m_volumeParameterInt[fileNo][0] * (m_volumeParameterInt[fileNo][1] * m_volumeParameterInt[fileNo][2] + 1);
  } else {
    TERMM(1, "Unknown volume data shape: " + m_volumeShape[fileNo]);
  }
 
  return noPoints;
}
 
 
template <MInt nDim, class SolverType>
void VolumeData<nDim, SolverType>::getSamplingPoint(const MInt fileNo, const MInt pointId, MFloat* const coordinates) {
  TRACE();
 
  if(pointId < 0 || pointId > noSamplingPoints(fileNo)) {
    TERMM(1, "Invalid point id");
  }
 
  // Generate coordinates depending on volume shape
  std::vector<MInt> indexNd(3);
  if(m_volumeShape[fileNo] == "box") {
    // Points for box are created with [x,y,z] = [i,j,k] ordering, i.e., k varies the fastest
    indexXD(pointId, &m_volumeParameterInt[fileNo][0], &indexNd[0]);
 
    for(MInt i = 0; i < nDim; i++) {
      const MFloat minCoord = m_volumeParameterFloat[fileNo][i];
      const MFloat maxCoord = m_volumeParameterFloat[fileNo][nDim + i];
      const MFloat delta = (maxCoord - minCoord) / static_cast<MFloat>(m_volumeParameterInt[fileNo][i] - 1);
      coordinates[i] = minCoord + indexNd[i] * delta;
    }
  } else if(m_volumeShape[fileNo] == "cylinder") {
    /* TERMM(1, "TODO: test volume data sampling with cylinder shape generation"); */
    // ordering: z, r, phi
    // Note: center point (r=0) only appears once
    const MInt nPointsPerCirc = m_volumeParameterInt[fileNo][1] * m_volumeParameterInt[fileNo][2] + 1;
    indexNd[0] = floor(pointId / nPointsPerCirc); // z
 
    const MInt tmpIndex = pointId - indexNd[0] * nPointsPerCirc;
    indexNd[1] = ceil(static_cast<MFloat>(tmpIndex) / static_cast<MFloat>(m_volumeParameterInt[fileNo][2]));
    indexNd[2] = std::max((tmpIndex - 1) % m_volumeParameterInt[fileNo][2], 0);
 
    const MInt dir = m_volumeParameterInt[fileNo][3];
    const MInt dir2 = (dir == 0) ? 1 : 0;
    const MInt dir3 = 3 - dir - dir2;
 
    const MFloat* startPos = &m_volumeParameterFloat[fileNo][0];
    const MFloat nz = static_cast<MFloat>(m_volumeParameterInt[fileNo][0]);
    const MFloat dist_z = static_cast<MFloat>(indexNd[0]) * (startPos[3] - startPos[dir]) / (nz - 1);
 
    const MFloat r = m_volumeParameterFloat[fileNo][4];
    const MFloat nr = m_volumeParameterInt[fileNo][1];
    const MFloat dist_r = static_cast<MFloat>(indexNd[1]) * r / nr;
 
    const MFloat nphi = m_volumeParameterInt[fileNo][2];
    const MFloat phi = static_cast<MFloat>(indexNd[2]) * 2.0 * PI / nphi;
 
    coordinates[dir] = startPos[dir] + dist_z;
    coordinates[dir2] = startPos[dir2] + dist_r * sin(phi);
    coordinates[dir3] = startPos[dir3] + dist_r * cos(phi);
  } else {
    TERMM(1, "Unknown volume data shape: " + m_volumeShape[fileNo]);
  }
}
 
#endif // SAMPLINGDATA_H_