CoarseGeometry.cpp

#include "M_IOFiles.hpp"
#include "CoarseGeometry.hpp"


M_ifstream CoarseGeometry::loadGeometryFile(const std::string &name) {
  std::string GeoPath = Config::getSearchPath() + "conf/geo/";
  config.get<std::string>("GeoPath", GeoPath);
  std::string GeoName = name + ".geo";
  if (!FileExists(GeoName))
    GeoName = GeoPath + GeoName;
  if (!FileExists(GeoName)) Exit("no mesh with name: " + name);
  return {GeoName.c_str()};
}

bool CoarseGeometry::readGeometryFile(std::istream &geoFileStream) {
  char L[len];
  geoFileStream.getline(L, len);
  if (strncasecmp("points", L, 6) == 0) geoFileStream.getline(L, len);
  while (insertCoordinate(L, coordinates)) geoFileStream.getline(L, len);
  if (strncasecmp("cells", L, 5) == 0) geoFileStream.getline(L, len);
  while (insertID(L, cellIds)) geoFileStream.getline(L, len);
  if (strncasecmp("faces", L, 5) == 0) geoFileStream.getline(L, len);
  while (insertID(L, faceIds)) geoFileStream.getline(L, len);
  if (strncasecmp("vdata", L, 5) == 0) geoFileStream.getline(L, len);
  while (insertData(L, vdataList)) geoFileStream.getline(L, len);
  if (strncasecmp("cdata", L, 5) == 0) geoFileStream.getline(L, len);
  while (insertData(L, cdataList)) geoFileStream.getline(L, len);
  if (strncasecmp("times", L, 5) == 0) {
    geoFileStream.getline(L, len);
    insertTimeSteps(L, timeSteps);
  }
  return true;
}

void CoarseGeometry::initIds(const vector<vector<int>> &values, vector<Ids> &ids) {
  for (const auto &value : values) {
    if (value.size() < 3) continue;
    int n = value[0];
    int f = value[1];
    int tmp[27];
    for (int j = 2; j < value.size(); ++j)
      tmp[j - 2] = value[j];
    ids.emplace_back(Ids(n, value.size() - 2, f, tmp));
  }
}

bool CoarseGeometry::insertCoordinate(char *L, vector<Point> &coordinateList) {
  dout(100) << L << "\n";
  double z[4] = {0, 0, 0, 0};
  int d = sscanf(L, "%lf %lf %lf %lf", z, z + 1, z + 2, z + 3);
  if (d < 1) return false;
  coordinates.emplace_back(Point(z));
  return true;
}

void CoarseGeometry::insertTimeSteps(char *L, std::vector<double> &timeSteps) {
  std::string content(L);
  std::vector<string> ts_strings = split(content, " ");
  ts_strings = ft::filter(ts_strings, [](auto &s) { return !s.empty(); });
  timeSteps.resize(ts_strings.size(), 0.0);
  std::transform(ts_strings.begin(), ts_strings.end(), timeSteps.begin(), ft::stod);
}

std::vector<double> CoarseGeometry::refineTimeSteps(const std::vector<double> &timeSteps) {
  // refining scale times the time steps
  Exit("Not tested!");
  vector<vector<double>> init_timesteps(refineTimeStepCount + 1);
  init_timesteps[0] = timeSteps;
  int SF_inv = 2.0;
  for (unsigned int S = 1; S < refineTimeStepCount + 1; ++S) {
    unsigned int m = init_timesteps[S - 1].size();
    init_timesteps[S].resize(2 * m - 1);
    init_timesteps[S][0] = init_timesteps[S - 1][0];
    init_timesteps[S][m + (m - 1) * (SF_inv - 1) - 1] = init_timesteps[S - 1][m - 1];
    for (unsigned int k = 0; k < m - 1; ++k) {
      init_timesteps[S][k * SF_inv + 1] = init_timesteps[S - 1][k + 1];
      double dt = init_timesteps[S - 1][k + 2] - init_timesteps[S - 1][k + 1];
      for (unsigned int l = 1; l < SF_inv; ++l)
        init_timesteps[S][k * SF_inv + l + 1] =
            init_timesteps[S - 1][k + 1] + dt * 1 / double(SF_inv) * l;
    }
  }
  return init_timesteps[refineTimeStepCount];
}

