SerialVtuPlot.cpp

#include "VtuPlot.hpp"
#include "Parallel.hpp"
#include "BasicDoFs.hpp"
#include "LagrangeDoF.hpp"
#include "MultiPartDoF.hpp"
#include "DGDoF.hpp"

using namespace std;

SerialVtuPlot::SerialVtuPlot() : VtuPlot() {

}

SerialVtuPlot::SerialVtuPlot(const Mesh &init) : VtuPlot(init) {
  initialize(plotMesh());
}

void SerialVtuPlot::PlotFile(const string &filename) {
  if (!plotMesh.IsInstatiated()) {
    Exit("Nothing to plot")
  }

  if (PPM->master()) {
    createVtu(filename, FileString());
  }
}

void SerialVtuPlot::initPoints(const Mesh &init) {
  pointData["ID"] = new std::unordered_map<Point, PointData>;
  auto pointIDs = pointData["ID"];

  // == Send Points to Master Proc ==
  ExchangeBuffer E;
  for (vertex v = init.vertices(); v != init.vertices_end(); ++v)
    E.Send(0) << v();
  E.Communicate();

  // == Recieve Points and insert into plot ==
  for (short q = 0; q < PPM->size(); ++q) {
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      Point x;
      E.Receive(q) >> x;

      PointData id(1);
      id.AddData(0, pointN);
      pointIDs->insert({x, id});
      points[pointN] = x;
      displacement[pointN++] = Origin;
    }
  }
}

void SerialVtuPlot::initCells(const Mesh &init) {
  ExchangeBuffer E;
  int totalOffset = 0;


  // == Send Cells to Master Proc ==
  for (cell c = init.cells(); c != init.cells_end(); ++c) {
    E.Send(0) << c.ReferenceType();
    E.Send(0) << c();
    E.Send(0) << short(c.Corners());
    for (int i = 0; i < c.Corners(); ++i) E.Send(0) << c[i];
  }
  E.Communicate();

  // == Recieve Cells and insert into plot ==
  for (short q = 0; q < PPM->size(); ++q) {
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      CELLTYPE cType;
      E.Receive(q) >> cType;
      cellTypes += plotLineBreak(cellN, lineBreakZ, 5) + to_string(vtkCellType(cType));

      Point z;
      E.Receive(q) >> z;
      cells[cellN] = z;

      short m;
      E.Receive(q) >> m;

      for (int i = 0; i < m; ++i) {
        E.Receive(q) >> z;
        cellConnectivity +=
            (i == 0 ? plotLineBreak(cellN, lineBreakC, 5) : " ") + to_string(pointID(z));
      }
      totalOffset += m;
      cellOffset += plotLineBreak(cellN++, lineBreakZ, 5) + to_string(totalOffset);
    }
  }
}

void SerialVtuPlot::addPointData(const string &name, const Vector &data) {
  checkMesh(data.GetMesh());

  int dim = data.sizeDoF(0);
  ExchangeBuffer E;
  for(row r = data.rows(); r!= data.rows_end(); ++r){
    if(plotMesh().find_vertex(r())!=plotMesh().vertices_end()){
      E.Send(0) << r();
      for (int j = 0; j < dim; ++j)
        E.Send(0) << data(r, j);
    }
  }
  E.Communicate();

  auto *newPointData = new std::unordered_map<Point, PointData>;
  for (short q = 0; q < PPM->size(); ++q)
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      Point z;
      E.Receive(q) >> z;

      PointData newData(dataDim(dim));
      for (int i = 0; i < dim; i++) {
        double a;
        E.Receive(q) >> a;
        newData.AddData(i, a);
      }
      for (int i = dim; i < dataDim(dim); i++) {
        newData.AddData(i, 0.0);
      }

      newPointData->insert({z, newData});
    }
  pointData[name] = newPointData;
  pointDataSize[name] = dataDim(dim);

  if (activeArrays.first.empty()) activeArrays.first = name;
}

int SerialVtuPlot::dataDim(int dim) const { return (dim > 1 ? 3 : 1); }

