Commit 12dd1b80 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch 'Refactor-VtuPlot' into 'feature'

Refactor vtu plot

See merge request mpp/mpp!392
parents b1ac1d8e 7336897e
......@@ -73,7 +73,7 @@ mpitest-mpp:
image: ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
script:
- cd /mpp/build
- python3 mppyrun.py --mpi_tests=1 --mute=1
- python3 mppyrun.py --mpi_tests=1 --mute=0
after_script:
- docker rmi ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_MPP}
dependencies: [ "build-mpp" ]
......
......@@ -96,13 +96,8 @@ void Distribution::Communicate(Mesh &mesh) {
vector<std::list<bnd_face>> MarkedBnd(PPM->Size(commSplit));
for (int q = 0; q < PPM->Size(commSplit); ++q) {
for (auto c = markedCells[q].begin(); c != markedCells[q].end(); ++c) {
for (int i = 0; i < c->SpaceCell().Corners(); ++i) {
if (!mesh.IsSTMesh()) {
for (int i = 0; i < c->Corners(); ++i) {
mesh.procSets.Add(c->Corner(i), q);
} else {
mesh.procSets.Add(Point((*c).SpaceCell().Corner(i), (*c).min()), q);
mesh.procSets.Add(Point((*c).SpaceCell().Corner(i), (*c).max()), q);
}
}
for (int i = 0; i < c->SpaceCell().Edges(); ++i)
mesh.procSets.Add(c->SpaceCell().Edge(i), q);
......
......@@ -38,10 +38,8 @@ set(fem_src
transfer/MultiPartTransfer.cpp
transfer/Transfers.cpp
transfer/TransferWrapper.cpp
#plot/ParallelVtuPlot.cpp
plot/Plotting.cpp
plot/Plot.cpp
plot/SerialVtuPlot.cpp
plot/VtuPlot.cpp
plot/VtuPointPlot.cpp
Interface.cpp
......
#include "VtuPlot.hpp"
#include "Config.hpp"
#include "Parallel.hpp"
#include "ParallelVtuPlot.hpp"
#include <fstream>
using namespace std;
ParallelVtuPlot::ParallelVtuPlot(const Mesh &init) : VtuPlot(init) {
pointData["ID"] = new std::unordered_map<Point, PointData>;
auto pointIDs = pointData["ID"];
for (vertex v = init.vertices(); v != init.vertices_end(); v++, pointN++) {
PointData id(1);
id.AddData(0, pointN);
pointIDs->insert({v->first, id});
points[pointN] = v->first;
displacement[pointN] = Point();
}
int totalOffset = 0;
for (cell c = init.cells(); c != init.cells_end(); c++, cellN++) {
cells[cellN] = c();
for (int i = 0; i < c.Corners(); i++) {
cellConnectivity +=
(i == 0 ? plotLineBreak(cellN, lineBreakC, 5) : " ") + to_string(pointID(c.Corner(i)));
}
totalOffset += c.Corners();
cellOffset += plotLineBreak(cellN, lineBreakZ, 5) + to_string(totalOffset) + " ";
cellTypes +=
plotLineBreak(cellN, lineBreakZ, 5) + to_string(vtkCellType(c.ReferenceType())) + " ";
}
}
void ParallelVtuPlot::AddPointData(const string &name, const Vector &data, int dim) {
if (&data.GetMesh() != &plotMesh()) {
THROW("Can't add Data: Meshes don't match!")
}
auto newPointData = new std::unordered_map<Point, PointData>;
for (vertex v = plotMesh().vertices(); v != plotMesh().vertices_end(); v++) {
PointData newData(dim);
for (int i = 0; i < dim; i++)
newData.AddData(i, data(v->first, i));
newPointData->insert({v->first, newData});
}
pointData[name] = newPointData;
pointDataSize[name] = dim;
if (activeArrays.first == "") activeArrays.first = name;
}
void ParallelVtuPlot::AddCellData(const string &name, const Vector &data, int dim) {
if (&data.GetMesh() != &plotMesh()) {
THROW("Can't add Data: Meshes don't match!")
}
auto newCellData = new std::unordered_map<Point, PointData>;
for (cell c = plotMesh().cells(); c != plotMesh().cells_end(); c++) {
PointData newData(dim);
double n = c.Corners();
for (int i = 0; i < dim; i++) {
double d = 0.0;
for (int k = 0; k < n; k++) {
d += data(c.Corner(k), i) / n;
}
newData.AddData(i, d);
}
newCellData->insert({c(), newData});
}
cellData[name] = newCellData;
cellDataSize[name] = dim;
}
void ParallelVtuPlot::PlotFile(const string &filename, bool flush) {
createVtu(filename + string("-Proc") + to_string(PPM->Proc(commSplit())));
if (PPM->Master(commSplit()))
createMasterFile(filename);
if (flush)
ClearData();
}
void ParallelVtuPlot::createMasterFile(const string &filename) {
string filepath = Config::getDataPath() + string("/vtk/") + filename + string(".pvtu");
std::ofstream out(filepath.c_str());
out
<< R"(<VTKFile type ="PUnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">)"
<< endl
<< plotIndent(1) << R"(<PUnstructuredGrid GhostLevel="0">)" << endl;
// Points
out << plotIndent(2) << "<PPoints>" << endl
<< plotIndent(3) << R"(<PDataArray type="Float32" NumberOfComponents="3" format="ascii"/>)"
<< endl
<< plotIndent(2) << "</PPoints>" << endl;
// PointData
out << plotIndent(2) << "<PPointData";
if (!activeArrays.first.empty()) {
string aArray = " ";
aArray += (pointDataSize.at(activeArrays.first) > 1 ? R"(Vectors=")" : R"(Scalars=")");
aArray += activeArrays.first;
aArray += R"(")";
out << aArray;
}
out << ">" << endl;
for (const auto &pData : pointData) {
if (pData.first == "ID") continue;
out << plotIndent(3)
<< R"(<PDataArray type="Float32" Name=")" << pData.first
<< R"(" NumberOfComponents=")" << pointDataSize.at(pData.first)
<< R"("/>)" << endl;
}
out << plotIndent(3) << "</PPointData>" << endl;
// CellData
out << plotIndent(2) << "<PCellData";
if (!activeArrays.second.empty()) {
string aArray = " ";
aArray += (pointDataSize.at(activeArrays.second) > 1 ? R"(Vectors=")" : R"(Scalars=")");
aArray += activeArrays.second;
aArray += R"(")";
out << aArray;
}
out << ">" << endl;
for (const auto &cData : cellData) {
if (cData.first == "ID") continue;
out << plotIndent(3)
<< R"(<PDataArray type="Float32" Name=")" << cData.first
<< R"(" NumberOfComponents=")" << pointDataSize.at(cData.first)
<< R"("/>)" << endl;
}
out << plotIndent(2) << "</PCellData>" << endl;
// Pieces
for (int p = 0; p < PPM->Size(commSplit()); p++) {
out << plotIndent(2) << R"(<Piece Source=")" << filename << "-Proc" << p << ".vtu" << R"("/>)"
<< endl;
}
// End
out << plotIndent(1) << "</PUnstructuredGrid>" << endl
<< "</VTKFile>";
}
void ParallelVtuPlot::DeformGrid(const Vector &deformation) {
if (&deformation.GetMesh() != &plotMesh()) {
THROW("Can't add Data: Meshes don't match!")
}
for (const auto &gridP : points) {
displacement[gridP.first] = deformation(gridP.second);
}
}
void ParallelVtuPlot::AddSubdomain(const Mesh &grid) {
if (&grid != &plotMesh()) {
THROW("Can't add Data: Meshes don't match!")
}
auto newCellData = new std::unordered_map<Point, PointData>;
for (cell c = plotMesh().cells(); c != plotMesh().cells_end(); c++) {
mout << c() << endl;
PointData newData(1);
newData.AddData(0, c.Subdomain());
newCellData->insert({c(), newData});
}
cellData["Subdomain"] = newCellData;
cellDataSize["Subdomain"] = 1;
}
void ParallelVtuPlot::AddProcLoad(const Mesh &grid) {
}
#ifndef PARALLELVTUPLOT_HPP
#define PARALLELVTUPLOT_HPP
#include "VtuPlot.hpp"
// TODO: Reimplement
/*
class ParallelVtuPlot : public VtuPlot {
void createMasterFile(const string &filename);
public:
explicit ParallelVtuPlot(const Mesh &init);
/// Creates a Plotter where the Vector init defines the Grid on which the VtuPlotter operates.
explicit ParallelVtuPlot(const Vector &init) : ParallelVtuPlot(init.GetMesh()) {}
void PlotFile(const string &filename, bool flush) override;
void DeformGrid(const Vector &diplacement) override;
void AddPointData(const string &name, const Vector &data, int dim) override;
void AddCellData(const string &name, const Vector &data, int dim) override;
void AddSubdomain(const Mesh &grid) override;
void AddProcLoad(const Mesh &grid) override;
};
*/
#endif //PARALLELVTUPLOT_HPP
......@@ -32,11 +32,11 @@ void Plotting::Precision(int p) {
}
void Plotting::AddPlot(const string &name) {
plots.try_emplace(name, std::make_unique<SerialVtuPlot>());
plots.try_emplace(name, std::make_unique<VtuPlot>());
}
void Plotting::AddPlot(const std::string &name, const Mesh &mesh) {
plots.try_emplace(name, std::make_unique<SerialVtuPlot>(mesh));
plots.try_emplace(name, std::make_unique<VtuPlot>(mesh));
}
void Plotting::FlushPlot(const std::string &name) {
......
#include "VtuPlot.hpp"
#include "Parallel.hpp"
#include "DGDiscretization.hpp"
#include "DivergenceFreeElement.hpp"
#include "DGVectorFieldElement.hpp"
#if USE_SPACETIME
#include "SpaceTimeDiscretizations.hpp"
#include "STDGDGViscoAcousticElement.hpp"
#endif
using namespace std;
SerialVtuPlot::SerialVtuPlot() : VtuPlot() {
}
SerialVtuPlot::SerialVtuPlot(const Mesh &init) : VtuPlot(init) {
initialize(plotMesh());
}
void SerialVtuPlot::PlotFile(const string &filename) {
if (!plotMesh.IsInstantiated()) {
THROW("Nothing to plot")
}
if (PPM->Master(commSplit())) {
createVtu(filename);
}
}
void SerialVtuPlot::initPoints(const Mesh &init) {
// == Send Points to Master Proc ==
ExchangeBuffer exBuffer(init.CommSplit());
for (vertex v = init.vertices(); v != init.vertices_end(); ++v) exBuffer.Send(0) << v();
exBuffer.Communicate();
// == Receive Points and insert into plot ==
points.clear();
displacements.clear();
for (int q = 0; q < PPM->Size(init.CommSplit()); ++q) {
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point v;
exBuffer.Receive(q) >> v;
PointData id(1);
id.AddData(0, (double) points[v.t()].size());
pointData[v.t()]["ID"].insert({v, id});
points[v.t()].push_back(v);
displacements[v.t()].push_back(Origin);
}
}
}
void SerialVtuPlot::initCells(const Mesh &init) {
ExchangeBuffer exBuffer(init.CommSplit());
std::map<double, int> totalOffset;
for (int i = 0; i < init.steps() - 1; i++)
totalOffset[init.t(i)] = 0;
// == Send Cells to Master Proc ==
for (cell c = init.cells(); c != init.cells_end(); ++c) {
exBuffer.Send(0) << c.ReferenceType();
exBuffer.Send(0) << c();
exBuffer.Send(0) << short(c.Corners());
for (int i = 0; i < c.Corners(); ++i) exBuffer.Send(0) << c[i];
}
exBuffer.Communicate();
// == Receive Cells and insert into plot ==
cells.clear();
for (int q = 0; q < PPM->Size(init.CommSplit()); ++q) {
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
CELLTYPE cType;
exBuffer.Receive(q) >> cType;
Point center;
exBuffer.Receive(q) >> center;
cells[center.t()].push_back(center);
short m;
exBuffer.Receive(q) >> m;
for (int i = 0; i < m; ++i) {
Point corner;
exBuffer.Receive(q) >> corner;
if (corner.t() <= center.t()) {
if (i == 0) {
cellConnectivity[center.t()] +=
plotLineBreak(cells[center.t()].size() - 1, lineBreakC, 5);
cellConnectivity[center.t()] +=
to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
} else {
cellConnectivity[center.t()] +=
" " + to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
}
}
}
cellTypes[center.t()] += plotLineBreak(cells[center.t()].size() - 1, lineBreakZ, 5);
cellTypes[center.t()] += to_string(VtkCellType(cType));
totalOffset[center.t()] += m;
cellOffset[center.t()] += plotLineBreak(cells[center.t()].size() - 1, lineBreakZ, 5);
if (init.IsSTMesh())
cellOffset[center.t()] += to_string(totalOffset[center.t()] / 2);
else
cellOffset[center.t()] += to_string(totalOffset[center.t()]);
}
}
}
void SerialVtuPlot::addPointData(const string &name, const Vector &data,
const std::vector<int> &indices) {
checkMesh(data.GetMesh());
ExchangeBuffer exBuffer(data.CommSplit());
for (row r = data.rows(); r != data.rows_end(); ++r) {
if (plotMesh().find_vertex(r()) != plotMesh().vertices_end()) {
exBuffer.Send(0) << r();
for (int j: indices)
exBuffer.Send(0) << data(r, j);
}
}
communicatePointData(name, exBuffer, indices.size() > 1 ? 3 : 1, indices.size());
}
void SerialVtuPlot::communicatePointData(const string &name, ExchangeBuffer &exBuffer,
int dataDim, int numIndices) {
exBuffer.Communicate();
for (short q = 0; q < PPM->Size(commSplit()); ++q)
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point z;
exBuffer.Receive(q) >> z;
PointData newData(dataDim);
for (int i = 0; i < numIndices; i++) {
double a{};
exBuffer.Receive(q) >> a;
newData.AddData(i, a);
}
for (int i = numIndices; i < dataDim; i++) {
newData.AddData(i, 0.0);
}
pointData[z.t()][name].insert({z, newData});
pointDataSize[z.t()][name] = dataDim;
}
if (activeArrays.first.empty()) activeArrays.first = name;
}
void SerialVtuPlot::addCellData(const string &name, const Vector &data,
const std::vector<int> &indices) {
checkMesh(data.GetMesh());
ExchangeBuffer exBuffer(data.CommSplit());
for (row r = data.rows(); r != data.rows_end(); ++r) {
if (plotMesh().find_cell(r()) != plotMesh().cells_end()) {
exBuffer.Send(0) << r();
for (int j: indices)
exBuffer.Send(0) << data(r, j);
}
}
communicateCellData(name, exBuffer, indices.size() > 1 ? 3 : 1, indices.size());
}
void SerialVtuPlot::communicateCellData(const string &name, ExchangeBuffer &exBuffer,
int dataDim, int numIndices) {
exBuffer.Communicate();
for (short q = 0; q < PPM->Size(commSplit()); ++q) {
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point z;
exBuffer.Receive(q) >> z;
PointData newData(dataDim);
for (int i = 0; i < numIndices; i++) {
double a;
exBuffer.Receive(q) >> a;
newData.AddData(i, a);
}
for (int i = numIndices; i < dataDim; i++) {
newData.AddData(i, 0.0);
}
cellData[z.t()][name].insert({z, newData});
cellDataSize[z.t()][name] = dataDim;
}
}
if (activeArrays.second.empty()) activeArrays.second = name;
}
void SerialVtuPlot::addArgyrisData(const string &name, const Vector &data,
const std::vector<int> &indices) {
checkMesh(data.GetMesh());
std::unordered_map<Point, std::vector<double>> pData{};
for (cell c = data.cells(); c != data.cells_end(); ++c) {
ArgyrisElement E(data, c);
for (int i = 0; i < c.Corners(); ++i) {
const Point &corner = c.Corner(i);
if (pData.find(corner) == pData.end()) {
std::vector<double> cornerData{};
for (int m: indices)
cornerData.push_back(E.Value(c.LocalCorner(i), data, m));
pData[corner] = cornerData;
}
}
}
ExchangeBuffer E;
for (const auto &[point, value]: pData) {
E.Send(0) << point;
for (int j = 0; j < indices.size(); ++j)
E.Send(0) << value[j];
}
communicatePointData(name, E, indices.size() > 1 ? 3 : 1, indices.size());
}
void SerialVtuPlot::addSpaceTimeData(const string &name,
const Vector &data,
const std::vector<int> indices) {
checkMesh(data.GetMesh());
addCellData(name, data, indices);
}
void SerialVtuPlot::addDivergenceFreeData(const string &name, const Vector &data) {
checkMesh(data.GetMesh());
std::unordered_map<Point, Velocity> pData{};
for (cell c = data.cells(); c != data.cells_end(); ++c) {
DivergenceFreeElement E(data, c);
for (int i = 0; i < c.Corners(); ++i) {
const Point &corner = c.Corner(i);
if (pData.find(corner) == pData.end())
pData[corner] = E.VelocityField(c.LocalCorner(i), data);
}
}
ExchangeBuffer E;
for (const auto &[point, value]: pData) {
E.Send(0) << point << value[0] << value[1];
}
communicatePointData(name, E, 3, 2);
}
void SerialVtuPlot::DeformGrid(const Vector &displacement) {
checkMesh(displacement.GetMesh());
const IDiscretization &discToPlot = displacement.GetDisc();
if (typeid(discToPlot) == typeid(DGDiscretization)) {
deformDGGrid(displacement);
} else {
deformLagrangeGrid(displacement);
}
}
void SerialVtuPlot::deformDGGrid(const Vector &displacement) {
ExchangeBuffer exBuffer(displacement.CommSplit());
int dim = displacement.dim();
for (row r = displacement.rows(); r != displacement.rows_end(); ++r) {
auto c = plotMesh().find_cell(r());
if (c != plotMesh().cells_end()) {
DGVectorFieldElement E(displacement, c);
for (int i = 0; i < c.Corners(); ++i) {
exBuffer.Send(0) << c.Corner(i);
VectorField u = E.VectorValue(c.Corner(i), displacement);
for (int j = 0; j < dim; ++j)
exBuffer.Send(0) << u[j];
}
}
}
exBuffer.Communicate();
unordered_map<Point, Point> deformationData;
unordered_map<Point, int> deformationCount;
for (short q = 0; q < PPM->Size(displacement.CommSplit()); ++q)
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point z;
exBuffer.Receive(q) >> z;
Point d;
for (int i = 0; i < dim; i++) {
exBuffer.Receive(q) >> d[i];
}
if (deformationData.find(z) != deformationData.end()) {
deformationData[z] += d;
deformationCount[z] += 1;
} else {
deformationData.insert({z, d});
deformationCount.insert({z, 1});
}
}
// for (auto &p: points) {
for (int i = 0; i < points.size(); i++) {
displacements[0.0][i] =
deformationData[points[0.0][i]] / deformationCount[points[0.0][i]];
}
// }
deformationData.clear();
}
void SerialVtuPlot::deformLagrangeGrid(const Vector &displacement) {
ExchangeBuffer exBuffer(displacement.CommSplit());
int dim = displacement.dim();
for (row r = displacement.rows(); r != displacement.rows_end(); ++r) {
if (plotMesh().find_vertex(r()) != plotMesh().vertices_end()) {
exBuffer.Send(0) << r();
for (int j = 0; j < dim; ++j)
exBuffer.Send(0) << displacement(r, j);
}
}
exBuffer.Communicate();
auto *deformationData = new unordered_map<Point, Point>;
for (short q = 0; q < PPM->Size(displacement.CommSplit()); ++q)
while (exBuffer.Receive(q).size() < exBuffer.ReceiveSize(q)) {
Point z;
exBuffer.Receive(q) >> z;
Point d;
for (int i = 0; i < dim; i++) {
exBuffer.Receive(q) >> d[i];
}
deformationData->insert({z, d});
}