bool CoarseGeometry::insertID(char *L, vector<Ids> &ids) {
  int n, f, pos;
  int m = sscanf(L, "%d %d %n", &n, &f, &pos);

  if (m < 2)
    return false;
  int cornerIndices[n];
  for (int i = 0; i < n; i++) {
    int cur;
    assert(sscanf(L + pos, "%d %n", &cornerIndices[i], &cur) == 1);
    pos += cur;
  }
  ids.emplace_back(Ids(n, n, f, cornerIndices));
  return true;
}

bool CoarseGeometry::insertData(char *L, vector<DataContainer> &vertexData) {
  int n, pos;
  int m = sscanf(L, "%d %n", &n, &pos);

  if (m < 1 || n < 1)
    return false;
  DataContainer pContainer(n);
  for (int i = 0; i < n; i++) {
    double d;
    int cur;

    assert(sscanf(L + pos, "%lf %n", &d, &cur) == 1);
    pContainer.Set(i, d);
    pos += cur;
  }
  vertexData.emplace_back(pContainer);
  return true;
}

CoarseGeometry::CoarseGeometry(const string &geometryFilename) {
  if (PPM->Master()) {
    M_ifstream geoFileStream = loadGeometryFile(geometryFilename);
    readGeometryFile(geoFileStream);
    if (scalingFactor != 1.0) {
      mout << "  scaling geometry by " << scalingFactor << endl;
      Scale(scalingFactor);
    }
  }
  CommunicateGeometry();
  geoName = geometryFilename;
}

void CoarseGeometry::CommunicateGeometry() {
  if (PPM->Size() != 1) {
    ExchangeBuffer exBuffer;
    if (PPM->Master()) {
      for (int q = 1; q < PPM->Size(); ++q) {
        exBuffer.Send(q) << int(timeSteps.size());
        for (double timeStep : timeSteps) {
          exBuffer.Send(q) << double(timeStep);
        }
      }
    }
    exBuffer.Communicate();
    if (!PPM->Master()) {
      int m;
      exBuffer.Receive(0) >> m;
      timeSteps.resize(m);
      for (int i = 0; i < m; i++) {
        exBuffer.Receive(0) >> timeSteps[i];

      }
    }
    exBuffer.ClearBuffers();
  }
};

CoarseGeometry::CoarseGeometry(const MeshPoints &points, const MeshCells &cells,
                               const MeshFaces &faces, const TimeSteps &timeSteps) :
    timeSteps(timeSteps) {
  for (const auto &point : points) {
    double z[]{0.0, 0.0, 0.0, 0.0};
    for (int j = 0; j < min(3, int(point.size())); ++j)
      z[j] = point[j];
    coordinates.emplace_back(Point(z));
  }
  initIds(cells, cellIds);
  initIds(faces, faceIds);
}

void CoarseGeometry::Scale(double scalingFactor) {
  for (auto &z : coordinates)
    z *= scalingFactor;
}

Logging &operator<<(Logging &s, const CoarseGeometry &M) {
  return s << "POINTS: " << M.coordinates.size() << endl << M.coordinates
           << "CELLS: " << M.cellIds.size() << endl << M.cellIds
           << "FACES: " << M.faceIds.size() << endl << M.faceIds
           #ifdef USE_DATAMESH
           << "VDATA: " << M.vdataList.size() << endl << M.vdataList
           << "CDATA: " << M.cdataList.size() << endl << M.cdataList
           #endif
           << endl;
}

std::string CoarseGeometry::Name() const {
  return geoName;
}