SerialVtuPlot.cpp 7.53 KB
Newer Older
1
#include "VtuPlot.hpp"
2
#include "Parallel.hpp"
3
4
#include "BasicDoFs.hpp"
#include "LagrangeDoF.hpp"
5
#include "MultiPartDoF.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
6
#include "DGDoF.hpp"
jonathan.froehlich's avatar
jonathan.froehlich committed
7

8
9
using namespace std;

10
11
12
13
SerialVtuPlot::SerialVtuPlot() : VtuPlot() {

}

jonathan.froehlich's avatar
jonathan.froehlich committed
14
SerialVtuPlot::SerialVtuPlot(const Mesh &init) : VtuPlot(init) {
15
  initialize(plotMesh());
jonathan.froehlich's avatar
jonathan.froehlich committed
16
17
}

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

23
  if (PPM->master()) {
24
    createVtu(filename, FileString());
25
  }
jonathan.froehlich's avatar
jonathan.froehlich committed
26
27
}

jonathan.froehlich's avatar
jonathan.froehlich committed
28
void SerialVtuPlot::initPoints(const Mesh &init) {
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
  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;
48
      displacement[pointN++] = Origin;
jonathan.froehlich's avatar
jonathan.froehlich committed
49
    }
50
  }
jonathan.froehlich's avatar
jonathan.froehlich committed
51
52
53
}

void SerialVtuPlot::initCells(const Mesh &init) {
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
  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;
72
      cellTypes += plotLineBreak(cellN, lineBreakZ, 5) + to_string(vtkCellType(cType));
73
74
75
76
77
78
79
80
81
82
83
84
85
86

      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;
87
      cellOffset += plotLineBreak(cellN++, lineBreakZ, 5) + to_string(totalOffset);
jonathan.froehlich's avatar
jonathan.froehlich committed
88
    }
89
  }
jonathan.froehlich's avatar
jonathan.froehlich committed
90
91
}

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

95
  int dim = data.sizeDoF(0);
96
  ExchangeBuffer E;
97
98
99
100
101
102
  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);
    }
103
104
105
  }
  E.Communicate();

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

112
      PointData newData(dataDim(dim));
113
114
115
116
117
      for (int i = 0; i < dim; i++) {
        double a;
        E.Receive(q) >> a;
        newData.AddData(i, a);
      }
118
      for (int i = dim; i < dataDim(dim); i++) {
119
120
        newData.AddData(i, 0.0);
      }
121

122
      newPointData->insert({z, newData});
jonathan.froehlich's avatar
jonathan.froehlich committed
123
    }
124
  pointData[name] = newPointData;
125
  pointDataSize[name] = dataDim(dim);
jonathan.froehlich's avatar
jonathan.froehlich committed
126

127
  if (activeArrays.first.empty()) activeArrays.first = name;
jonathan.froehlich's avatar
jonathan.froehlich committed
128
129
}

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

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

135
  int dim = data.sizeDoF(0);
136
  ExchangeBuffer E;
137
138
139
140
141
142
  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);
    }
143
144
145
  }
  E.Communicate();

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

152
      PointData newData(dataDim(dim));
153
154
155
156
      for (int i = 0; i < dim; i++) {
        double a;
        E.Receive(q) >> a;
        newData.AddData(i, a);
157
      }
158
      for(int i = dim; i < dataDim(dim); i++){
159
160
        newData.AddData(i, 0.0);
      }
161
      newCellData->insert({z, newData});
jonathan.froehlich's avatar
jonathan.froehlich committed
162
    }
163
  }
jonathan.froehlich's avatar
jonathan.froehlich committed
164

165
  cellData[name] = newCellData;
166
  cellDataSize[name] = dataDim(dim);
jonathan.froehlich's avatar
jonathan.froehlich committed
167

168
  if (activeArrays.second.empty()) activeArrays.second = name;
jonathan.froehlich's avatar
jonathan.froehlich committed
169
170
}

jonathan.froehlich's avatar
jonathan.froehlich committed
171
void SerialVtuPlot::DeformGrid(const Vector &deformation) {
172
  checkMesh(deformation.GetMesh());
173
174
175

  ExchangeBuffer E;
  int dim = deformation.dim();
176
177
178
179
180
181
  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);
    }
182
183
184
  }
  E.Communicate();

185
  auto *deformationData = new std::unordered_map<Point, Point>;
186
187
188
189
190
191
192
193
194
195
  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});
jonathan.froehlich's avatar
jonathan.froehlich committed
196
197
    }

198
199
200
  for (const auto &gridP : points) {
    displacement[gridP.first] = deformationData->at(gridP.second);
  }
jonathan.froehlich's avatar
jonathan.froehlich committed
201

202
203
  deformationData->clear();
  delete deformationData;
jonathan.froehlich's avatar
jonathan.froehlich committed
204
205
206
}

void SerialVtuPlot::AddSubdomain(const Mesh &grid) {
207
  checkMesh(grid);
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

  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});
jonathan.froehlich's avatar
jonathan.froehlich committed
228
    }
229
  }
jonathan.froehlich's avatar
jonathan.froehlich committed
230

231
232
  cellData["Subdomain"] = newCellData;
  cellDataSize["Subdomain"] = 1;
jonathan.froehlich's avatar
jonathan.froehlich committed
233
234
}

235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251

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;

252
253
      int proc;
      E.Receive(q) >> proc;
254
255

      PointData newData(1);
256
      newData.AddData(0, proc);
257
258
259
260
261
262
263
264
265
266
267
      newCellData->insert({center, newData});
    }
  }

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

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

  AddSubdomain(init);
  AddProcLoad(init);
271
272
273
274
275
276
277
278
279
280
281
282
}

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

283
void SerialVtuPlot::AddData(const string &name, const Vector &data) {
284
  const DoF& dofToPlot = data.value(0);
285
  if(typeid(dofToPlot) == typeid(LagrangeDoF) || typeid(dofToPlot) == typeid(MultiPartDoF)){
286
    int d = ((const LagrangeDoF&) dofToPlot).Degree();
287
288
289
290
291
292
293
294
295
296
297
    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);
  }
niklas.baumgarten's avatar
niklas.baumgarten committed
298
299
300
  else if(typeid(dofToPlot) == typeid(DGDoF)) {
    addCellData(name, data);
  }
301
302
303
  else{
    mout << "Warning: Plot type for DoF " << typeid(dofToPlot).name() << " not defined." << endl;
  }
304
305
306

}

307