void SerialVtuPlot::addCellData(const string &name, const Vector &data) {
  checkMesh(data.GetMesh());

  int dim = data.sizeDoF(0);
  ExchangeBuffer E;
  for(row r = data.rows(); r!= data.rows_end(); ++r){
    if(plotMesh().find_cell(r())!=plotMesh().cells_end()){
      E.Send(0) << r();
      for (int j = 0; j < dim; ++j)
        E.Send(0) << data(r, j);
    }
  }
  E.Communicate();

  auto *newCellData = new std::unordered_map<Point, PointData>;
  for (short q = 0; q < PPM->size(); ++q) {
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      Point z;
      E.Receive(q) >> z;

      PointData newData(dataDim(dim));
      for (int i = 0; i < dim; i++) {
        double a;
        E.Receive(q) >> a;
        newData.AddData(i, a);
      }
      for(int i = dim; i < dataDim(dim); i++){
        newData.AddData(i, 0.0);
      }
      newCellData->insert({z, newData});
    }
  }

  cellData[name] = newCellData;
  cellDataSize[name] = dataDim(dim);

  if (activeArrays.second.empty()) activeArrays.second = name;
}

void SerialVtuPlot::DeformGrid(const Vector &deformation) {
  checkMesh(deformation.GetMesh());

  ExchangeBuffer E;
  int dim = deformation.dim();
  for(row r = deformation.rows(); r!= deformation.rows_end(); ++r){
    if(plotMesh().find_vertex(r())!=plotMesh().vertices_end()){
      E.Send(0) << r();
      for (int j = 0; j < dim; ++j)
        E.Send(0) << deformation(r, j);
    }
  }
  E.Communicate();

  auto *deformationData = new std::unordered_map<Point, Point>;
  for (short q = 0; q < PPM->size(); ++q)
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      Point z;
      E.Receive(q) >> z;

      Point d;
      for (int i = 0; i < dim; i++) {
        E.Receive(q) >> d[i];
      }
      deformationData->insert({z, d});
    }

  for (const auto &gridP : points) {
    displacement[gridP.first] = deformationData->at(gridP.second);
  }

  deformationData->clear();
  delete deformationData;
}

void SerialVtuPlot::AddSubdomain(const Mesh &grid) {
  checkMesh(grid);

  ExchangeBuffer E;
  for (cell c = grid.cells(); c != grid.cells_end(); ++c) {
    E.Send(0) << c();
    E.Send(0) << short(c.Subdomain());
  }
  E.Communicate();

  auto newCellData = new std::unordered_map<Point, PointData>;
  for (short q = 0; q < PPM->size(); ++q) {
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      Point center;
      E.Receive(q) >> center;

      short sd;
      E.Receive(q) >> sd;

      PointData newData(1);
      newData.AddData(0, sd);
      newCellData->insert({center, newData});
    }
  }

  cellData["Subdomain"] = newCellData;
  cellDataSize["Subdomain"] = 1;
}


void SerialVtuPlot::AddProcLoad(const Mesh &grid) {
  checkMesh(grid);

  ExchangeBuffer E;
  for (cell c = grid.cells(); c != grid.cells_end(); ++c) {
    E.Send(0) << c();
    E.Send(0) << PPM->Proc();
  }
  E.Communicate();

  auto newCellData = new std::unordered_map<Point, PointData>;
  for (short q = 0; q < PPM->size(); ++q) {
    while (E.Receive(q).size() < E.ReceiveSize(q)) {
      Point center;
      E.Receive(q) >> center;

      int proc;
      E.Receive(q) >> proc;

      PointData newData(1);
      newData.AddData(0, proc);
      newCellData->insert({center, newData});
    }
  }

  cellData["ProcLoad"] = newCellData;
  cellDataSize["ProcLoad"] = 1;
}

void SerialVtuPlot::initialize(const Mesh &init) {
  initPoints(init);
  initCells(init);

  AddSubdomain(init);
  AddProcLoad(init);
}

void SerialVtuPlot::checkMesh(const Mesh &other) {
  if(!plotMesh.IsInstatiated()){
    plotMesh(other);
    initialize(other);
  }
  else if (&other != &plotMesh()) {
    Exit("Plotting: Meshes don't match!")
  }
}

void SerialVtuPlot::AddData(const string &name, const Vector &data) {
  const DoF& dofToPlot = data.value(0);
  if(typeid(dofToPlot) == typeid(LagrangeDoF) || typeid(dofToPlot) == typeid(MultiPartDoF)){
    int d = ((const LagrangeDoF&) dofToPlot).Degree();
    if(d>0)
      addPointData(name, data);
    else
      addCellData(name, data);
  }
  else if(typeid(dofToPlot) == typeid(CellDoF)) {
    addCellData(name, data);
  }
  else if(typeid(dofToPlot) == typeid(VertexDoF)) {
    addPointData(name, data);
  }
  else if(typeid(dofToPlot) == typeid(DGDoF)) {
    addCellData(name, data);
  }
  else{
    mout << "Warning: Plot type for DoF " << typeid(dofToPlot).name() << " not defined." << endl;
  }

}