DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Mesh_Plts.cpp
Go to the documentation of this file.
1#include "Geom/Quadrature.hpp"
2#include "Mesh.hpp"
3#include "CGNS.hpp"
4#include <thread>
5#include <filesystem>
6#include <base64_rfc4648.hpp>
7#include <zlib.h>
8
9#include <nanoflann.hpp>
10#include "PointCloud.hpp"
11
12#include <hdf5.h>
13#include <cgnslib.h>
14#include <pcgnslib.h>
15
16namespace DNDS::Geom
17{
18
19 void UnstructuredMeshSerialRW::
20 GetCurrentOutputArrays(int flag,
21 tCoordPair &coordOut,
22 tAdjPair &cell2nodeOut,
23 tPbiPair &cell2nodePbiOut,
24 tElemInfoArrayPair &cellElemInfoOut,
25 index &nCell, index &nNode)
26 {
27 /*************************************************/
28 // get the output arrays: serial or parallel
29 tCoord coordSerialDummy;
30 tAdj cell2nodeSerialDummy;
31 tPbi cell2nodePbiSerialDummy;
32 tElemInfoArray cellElemInfoSerialDummy;
33 coordSerialDummy = make_ssp<decltype(coordSerialDummy)::element_type>(ObjName{"GetCurrentOutputArrays::coordSerialDummy"}, mesh->getMPI());
34 cell2nodeSerialDummy = make_ssp<decltype(cell2nodeSerialDummy)::element_type>(ObjName{"GetCurrentOutputArrays::cell2nodeSerialDummy"}, mesh->getMPI());
35 cell2nodePbiSerialDummy = make_ssp<decltype(cell2nodePbiSerialDummy)::element_type>(ObjName{"GetCurrentOutputArrays::cell2nodePbiSerialDummy"}, NodePeriodicBits::CommType(), NodePeriodicBits::CommMult(), mesh->getMPI());
36 cellElemInfoSerialDummy = make_ssp<decltype(cellElemInfoSerialDummy)::element_type>(ObjName{"GetCurrentOutputArrays::cellElemInfoSerialDummy"}, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
37
38 if (flag == 0)
39 {
40 coordOut.father = coordSerial;
41 coordOut.son = coordSerialDummy;
42 cell2nodeOut.father = cell2nodeSerial;
43 cell2nodeOut.son = cell2nodeSerialDummy;
44 if (mesh->isPeriodic)
45 {
46 cell2nodePbiOut.father = cell2nodePbiSerial;
47 cell2nodePbiOut.son = cell2nodePbiSerialDummy;
48 }
49 cellElemInfoOut.father = cellElemInfoSerial;
50 cellElemInfoOut.son = cellElemInfoSerialDummy;
51 nCell = cell2nodeOut.Size();
52 nNode = coordOut.Size();
53 }
54 else if (flag == 1)
55 {
56 coordOut.father = mesh->coords.father;
57 coordOut.son = mesh->coords.son;
58 cell2nodeOut.father = mesh->cell2node.father;
59 cell2nodeOut.son = mesh->cell2node.son;
60 if (mesh->isPeriodic)
61 {
62 cell2nodePbiOut.father = mesh->cell2nodePbi.father;
63 cell2nodePbiOut.son = mesh->cell2nodePbi.son;
64 }
65 cellElemInfoOut.father = mesh->cellElemInfo.father;
66 cellElemInfoOut.son = mesh->cellElemInfo.son;
67 nCell = cell2nodeOut.father->Size(); //! only non-ghost cells are output
68 nNode = coordOut.Size(); //! need all the nodes with ghost
69 }
70 /*************************************************/
71 }
72
73 static void GetViolentNodeDeduplication(
74 index nNode, const std::vector<tPoint> &nodesExtra, tCoordPair &coordOut,
75 std::vector<index> &nodeDedu2Old, std::vector<index> &nodeOld2Dedu)
76 {
77 PointCloudFunctional nodesCloud(
78 [&](size_t idx) -> tPoint
79 {
80 if (idx < nNode)
81 return coordOut[idx];
82 else
83 return nodesExtra.at(idx - nNode);
84 },
85 nNode + nodesExtra.size());
86 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
87 nanoflann::L2_Simple_Adaptor<real, PointCloudFunctional>,
88 PointCloudFunctional,
89 3,
90 index>;
91 kdtree_t kdTree(3, nodesCloud);
92 nodeOld2Dedu.resize(nodesCloud.size(), -1);
93 index iNodeDeduTop{0};
94 for (index i = 0; i < nodesCloud.size(); i++)
95 {
96 if (nodeOld2Dedu[i] != -1)
97 continue;
98 tPoint cn = nodesCloud[i];
99#if NANOFLANN_VERSION < 0x150
100 std::vector<std::pair<DNDS::index, DNDS::real>> IndicesDists;
101#else
102 std::vector<nanoflann::ResultItem<DNDS::index, DNDS::real>> IndicesDists;
103#endif
104 IndicesDists.reserve(5);
105#if NANOFLANN_VERSION < 0x150
106 nanoflann::SearchParams params{}; // default params
107#else
108 nanoflann::SearchParameters params{};
109#endif
110 index nFound = kdTree.radiusSearch(cn.data(), 1e-20, IndicesDists, params);
111 int nSelf{0};
112 for (auto &v : IndicesDists)
113 {
114 if (v.first == i)
115 nSelf++;
116 nodeOld2Dedu.at(v.first) = iNodeDeduTop;
117 }
118 iNodeDeduTop++;
119 DNDS_assert(nSelf == 1);
120 }
121 nodeDedu2Old.resize(iNodeDeduTop, -1);
122 for (index i = 0; i < nodeOld2Dedu.size(); i++)
123 nodeDedu2Old.at(nodeOld2Dedu.at(i)) = i;
124 }
125
126 void UnstructuredMeshSerialRW::
127 PrintSerialPartPltBinaryDataArray(std::string fname,
128 int arraySiz, int arraySizPoint,
129 const tFGetName &names,
130 const tFGetData &data,
131 const tFGetName &namesPoint,
132 const tFGetData &dataPoint,
133 double t, int flag)
134 {
135 auto mpi = mesh->getMPI();
136
137 if (mpi.rank != mRank && flag == 0) //* now only operating on mRank if serial
138 return;
139
140 if (flag == 0)
141 {
142 fname += ".plt";
143 std::filesystem::path outFile{fname};
144 std::filesystem::create_directories(outFile.parent_path() / ".");
145 DNDS_assert(mode == SerialOutput && dataIsSerialOut);
146 }
147 if (flag == 1)
148 {
149
150 std::filesystem::path outPath{fname + ".dir"};
151 std::filesystem::create_directories(outPath);
152 char BUF[512];
153 std::sprintf(BUF, "%06d", mpi.rank);
154 fname = getStringForcePath(outPath / (std::string(BUF) + ".plt"));
155 }
156
157 std::ofstream fout(fname, std::ios::binary);
158 if (!fout)
159 {
160 DNDS::log() << "Error: WriteMeshDebugTecASCII open \"" << fname << "\" failure" << std::endl;
161 DNDS_assert(false);
162 }
163 const char magic_word[] = "#!TDV112";
164 const int b_magic_word = sizeof(magic_word) - 1;
165 int32_t intBuf;
166 double_t doubleBuf;
167 float_t floatBuf;
168
169 auto writeInt = [&](int d) -> void
170 {
171 intBuf = d;
172 fout.write((char *)(&intBuf), sizeof(intBuf));
173 };
174 auto writeFloat = [&](float_t d) -> void
175 {
176 floatBuf = d;
177 fout.write((char *)(&floatBuf), sizeof(floatBuf));
178 };
179 auto writeDouble = [&](double_t d) -> void
180 {
181 doubleBuf = d;
182 fout.write((char *)(&doubleBuf), sizeof(doubleBuf));
183 };
184 auto writeString = [&](const std::string &title) -> void
185 {
186 for (auto i : title)
187 {
188 intBuf = i;
189 fout.write((char *)(&intBuf), sizeof(intBuf));
190 }
191 intBuf = 0;
192 fout.write((char *)(&intBuf), sizeof(intBuf));
193 };
194 fout.write(magic_word, b_magic_word);
195 writeInt(1);
196 writeInt(0); //! full: write both grid and data
197 writeString("Title Here");
198 writeInt(arraySiz + 3 + 1 + arraySizPoint); // nvars
199 writeString("X");
200 writeString("Y");
201 writeString("Z");
202 for (int idata = 0; idata < arraySizPoint; idata++)
203 writeString(namesPoint(idata));
204 for (int idata = 0; idata < arraySiz; idata++)
205 writeString(names(idata));
206 writeString("iPart");
207
208 /********************************/
209 // cellZone in header
210 /********************************/
211 writeFloat(299.0F); // 299.0 indicates a v112 zone header, available in v191
212 writeString("zone_0");
213 writeInt(-1); // ParentZone: No longer used.
214 writeInt(-1); // StrandID: static strand ID
215 writeDouble(t); // solution time
216 writeInt(-1); // default zone color
217 if (mesh->dim == 2)
218 writeInt(3); // 2d: quad zone
219 else if (mesh->dim == 3)
220 writeInt(5); // 3d: brick zone
221 else if (mesh->dim == 1)
222 writeInt(1); // 1d: line zone
223 writeInt(1); // specifyVarLocation
224
225 for (int idim = 0; idim < 3; idim++)
226 writeInt(0); // xyz at node
227 for (int idata = 0; idata < arraySizPoint; idata++)
228 writeInt(0); // data at point
229 for (int idata = 0; idata < arraySiz; idata++)
230 writeInt(1); // data at center
231 writeInt(1); // iPart
232
233 writeInt(0); // Are raw local 1-to-1 face neighbors supplied?
234 writeInt(0); // Number of miscellaneous user-defined face
235
236 /*******************************************************/
237 // output preparation
238 tCoordPair coordOut;
239 tAdjPair cell2nodeOut;
240 tPbiPair cell2nodePbiOut;
241 tElemInfoArrayPair cellElemInfoOut;
242 index nCell{-1}, nNode{-1};
243 this->GetCurrentOutputArrays(flag, coordOut, cell2nodeOut, cell2nodePbiOut, cellElemInfoOut,
244 nCell, nNode);
245
246 std::vector<Geom::tPoint>
247 nodesExtra;
248 std::vector<index> nodesExtraAtOriginal;
249 if (mesh->isPeriodic)
250 {
251 for (index iCell = 0; iCell < nCell; iCell++)
252 for (rowsize ic2n = 0; ic2n < cell2nodePbiOut.RowSize(iCell); ic2n++)
253 if (cell2nodePbiOut[iCell][ic2n])
254 nodesExtra.push_back(mesh->periodicInfo.GetCoordByBits(coordOut[cell2nodeOut[iCell][ic2n]], cell2nodePbiOut[iCell][ic2n])),
255 nodesExtraAtOriginal.push_back(cell2nodeOut(iCell, ic2n));
256 }
257 // if (mpi.rank == mRank)
258 // std::cout << "PrintSerialPartPltBinaryDataArray === " << std::endl;
259 // auto printSize = [&]()
260 // { std::cout << "Rank [" << mpi.rank << "]"
261 // << " Size: " << nNode << " " << nCell << " " << nodesExtra.size() << std::endl; };
262 // if (flag == 1)
263 // MPISerialDo(mpi, printSize);
264 // else
265 // printSize();
266 std::vector<index> nodeDedu2Old;
267 std::vector<index> nodeOld2Dedu;
268 GetViolentNodeDeduplication(nNode, nodesExtra, coordOut, nodeDedu2Old, nodeOld2Dedu);
269 /*******************************************************/
270
271 writeInt(nodeDedu2Old.size()); // node number
272 writeInt(nCell); // cell number
273
274 writeInt(0); // I dim
275 writeInt(0); // J dim
276 writeInt(0); // K dim
277 writeInt(0); // No more Auxiliary name/value pairs.
278
279 writeFloat(357.0F); // end of header, EOH marker
280
281 writeFloat(299.0F); // 299.0 indicates a v112 zone header, available in v191
282
283 /********************************/
284 // cellZone data
285 /********************************/
286
287 for (int idim = 0; idim < 3; idim++)
288 writeInt(2); // double for node
289 for (int idata = 0; idata < arraySizPoint; idata++)
290 writeInt(2); // double for data pint
291 for (int idata = 0; idata < arraySiz; idata++)
292 writeInt(2); // double for data
293 writeInt(2); // double for iPart
294
295 writeInt(0); // no passive
296 writeInt(0); // no sharing
297 writeInt(-1); // no sharing
298
299 std::vector<double_t> minVal(3 + arraySiz, DNDS::veryLargeReal);
300 std::vector<double_t> maxVal(3 + arraySiz, -DNDS::veryLargeReal); // for all non-shared non-passive
301 std::vector<double_t> minValPoint(arraySizPoint, DNDS::veryLargeReal);
302 std::vector<double_t> maxValPoint(arraySizPoint, -DNDS::veryLargeReal); //! Tecplot is sensitive to the correctness of min/max val
303 for (int idim = 0; idim < 3; idim++)
304 for (DNDS::index i = 0; i < nodeDedu2Old.size(); i++)
305 {
306 index iN = nodeDedu2Old.at(i);
307 if (iN < nNode)
308 {
309 minVal[idim] = std::min(coordOut[iN](idim), minVal[idim]);
310 maxVal[idim] = std::max(coordOut[iN](idim), maxVal[idim]);
311 }
312 else
313 {
314 minVal[idim] = std::min(nodesExtra.at(iN - nNode)(idim), minVal[idim]);
315 maxVal[idim] = std::max(nodesExtra.at(iN - nNode)(idim), maxVal[idim]);
316 }
317 }
318
319 for (int idata = 0; idata < arraySiz; idata++)
320 for (DNDS::index iv = 0; iv < nCell; iv++)
321 {
322 minVal[3 + idata] = std::min(data(idata, iv), minVal[3 + idata]);
323 maxVal[3 + idata] = std::max(data(idata, iv), maxVal[3 + idata]);
324 }
325
326 for (int idata = 0; idata < arraySizPoint; idata++)
327 for (DNDS::index iv = 0; iv < nNode; iv++)
328 {
329 minValPoint[idata] = std::min(dataPoint(idata, iv), minValPoint[idata]);
330 maxValPoint[idata] = std::max(dataPoint(idata, iv), maxValPoint[idata]);
331 }
332
333 for (int idim = 0; idim < 3; idim++)
334 {
335 writeDouble(minVal[idim]);
336 writeDouble(maxVal[idim]);
337 }
338 for (int idata = 0; idata < arraySizPoint; idata++)
339 {
340 writeDouble(minValPoint[idata]);
341 writeDouble(maxValPoint[idata]);
342 }
343 for (int idata = 0; idata < arraySiz; idata++)
344 {
345 writeDouble(minVal[3 + idata]);
346 writeDouble(maxVal[3 + idata]);
347 }
348 writeDouble(0);
349 writeDouble(mpi.size);
350
351 for (int idim = 0; idim < 3; idim++)
352 {
353 for (DNDS::index i = 0; i < nodeDedu2Old.size(); i++)
354 {
355 DNDS::index iN = nodeDedu2Old.at(i);
356 if (iN < nNode)
357 writeDouble(coordOut[iN](idim));
358 else
359 writeDouble(nodesExtra.at(iN - nNode)(idim));
360 }
361 }
362
363 for (int idata = 0; idata < arraySizPoint; idata++)
364 {
365 for (DNDS::index i = 0; i < nodeDedu2Old.size(); i++)
366 {
367 DNDS::index iN = nodeDedu2Old.at(i);
368 if (iN < nNode)
369 writeDouble(dataPoint(idata, iN));
370 else
371 writeDouble(dataPoint(idata, nodesExtraAtOriginal.at(iN - nNode)));
372 }
373 }
374
375 for (int idata = 0; idata < arraySiz; idata++)
376 for (DNDS::index iv = 0; iv < nCell; iv++)
377 {
378 writeDouble(data(idata, iv));
379 }
380
381 for (DNDS::index iv = 0; iv < nCell; iv++)
382 {
383 DNDS::MPI_int r = -1;
384 DNDS::index v = -1;
385 if (flag == 0)
386 {
387 bool ret = cell2nodeSerialOutTrans.pLGlobalMapping->search(iv, r, v);
388 DNDS_assert_info(ret, "search failed");
389 }
390 else if (flag == 1)
391 {
392 r = mesh->getMPI().rank;
393 // if (mesh->localPartitionStarts.size()) // debug: for inspection of sub-partitions
394 // {
395 // r *= 10000;
396 // auto &lps = mesh->localPartitionStarts;
397 // r += std::upper_bound(lps.begin(), lps.end(), iv) - lps.begin();
398 // }
399 }
400 writeDouble(r);
401 }
402
403 index nExtra{0};
404
405 for (DNDS::index iv = 0; iv < nCell; iv++)
406 {
407 auto elem = Elem::Element{cellElemInfoOut[iv]->getElemType()};
408 std::vector<index> c2n = cell2nodeOut[iv];
409 if (mesh->isPeriodic)
410 for (rowsize ic2n = 0; ic2n < c2n.size(); ic2n++)
411 {
412 if (cell2nodePbiOut[iv][ic2n])
413 c2n[ic2n] = (nExtra++) + nNode;
414 }
415
416 for (auto &v : c2n)
417 v = nodeOld2Dedu.at(v);
418 switch (elem.GetParamSpace())
419 {
420 case Elem::ParamSpace::LineSpace:
421 writeInt(c2n[0] + 0LL);
422 writeInt(c2n[1] + 0LL);
423 break;
424 case Elem::ParamSpace::TriSpace:
425 writeInt(c2n[0] + 0LL);
426 writeInt(c2n[1] + 0LL);
427 writeInt(c2n[2] + 0LL);
428 writeInt(c2n[2] + 0LL);
429 break;
430 case Elem::ParamSpace::QuadSpace:
431 writeInt(c2n[0] + 0LL);
432 writeInt(c2n[1] + 0LL);
433 writeInt(c2n[2] + 0LL);
434 writeInt(c2n[3] + 0LL); // ! note that tis is zero based
435 break;
436 case Elem::ParamSpace::TetSpace:
437 writeInt(c2n[0] + 0LL);
438 writeInt(c2n[1] + 0LL);
439 writeInt(c2n[2] + 0LL);
440 writeInt(c2n[2] + 0LL);
441 writeInt(c2n[3] + 0LL);
442 writeInt(c2n[3] + 0LL);
443 writeInt(c2n[3] + 0LL);
444 writeInt(c2n[3] + 0LL);
445 break;
446 case Elem::ParamSpace::HexSpace:
447 writeInt(c2n[0] + 0LL);
448 writeInt(c2n[1] + 0LL);
449 writeInt(c2n[2] + 0LL);
450 writeInt(c2n[3] + 0LL);
451 writeInt(c2n[4] + 0LL);
452 writeInt(c2n[5] + 0LL);
453 writeInt(c2n[6] + 0LL);
454 writeInt(c2n[7] + 0LL);
455 break;
456 case Elem::ParamSpace::PrismSpace:
457 writeInt(c2n[0] + 0LL);
458 writeInt(c2n[1] + 0LL);
459 writeInt(c2n[2] + 0LL);
460 writeInt(c2n[2] + 0LL);
461 writeInt(c2n[3] + 0LL);
462 writeInt(c2n[4] + 0LL);
463 writeInt(c2n[5] + 0LL);
464 writeInt(c2n[5] + 0LL);
465 break;
466 case Elem::ParamSpace::PyramidSpace:
467 writeInt(c2n[0] + 0LL);
468 writeInt(c2n[1] + 0LL);
469 writeInt(c2n[2] + 0LL);
470 writeInt(c2n[3] + 0LL);
471 writeInt(c2n[4] + 0LL);
472 writeInt(c2n[4] + 0LL);
473 writeInt(c2n[4] + 0LL);
474 writeInt(c2n[4] + 0LL);
475 break;
476 default:
477 DNDS_assert(false); //! 2d or 3d elems
478 }
479 };
480 fout.close();
481 }
482
483 static void updateVTKSeries(std::string seriesName, std::string fname, real tSimu)
484 {
485 using json = nlohmann::json;
486 json j;
487 if (std::filesystem::exists(seriesName))
488 {
489 std::ifstream fin(seriesName);
490 fin >> j;
491 fin.close();
492 DNDS_assert(j["files"].is_array());
493 }
494 else
495 {
496 j["file-series-version"] = "1.0";
497 j["files"] = json::array();
498 }
499 std::map<std::string, real> series;
500 for (auto &f : j["files"])
501 series[f["name"].get<std::string>()] = f["time"].get<real>();
502 series[fname] = tSimu;
503 j["files"].clear();
504 for (auto &[fn, tS] : series)
505 {
506 j["files"].emplace_back(json::object());
507 j["files"].back()["name"] = fn;
508 j["files"].back()["time"] = tS;
509 }
510 std::ofstream fout(seriesName);
511 fout << j.dump(4);
512 fout.close();
513 }
514
515 // =================================================================
516 // File-local helpers for PrintSerialPartVTKDataArray
517 // =================================================================
518 namespace
519 {
520 /**
521 * \brief Result of VTK cell topology construction.
522 */
523 struct VTKCellTopology
524 {
525 std::vector<int64_t> connectivity; ///< Flat node-index list for all cells
526 std::vector<int64_t> offsets; ///< Per-cell cumulative offset into connectivity
527 std::vector<uint8_t> cellTypes; ///< VTK cell type per cell
528 };
529
530 /**
531 * \brief Build VTK connectivity, offsets, and cell-type arrays from
532 * the output mesh arrays.
533 *
534 * For periodic meshes, replaces periodic-boundary node references with
535 * newly created "extra" node indices, then applies node deduplication.
536 * Converts element node ordering to VTK convention.
537 *
538 * \param[in] nCell Number of cells.
539 * \param[in] nNode Number of original (non-extra) nodes.
540 * \param[in] cell2nodeOut Cell-to-node adjacency.
541 * \param[in] cell2nodePbiOut Cell-to-node periodic bits.
542 * \param[in] cellElemInfoOut Cell element info (for element type).
543 * \param[in] nodeOld2Dedu Old-to-deduplicated node index mapping.
544 * \param[in] isPeriodic Whether periodic mesh handling is enabled.
545 */
546 VTKCellTopology BuildVTKCellTopology(
547 DNDS::index nCell,
548 DNDS::index nNode,
549 const tAdjPair &cell2nodeOut,
550 const tPbiPair &cell2nodePbiOut,
551 const tElemInfoArrayPair &cellElemInfoOut,
552 const std::vector<DNDS::index> &nodeOld2Dedu,
553 bool isPeriodic)
554 {
555 VTKCellTopology result;
556 result.connectivity.reserve(cell2nodeOut.father->DataSize() * 2);
557 result.offsets.resize(nCell);
558 result.cellTypes.resize(nCell);
559 DNDS::index cellEnd = 0;
560 DNDS::index nNodeExtra{0};
561 for (DNDS::index iCell = 0; iCell < nCell; iCell++)
562 {
563 auto elem = Elem::Element{cellElemInfoOut[iCell]->getElemType()};
564 std::vector<DNDS::index> c2n = cell2nodeOut[iCell];
565 if (isPeriodic) // alter the pointing
566 for (DNDS::rowsize ic2n = 0; ic2n < c2n.size(); ic2n++)
567 if (cell2nodePbiOut(iCell, ic2n))
568 c2n[ic2n] = (nNodeExtra++) + nNode;
569 for (auto &v : c2n)
570 v = nodeOld2Dedu.at(v);
571 auto vtkCell = Elem::ToVTKVertsAndData(elem, c2n);
572 result.cellTypes[iCell] = vtkCell.first;
573 for (auto in : vtkCell.second)
574 result.connectivity.push_back(in);
575 cellEnd += vtkCell.second.size();
576 result.offsets[iCell] = cellEnd;
577 }
578 return result;
579 }
580
581 /**
582 * \brief Write the .pvtu parallel VTK master file referencing all
583 * per-rank .vtu piece files.
584 *
585 * Only called by rank 0 (mRank) when flag==1 (parallel output).
586 */
587 void WritePVTUMasterFile(
588 const std::string &fnameIN,
589 const std::filesystem::path &outPath,
590 const std::string &endianName,
591 const std::string &seriesName,
592 double t,
593 DNDS::MPI_int nMPIRanks,
594 int arraySiz, int vecArraySiz,
595 int arraySizPoint, int vecArraySizPoint,
596 const tFGetName &names,
597 const tFGetName &vectorNames,
598 const tFGetName &namesPoint,
599 const tFGetName &vectorNamesPoint,
600 const std::string &indentV,
601 const std::string &newlineV)
602 {
603 // Replicate writeXMLEntity locally (thin wrapper for XML output)
604 auto writeXMLEntity = [&](std::ostream &out, int level, const auto &name,
605 const std::vector<std::pair<std::string, std::string>> &attr,
606 auto &&writeContent) -> void
607 {
608 for (int i = 0; i < level; i++)
609 out << indentV;
610 out << "<" << name << " ";
611 for (const auto &a : attr)
612 out << a.first << "=\"" << a.second << "\" ";
613 out << ">" << newlineV;
614 writeContent(out, level + 1);
615 for (int i = 0; i < level; i++)
616 out << indentV;
617 out << "</" << name << ">" << newlineV;
618 };
619
620 std::filesystem::path foutPP{fnameIN + ".pvtu"};
621 std::ofstream foutP{foutPP};
622 DNDS_assert(foutP);
623 if (seriesName.size())
624 updateVTKSeries(seriesName + ".pvtu.series", getStringForcePath(foutPP.filename()), t);
625
626 writeXMLEntity(
627 foutP, 0, "VTKFile",
628 {{"type", "PUnstructuredGrid"},
629 {"byte_order", endianName},
630 {"header_type", "UInt64"}},
631 [&](auto &out, int level)
632 {
633 writeXMLEntity(
634 out, level, "PUnstructuredGrid",
635 {{"GhostLevel", "0"}},
636 [&](auto &out, int level)
637 {
638 {
639 std::string namesAll;
640 for (int i = 0; i < arraySiz; i++)
641 namesAll.append(names(i) + " ");
642 std::string vectorsNameAll;
643 for (int i = 0; i < vecArraySiz; i++)
644 vectorsNameAll.append(vectorNames(i) + " ");
645 writeXMLEntity(
646 out, level, "PCellData",
647 {{"Scalars", namesAll},
648 {"Vectors", vectorsNameAll}},
649 [&](auto &out, int level)
650 {
651 for (int i = 0; i < arraySiz; i++)
652 {
653 writeXMLEntity(
654 out, level, "PDataArray",
655 {{"type", "Float64"},
656 {"Name", names(i)},
657 {"format", "ascii"}},
658 [&](auto &out, int level) {});
659 }
660 for (int i = 0; i < vecArraySiz; i++)
661 {
662 writeXMLEntity(
663 out, level, "PDataArray",
664 {{"type", "Float64"},
665 {"Name", vectorNames(i)},
666 {"NumberOfComponents", "3"},
667 {"format", "ascii"}},
668 [&](auto &out, int level) {});
669 }
670 });
671 }
672
673 {
674 std::string namesAll;
675 for (int i = 0; i < arraySizPoint; i++)
676 namesAll.append(namesPoint(i) + " ");
677 std::string vectorsNameAll;
678 for (int i = 0; i < vecArraySizPoint; i++)
679 vectorsNameAll.append(vectorNamesPoint(i) + " ");
680 writeXMLEntity(
681 out, level, "PPointData",
682 {{"Scalars", namesAll},
683 {"Vectors", vectorsNameAll}},
684 [&](auto &out, int level)
685 {
686 for (int i = 0; i < arraySizPoint; i++)
687 {
688 writeXMLEntity(
689 out, level, "PDataArray",
690 {{"type", "Float64"},
691 {"Name", namesPoint(i)},
692 {"format", "ascii"}},
693 [&](auto &out, int level) {});
694 }
695 for (int i = 0; i < vecArraySizPoint; i++)
696 {
697 writeXMLEntity(
698 out, level, "PDataArray",
699 {{"type", "Float64"},
700 {"Name", vectorNamesPoint(i)},
701 {"NumberOfComponents", "3"},
702 {"format", "ascii"}},
703 [&](auto &out, int level) {});
704 }
705 });
706 }
707 writeXMLEntity(
708 out, level, "PPoints",
709 {},
710 [&](auto &out, int level)
711 {
712 writeXMLEntity(
713 out, level, "PDataArray",
714 {{"type", "Float64"},
715 {"NumberOfComponents", "3"}},
716 [&](auto &out, int level) {});
717 });
718 for (DNDS::MPI_int iRank = 0; iRank < nMPIRanks; iRank++)
719 {
720 char BUF[512];
721 std::sprintf(BUF, "%04d", iRank);
722 std::string cFileNameRelPVTU = getStringForcePath(
723 outPath.lexically_relative(outPath.parent_path()) / (std::string(BUF) + ".vtu"));
724 writeXMLEntity(
725 out, level, "Piece",
726 {{"Source", cFileNameRelPVTU}},
727 [&](auto &out, int level) {});
728 }
729 });
730 });
731 }
732 } // anonymous namespace
733
734 /**
735 * @brief referencing https://docs.vtk.org/en/latest/design_documents/VTKFileFormats.html
736 *
737 * @param fname
738 * @param arraySiz
739 * @param vecArraySiz
740 * @param arraySizPoint
741 * @param vecArraySizPoint
742 * @param names
743 * @param data
744 * @param vectorNames
745 * @param vectorData
746 * @param namesPoint
747 * @param dataPoint
748 * @param vectorNamesPoint
749 * @param vectorDataPoint
750 * @param t
751 * @param flag
752 */
753 void UnstructuredMeshSerialRW::PrintSerialPartVTKDataArray(
754 std::string fname, std::string seriesName,
755 int arraySiz, int vecArraySiz, int arraySizPoint, int vecArraySizPoint,
756 const tFGetName &names,
757 const tFGetData &data,
758 const tFGetName &vectorNames,
759 const tFGetVecData &vectorData,
760 const tFGetName &namesPoint,
761 const tFGetData &dataPoint,
762 const tFGetName &vectorNamesPoint,
763 const tFGetVecData &vectorDataPoint,
764 double t, int flag)
765 {
766 auto mpi = mesh->getMPI();
767 std::string fnameIN = fname;
768
769 if (mpi.rank != mRank && flag == 0) //* now only operating on mRank if serial
770 return;
771
772 if (flag == 0)
773 {
774 fname += ".vtu";
775 std::filesystem::path outFile{fname};
776 std::filesystem::create_directories(outFile.parent_path() / ".");
777 DNDS_assert(mode == SerialOutput && dataIsSerialOut);
778 if (seriesName.size())
779 updateVTKSeries(seriesName + ".vtu.series", getStringForcePath(outFile.filename()), t);
780 }
781 std::filesystem::path outPath; // only valid if parallel out
782 if (flag == 1)
783 {
784 outPath = {fname + ".vtu.dir"};
785 std::filesystem::create_directories(outPath);
786 char BUF[512];
787 std::sprintf(BUF, "%04d", mpi.rank);
788 fname = getStringForcePath(outPath / (std::string(BUF) + ".vtu"));
789 }
790
791 std::ofstream fout(fname);
792 if (!fout)
793 {
794 DNDS::log() << "Error: PrintSerialPartVTKDataArray open \"" << fname << "\" failure" << std::endl;
795 DNDS_assert(false);
796 }
797
798 /*******************************************************/
799 // output preparation
800 tCoordPair coordOut;
801 tAdjPair cell2nodeOut;
802 tPbiPair cell2nodePbiOut;
803 tElemInfoArrayPair cellElemInfoOut;
804 index nCell{-1}, nNode{-1};
805 this->GetCurrentOutputArrays(flag, coordOut, cell2nodeOut, cell2nodePbiOut, cellElemInfoOut,
806 nCell, nNode);
807
808 std::vector<Geom::tPoint>
809 nodesExtra;
810 std::vector<index> nodesExtraAtOriginal;
811 if (mesh->isPeriodic)
812 {
813 for (index iCell = 0; iCell < nCell; iCell++)
814 for (rowsize ic2n = 0; ic2n < cell2nodePbiOut.RowSize(iCell); ic2n++)
815 if (cell2nodePbiOut[iCell][ic2n])
816 nodesExtra.push_back(mesh->periodicInfo.GetCoordByBits(coordOut[cell2nodeOut[iCell][ic2n]], cell2nodePbiOut[iCell][ic2n])),
817 nodesExtraAtOriginal.push_back(cell2nodeOut(iCell, ic2n));
818 }
819 // if (mpi.rank == mRank)
820 // std::cout << "PrintSerialPartVTKDataArray === " << std::endl;
821 // auto printSize = [&]()
822 // { std::cout << "Rank [" << mpi.rank << "]"
823 // << " Size: " << nNode << " " << nCell << " " << nodesExtra.size() << std::endl; };
824 // if (flag == 1)
825 // MPISerialDo(mpi, printSize);
826 // else
827 // printSize();
828 std::vector<index> nodeDedu2Old;
829 std::vector<index> nodeOld2Dedu;
830 GetViolentNodeDeduplication(nNode, nodesExtra, coordOut, nodeDedu2Old, nodeOld2Dedu);
831 /*******************************************************/
832
833 std::string indentV = " ";
834 std::string newlineV = "\n";
835
836 auto writeXMLEntity = [&](std::ostream &out, int level, const auto &name, const std::vector<std::pair<std::string, std::string>> &attr, auto &&writeContent) -> void
837 {
838 for (int i = 0; i < level; i++)
839 out << indentV;
840 out << "<" << name << " ";
841 for (const auto &a : attr)
842 out << a.first << "=\"" << a.second << "\" ";
843 out << ">" << newlineV;
844
845 writeContent(out, level + 1);
846
847 for (int i = 0; i < level; i++)
848 out << indentV;
849 out << "</" << name << ">" << newlineV;
850 };
851
852 auto writeCoords = [&](auto &out, int level)
853 {
854 writeXMLEntity(
855 out, level, "Points",
856 {},
857 [&](auto &out, int level)
858 {
859 writeXMLEntity(
860 out, level, "DataArray",
861 {{"type", "Float64"},
862 {"NumberOfComponents", "3"},
863 {"format", vtuFloatEncodeMode}},
864 [&](auto &out, int level)
865 {
866 if (vtuFloatEncodeMode == "binary")
867 {
868 /************************/
869 // base64 inline
870 uint64_t binSize = (nodeDedu2Old.size() * 3) * sizeof(double);
871 std::vector<uint8_t> dataOutBytes;
872 dataOutBytes.resize(binSize + sizeof(binSize), 0);
873 size_t top{0};
874 *(uint64_t *)(dataOutBytes.data() + top) = binSize, top += sizeof(uint64_t);
875 for (index i = 0; i < nodeDedu2Old.size(); i++)
876 {
877 index iN = nodeDedu2Old.at(i);
878 if (iN < nNode)
879 {
880 *(double *)(dataOutBytes.data() + top) = double(coordOut[iN](0)), top += sizeof(double);
881 *(double *)(dataOutBytes.data() + top) = double(coordOut[iN](1)), top += sizeof(double);
882 *(double *)(dataOutBytes.data() + top) = double(coordOut[iN](2)), top += sizeof(double);
883 }
884 else
885 {
886 *(double *)(dataOutBytes.data() + top) = double(nodesExtra.at(iN - nNode)(0)), top += sizeof(double);
887 *(double *)(dataOutBytes.data() + top) = double(nodesExtra.at(iN - nNode)(1)), top += sizeof(double);
888 *(double *)(dataOutBytes.data() + top) = double(nodesExtra.at(iN - nNode)(2)), top += sizeof(double);
889 }
890 }
891 out << cppcodec::base64_rfc4648::encode(dataOutBytes);
892 }
893 else if (vtuFloatEncodeMode == "ascii")
894 {
895 std::vector<double> coordsOutData(nodeDedu2Old.size() * 3);
896 for (index i = 0; i < nodeDedu2Old.size(); i++)
897 {
898 index iN = nodeDedu2Old.at(i);
899 if (iN < nNode)
900 {
901 coordsOutData[i * 3 + 0LL] = coordOut[iN](0);
902 coordsOutData[i * 3 + 1LL] = coordOut[iN](1);
903 coordsOutData[i * 3 + 2LL] = coordOut[iN](2);
904 }
905 else
906 {
907 coordsOutData[i * 3 + 0LL] = nodesExtra.at(iN - nNode)(0);
908 coordsOutData[i * 3 + 1LL] = nodesExtra.at(iN - nNode)(1);
909 coordsOutData[i * 3 + 2LL] = nodesExtra.at(iN - nNode)(2);
910 }
911 }
912 // out << cppcodec::base64_rfc4648::encode(
913 // (uint8_t *)coordsOutData.data(),
914 // nNode * 3 * sizeof(double));
915 for (auto v : coordsOutData)
916 out << std::setprecision(ascii_precision) << v << " ";
917 }
918 else
919 DNDS_assert(false);
920 out << newlineV;
921 });
922 });
923 };
924
925 auto vtkTopo = BuildVTKCellTopology(
926 nCell, nNode, cell2nodeOut, cell2nodePbiOut,
927 cellElemInfoOut, nodeOld2Dedu, mesh->isPeriodic);
928 auto &cell2nodeOutData = vtkTopo.connectivity;
929 auto &cell2nodeOffsetData = vtkTopo.offsets;
930 auto &cellTypeData = vtkTopo.cellTypes;
931
932 auto writeCells = [&](auto &out, int level)
933 {
934 writeXMLEntity(
935 out, level, "Cells",
936 {},
937 [&](auto &out, int level)
938 {
939 writeXMLEntity(
940 out, level, "DataArray",
941 {{"type", "Int64"},
942 {"Name", "connectivity"},
943 {"format", "ascii"}},
944 [&](auto &out, int level)
945 {
946 // out << cppcodec::base64_rfc4648::encode(
947 // (uint8_t *)cell2nodeOutData.data(),
948 // cell2nodeOutData.size() * sizeof(int64_t));
949 for (auto v : cell2nodeOutData)
950 out << v << " ";
951 out << newlineV;
952 });
953 writeXMLEntity(
954 out, level, "DataArray",
955 {{"type", "Int64"},
956 {"Name", "offsets"},
957 {"format", "ascii"}},
958 [&](auto &out, int level)
959 {
960 // out << cppcodec::base64_rfc4648::encode(
961 // (uint8_t *)cell2nodeOffsetData.data(),
962 // cell2nodeOffsetData.size() * sizeof(int64_t));
963 for (auto v : cell2nodeOffsetData)
964 out << v << " ";
965 out << newlineV;
966 });
967 writeXMLEntity(
968 out, level, "DataArray",
969 {{"type", "UInt8"},
970 {"Name", "types"},
971 {"format", "ascii"}},
972 [&](auto &out, int level)
973 {
974 // out << cppcodec::base64_rfc4648::encode(
975 // (uint8_t *)cellTypeData.data(),
976 // cellTypeData.size() * sizeof(uint8_t));
977 for (int v : cellTypeData) //! cout << uint8_t(v) is ill-posed
978 out << v << " ";
979 out << newlineV;
980 });
981 });
982 };
983
984 auto writeCellData = [&](auto &out, int level)
985 {
986 std::string namesAll;
987 for (int i = 0; i < arraySiz; i++)
988 namesAll.append(names(i) + " ");
989 std::string vectorsNameAll;
990 for (int i = 0; i < vecArraySiz; i++)
991 vectorsNameAll.append(vectorNames(i) + " ");
992
993 writeXMLEntity(
994 out, level, "CellData",
995 {{"Scalars", namesAll},
996 {"Vectors", vectorsNameAll}},
997 [&](auto &out, int level)
998 {
999 for (int i = 0; i < arraySiz; i++)
1000 {
1001 writeXMLEntity(
1002 out, level, "DataArray",
1003 {{"type", "Float64"},
1004 {"Name", names(i)},
1005 {"format", vtuFloatEncodeMode}},
1006 [&](auto &out, int level)
1007 {
1008 if (vtuFloatEncodeMode == "binary")
1009 {
1010 /************************/
1011 // base64 inline
1012 uint64_t binSize = nCell * sizeof(double);
1013 std::vector<uint8_t> dataOutBytes;
1014 dataOutBytes.resize(binSize + sizeof(binSize), 0);
1015 size_t top{0};
1016 *(uint64_t *)(dataOutBytes.data() + top) = binSize, top += sizeof(uint64_t);
1017 for (index iCell = 0; iCell < nCell; iCell++)
1018 *(double *)(dataOutBytes.data() + top) = double(data(i, iCell)), top += sizeof(double);
1019 out << cppcodec::base64_rfc4648::encode(dataOutBytes);
1020 }
1021 else if (vtuFloatEncodeMode == "ascii")
1022 {
1023 /************************/
1024 // ascii
1025 std::vector<double> dataOutC(nCell);
1026 for (index iCell = 0; iCell < nCell; iCell++)
1027 dataOutC[iCell] = data(i, iCell);
1028 // auto dataOutCompressed = zlibCompressData((uint8_t *)dataOutC.data(), dataOutC.size() * sizeof(double));
1029 // out << cppcodec::base64_rfc4648::encode((char *)dataOutC.data(), binSize);
1030 out << std::setprecision(ascii_precision);
1031 for (auto v : dataOutC)
1032 out << v << " ";
1033 }
1034 else
1035 DNDS_assert(false);
1036 out << newlineV;
1037 });
1038 }
1039 writeXMLEntity(
1040 out, level, "DataArray",
1041 {{"type", "Int32"},
1042 {"Name", "iPart"},
1043 {"format", "ascii"}},
1044 [&](auto &out, int level)
1045 {
1046 std::vector<double> dataOutC(nCell);
1047 for (index iCell = 0; iCell < nCell; iCell++)
1048 {
1049 MPI_int r{0};
1050 index v{0};
1051 if (flag == 0)
1052 {
1053 bool ret = cell2nodeSerialOutTrans.pLGlobalMapping->search(iCell, r, v);
1054 DNDS_assert_info(ret, "search failed");
1055 }
1056 else if (flag == 1)
1057 r = mesh->getMPI().rank;
1058 dataOutC[iCell] = r;
1059 }
1060 out << std::setprecision(ascii_precision);
1061 for (auto v : dataOutC)
1062 out << v << " ";
1063 out << newlineV;
1064 });
1065 for (int i = 0; i < vecArraySiz; i++)
1066 {
1067 writeXMLEntity(
1068 out, level, "DataArray",
1069 {{"type", "Float64"},
1070 {"Name", vectorNames(i)},
1071 {"NumberOfComponents", "3"},
1072 {"format", vtuFloatEncodeMode}},
1073 [&](auto &out, int level)
1074 {
1075 if (vtuFloatEncodeMode == "binary")
1076 {
1077 /************************/
1078 // base64 inline
1079 uint64_t binSize = (nCell * 3) * sizeof(double);
1080 std::vector<uint8_t> dataOutBytes;
1081 dataOutBytes.resize(binSize + sizeof(binSize), 0);
1082 size_t top{0};
1083 *(uint64_t *)(dataOutBytes.data() + top) = binSize, top += sizeof(uint64_t);
1084 for (index iCell = 0; iCell < nCell; iCell++)
1085 {
1086 *(double *)(dataOutBytes.data() + top) = vectorData(i, iCell, 0), top += sizeof(double);
1087 *(double *)(dataOutBytes.data() + top) = vectorData(i, iCell, 1), top += sizeof(double);
1088 *(double *)(dataOutBytes.data() + top) = vectorData(i, iCell, 2), top += sizeof(double);
1089 }
1090 out << cppcodec::base64_rfc4648::encode(dataOutBytes);
1091 }
1092 else if (vtuFloatEncodeMode == "ascii")
1093 {
1094 /************************/
1095 // ascii
1096 std::vector<double> dataOutC(nCell * 3);
1097 for (index iCell = 0; iCell < nCell; iCell++)
1098 {
1099 dataOutC[iCell * 3 + 0LL] = vectorData(i, iCell, 0);
1100 dataOutC[iCell * 3 + 1] = vectorData(i, iCell, 1);
1101 dataOutC[iCell * 3 + 2] = vectorData(i, iCell, 2);
1102 }
1103 // out << cppcodec::base64_rfc4648::encode(
1104 // (uint8_t *)dataOutC.data(),
1105 // dataOutC.size() * sizeof(double));
1106 for (auto v : dataOutC)
1107 out << std::setprecision(ascii_precision) << v << " ";
1108 }
1109 else
1110 DNDS_assert(false);
1111 out << newlineV;
1112 });
1113 }
1114 });
1115 };
1116
1117 auto writePointData = [&](auto &out, int level)
1118 {
1119 std::string namesAll;
1120 for (int i = 0; i < arraySizPoint; i++)
1121 namesAll.append(namesPoint(i) + " ");
1122 std::string vectorsNameAll;
1123 for (int i = 0; i < vecArraySizPoint; i++)
1124 vectorsNameAll.append(vectorNamesPoint(i) + " ");
1125
1126 writeXMLEntity(
1127 out, level, "PointData",
1128 {{"Scalars", namesAll},
1129 {"Vectors", vectorsNameAll}},
1130 [&](auto &out, int level)
1131 {
1132 for (int i = 0; i < arraySizPoint; i++)
1133 {
1134 writeXMLEntity(
1135 out, level, "DataArray",
1136 {{"type", "Float64"},
1137 {"Name", namesPoint(i)},
1138 {"format", vtuFloatEncodeMode}},
1139 [&](auto &out, int level)
1140 {
1141 if (vtuFloatEncodeMode == "binary")
1142 {
1143 /************************/
1144 // base64 inline
1145 uint64_t binSize = (nodeDedu2Old.size()) * sizeof(double);
1146 std::vector<uint8_t> dataOutBytes;
1147 dataOutBytes.resize(binSize + sizeof(binSize), 0);
1148 size_t top{0};
1149 *(uint64_t *)(dataOutBytes.data() + top) = binSize, top += sizeof(uint64_t);
1150 for (index ii = 0; ii < nodeDedu2Old.size(); ii++)
1151 {
1152 index iN = nodeDedu2Old.at(ii);
1153 if (iN < nNode)
1154 *(double *)(dataOutBytes.data() + top) = dataPoint(i, iN), top += sizeof(double);
1155 else
1156 *(double *)(dataOutBytes.data() + top) = dataPoint(i, nodesExtraAtOriginal.at(iN - nNode)), top += sizeof(double);
1157 }
1158 out << cppcodec::base64_rfc4648::encode(dataOutBytes);
1159 }
1160 else if (vtuFloatEncodeMode == "ascii")
1161 {
1162 /************************/
1163 // ascii
1164 out << std::setprecision(ascii_precision);
1165 for (index ii = 0; ii < nodeDedu2Old.size(); ii++)
1166 {
1167 index iN = nodeDedu2Old.at(ii);
1168 if (iN < nNode)
1169 out << dataPoint(i, iN) << " ";
1170 else
1171 out << dataPoint(i, nodesExtraAtOriginal.at(iN - nNode)) << " ";
1172 }
1173 }
1174 else
1175 DNDS_assert(false);
1176 out << newlineV;
1177 });
1178 }
1179 for (int i = 0; i < vecArraySizPoint; i++)
1180 {
1181 writeXMLEntity(
1182 out, level, "DataArray",
1183 {{"type", "Float64"},
1184 {"Name", vectorNamesPoint(i)},
1185 {"NumberOfComponents", "3"},
1186 {"format", vtuFloatEncodeMode}},
1187 [&](auto &out, int level)
1188 {
1189 if (vtuFloatEncodeMode == "binary")
1190 {
1191 /************************/
1192 // base64 inline
1193 uint64_t binSize = (nodeDedu2Old.size()) * 3 * sizeof(double);
1194 std::vector<uint8_t> dataOutBytes;
1195 dataOutBytes.resize(binSize + sizeof(binSize), 0);
1196 size_t top{0};
1197 *(uint64_t *)(dataOutBytes.data() + top) = binSize, top += sizeof(uint64_t);
1198 for (index ii = 0; ii < nodeDedu2Old.size(); ii++)
1199 {
1200 index iN = nodeDedu2Old.at(ii);
1201 if (iN < nNode)
1202 {
1203 *(double *)(dataOutBytes.data() + top) = vectorDataPoint(i, iN, 0), top += sizeof(double);
1204 *(double *)(dataOutBytes.data() + top) = vectorDataPoint(i, iN, 1), top += sizeof(double);
1205 *(double *)(dataOutBytes.data() + top) = vectorDataPoint(i, iN, 2), top += sizeof(double);
1206 }
1207 else
1208 {
1209 *(double *)(dataOutBytes.data() + top) = vectorDataPoint(i, nodesExtraAtOriginal.at(iN - nNode), 0), top += sizeof(double);
1210 *(double *)(dataOutBytes.data() + top) = vectorDataPoint(i, nodesExtraAtOriginal.at(iN - nNode), 1), top += sizeof(double);
1211 *(double *)(dataOutBytes.data() + top) = vectorDataPoint(i, nodesExtraAtOriginal.at(iN - nNode), 2), top += sizeof(double);
1212 }
1213 }
1214 out << cppcodec::base64_rfc4648::encode(dataOutBytes);
1215 }
1216 else if (vtuFloatEncodeMode == "ascii")
1217 {
1218 /****************************/
1219 // ascii
1220 out << std::setprecision(ascii_precision);
1221 for (index ii = 0; ii < nodeDedu2Old.size(); ii++)
1222 {
1223 index iN = nodeDedu2Old.at(ii);
1224 if (iN < nNode)
1225 {
1226 out << vectorDataPoint(i, iN, 0) << " ";
1227 out << vectorDataPoint(i, iN, 1) << " ";
1228 out << vectorDataPoint(i, iN, 2) << " ";
1229 }
1230 else
1231 {
1232 out << vectorDataPoint(i, nodesExtraAtOriginal.at(iN - nNode), 0) << " ";
1233 out << vectorDataPoint(i, nodesExtraAtOriginal.at(iN - nNode), 1) << " ";
1234 out << vectorDataPoint(i, nodesExtraAtOriginal.at(iN - nNode), 2) << " ";
1235 }
1236 }
1237 }
1238 else
1239 DNDS_assert(false);
1240 out << newlineV;
1241 });
1242 }
1243 });
1244 };
1245
1246 std::string endianName{"BigEndian"};
1247 uint32_t beTest = 1; // need c++20 for STL support
1248 uint8_t beTestS = *(uint8_t *)(&beTest);
1249 if (beTestS == 1)
1250 endianName = "LittleEndian";
1251
1252 writeXMLEntity(
1253 fout, 0, "VTKFile",
1254 {
1255 {"type", "UnstructuredGrid"},
1256 {"byte_order", endianName},
1257 {"header_type", "UInt64"},
1258 // {"compressor", "vtkZLibDataCompressor"}
1259 },
1260 [&](auto &out, int level)
1261 {
1262 writeXMLEntity(
1263 out, level, "UnstructuredGrid",
1264 {},
1265 [&](auto &out, int level)
1266 {
1267 writeXMLEntity( // https://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
1268 out, level, "FieldData",
1269 {},
1270 [&](auto &out, int level)
1271 {
1272 writeXMLEntity(
1273 out, level, "DataArray",
1274 {{"type", "Float64"},
1275 {"name", "TIME"},
1276 {"NumberOfTuples", "1"},
1277 {"format", "ascii"}},
1278 [&](auto &out, int level)
1279 {
1280 out << std::setprecision(16);
1281 out << t << '\n';
1282 });
1283 });
1284 writeXMLEntity(
1285 out, level, "Piece",
1286 {{"NumberOfPoints", std::to_string(nodeDedu2Old.size())},
1287 {"NumberOfCells", std::to_string(nCell)}},
1288 [&](auto &out, int level)
1289 {
1290 /************************/
1291 // coord data
1292 writeCoords(out, level);
1293 // std::cout << "CoordDone" << std::endl;
1294 writeCells(out, level);
1295 // std::cout << "CellDone" << std::endl;
1296 writeCellData(out, level);
1297 // std::cout << "CellDataDone" << std::endl;
1298 writePointData(out, level);
1299 // std::cout << "PointDataDone" << std::endl;
1300 });
1301 });
1302 });
1303
1304 if (mpi.rank == mRank && flag == 1)
1305 {
1306 WritePVTUMasterFile(
1307 fnameIN, outPath, endianName, seriesName, t, mpi.size,
1308 arraySiz, vecArraySiz, arraySizPoint, vecArraySizPoint,
1309 names, vectorNames, namesPoint, vectorNamesPoint,
1310 indentV, newlineV);
1311 }
1312 }
1313
1314 static void H5_WriteDataset(hid_t loc, const char *name, index nGlobal, index nOffset, index nLocal,
1315 hid_t file_dataType, hid_t mem_dataType, hid_t plist_id, hid_t dcpl_id, void *buf, int dim2 = -1)
1316 {
1317 int herr{0};
1318 DNDS_assert_info(nGlobal >= 0 && nLocal >= 0 && nOffset >= 0,
1319 fmt::format("{},{},{}", nGlobal, nLocal, nOffset));
1320 int rank = dim2 >= 0 ? 2 : 1;
1321 std::array<hsize_t, 2> ranksFull{hsize_t(nGlobal), hsize_t(dim2)};
1322 std::array<hsize_t, 2> ranksFullUnlim{H5S_UNLIMITED, hsize_t(dim2)};
1323 std::array<hsize_t, 2> offset{hsize_t(nOffset), 0};
1324 std::array<hsize_t, 2> siz{hsize_t(nLocal), hsize_t(dim2)};
1325 hid_t memSpace = H5Screate_simple(rank, siz.data(), NULL);
1326 hid_t fileSpace = H5Screate_simple(rank, ranksFull.data(), ranksFullUnlim.data());
1327 hid_t dset_id = H5Dcreate(loc, name, file_dataType, fileSpace, H5P_DEFAULT, dcpl_id, H5P_DEFAULT);
1328 DNDS_assert_info(H5I_INVALID_HID != dset_id, "dataset create failed");
1329 herr = H5Sclose(fileSpace);
1330 fileSpace = H5Dget_space(dset_id);
1331 herr |= H5Sselect_hyperslab(fileSpace, H5S_SELECT_SET, offset.data(), NULL, siz.data(), NULL);
1332 herr |= H5Dwrite(dset_id, mem_dataType, memSpace, fileSpace, plist_id, buf);
1333 herr |= H5Dclose(dset_id);
1334 herr |= H5Sclose(fileSpace);
1335 herr |= H5Sclose(memSpace);
1336 DNDS_assert_info(herr >= 0,
1337 "h5 error " + fmt::format(
1338 "nGlobal {}, nOffset {}, nLocal {}",
1339 nGlobal, nOffset, nLocal));
1340 }
1341
1342 void UnstructuredMesh::PrintParallelVTKHDFDataArray(
1343 std::string fname, std::string seriesName,
1344 int arraySiz, int vecArraySiz, int arraySizPoint, int vecArraySizPoint,
1345 const tFGetName &names,
1346 const tFGetData &data,
1347 const tFGetName &vectorNames,
1348 const tFGetVecData &vectorData,
1349 const tFGetName &namesPoint,
1350 const tFGetData &dataPoint,
1351 const tFGetName &vectorNamesPoint,
1352 const tFGetVecData &vectorDataPoint,
1353 double t, MPI_Comm commDup)
1354 {
1355 fname += ".vtkhdf";
1356 std::filesystem::path outFile{fname};
1357 // if (mpi.rank == mRank) // only mRank creates could be faulty?
1358 std::filesystem::create_directories(outFile.parent_path() / ".");
1359 if (mpi.rank == mRank)
1360 if (seriesName.size())
1361 updateVTKSeries(seriesName + ".vtkhdf.series", getStringForcePath(outFile.filename()), t);
1362
1363 if (isPeriodic)
1364 DNDS_assert(coordsPeriodicRecreated.father);
1365 else
1366 DNDS_assert(coords.father);
1367 tCoord coordOut = isPeriodic ? coordsPeriodicRecreated.father : coords.father;
1368
1369 // MPI::AllreduceOneIndex(nNodesLocal, MPI_SUM, mpi);
1370 // long long numberOfNodes = coordOut->globalSize();
1371 long long numberOfNodes = vtkNNodeGlobal;
1372 DNDS_assert(cell2node.father);
1373
1374 // long long numberOfCells = cell2node.father->globalSize();
1375 long long numberOfCells = vtkNCellGlobal;
1376 long long numberOfConnectivity = vtkCell2NodeGlobalSiz;
1377
1378 index gSize = 0;
1379 index cSize = coordOut->Size();
1380 MPI::Allreduce(&cSize, &gSize, 1, DNDS_MPI_INDEX, MPI_SUM, commDup);
1381 long long numberOfNodes1 = gSize;
1382
1383 DNDS_assert(cell2node.father);
1384 gSize = 0;
1385 cSize = cell2node.father->Size();
1386 MPI::Allreduce(&cSize, &gSize, 1, DNDS_MPI_INDEX, MPI_SUM, commDup);
1387 long long numberOfCells1 = gSize;
1388
1389 DNDS_assert(numberOfNodes1 == numberOfNodes);
1390 DNDS_assert(numberOfCells1 == numberOfCells);
1391 // std::cout << vtkNCellGlobal << " " << vtkNNodeGlobal << " ";
1392 // std::cout << numberOfCells << " " << numberOfNodes << " ";
1393 // std::cout << std::endl;
1394 // return;
1395 // std::this_thread::sleep_for(std::chrono::seconds(1));
1396
1397 herr_t herr{0};
1398#define H5SS DNDS_assert_info(herr >= 0, "H5 setting err")
1399#define H5S_Close DNDS_assert_info(herr >= 0, "H5 closing err")
1400 hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
1401
1402 herr = H5Pset_fapl_mpio(plist_id, commDup, MPI_INFO_NULL), H5SS; // Set up file access property list with parallel I/O access
1403 herr = H5Pset_all_coll_metadata_ops(plist_id, hdf5OutSetting.coll_on_meta), H5SS;
1404 herr = H5Pset_coll_metadata_write(plist_id, hdf5OutSetting.coll_on_meta), H5SS;
1405 hid_t file_id = H5Fcreate(fname.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
1406 DNDS_assert_info(H5I_INVALID_HID != file_id, "file open failed");
1407 herr = H5Pclose(plist_id), H5S_Close;
1408 DNDS_assert(herr >= 0 && file_id);
1409
1410 hid_t VTKHDF_group_id = H5Gcreate(file_id, "/VTKHDF", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1411 DNDS_assert_info(H5I_INVALID_HID != file_id, "group create failed");
1412
1413 {
1414 hid_t scalar_space = H5Screate(H5S_SCALAR);
1415 hid_t string_type = H5Tcreate(H5T_STRING, sizeof("UnstructuredGrid") - 1); // this is weird
1416 herr |= H5Tset_strpad(string_type, H5T_STR_NULLPAD);
1417 hid_t type_attr_id = H5Acreate(VTKHDF_group_id, "Type", string_type, scalar_space, H5P_DEFAULT, H5P_DEFAULT);
1418 herr |= H5Awrite(type_attr_id, string_type, "UnstructuredGrid");
1419 H5Aclose(type_attr_id);
1420 H5Tclose(string_type);
1421 H5Sclose(scalar_space);
1422 }
1423 DNDS_assert_info(herr >= 0, "h5 error");
1424
1425 // H5LTset_attribute_string(file_id, "VTKHDF", "Type", "UnstructuredGrid");
1426 std::array<long long, 2> version{1, 0};
1427 // herr |= H5LTset_attribute_long_long(file_id, "VTKHDF", "Version", version.data(), 2);
1428 {
1429 std::array<hsize_t, 2> dims{2, 1};
1430 hid_t space = H5Screate_simple(1, dims.data(), NULL);
1431 hid_t type_attr_id = H5Acreate(VTKHDF_group_id, "Version", H5T_NATIVE_LLONG, space, H5P_DEFAULT, H5P_DEFAULT);
1432 herr |= H5Awrite(type_attr_id, H5T_NATIVE_LLONG, version.data());
1433 H5Aclose(type_attr_id);
1434 H5Sclose(space);
1435 }
1436
1437 std::array<hsize_t, 1> numberSiz{1};
1438 // herr |= H5LTmake_dataset(VTKHDF_group_id, "NumberOfCells", 1, numberSiz.data(), H5T_NATIVE_LLONG, &numberOfCells);
1439 // herr |= H5LTmake_dataset(VTKHDF_group_id, "NumberOfPoints", 1, numberSiz.data(), H5T_NATIVE_LLONG, &numberOfNodes);
1440 // herr |= H5LTmake_dataset(VTKHDF_group_id, "NumberOfConnectivityIds", 1, numberSiz.data(), H5T_NATIVE_LLONG, &numberOfConnectivity);
1441
1442 DNDS_assert(herr >= 0);
1443
1444 plist_id = H5Pcreate(H5P_DATASET_XFER);
1445 if (hdf5OutSetting.coll_on_data || hdf5OutSetting.deflateLevel > 0)
1446 herr |= H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
1447 else
1448 herr |= H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT);
1449 std::array<hsize_t, 2> chunk_dims{hdf5OutSetting.chunkSize, 3};
1450 hid_t dcpl_id = H5Pcreate(H5P_DATASET_CREATE);
1451 herr |= H5Pset_chunk(dcpl_id, 1, chunk_dims.data());
1452 hid_t dcpl_id_3arr = H5Pcreate(H5P_DATASET_CREATE);
1453 herr |= H5Pset_chunk(dcpl_id_3arr, 2, chunk_dims.data());
1454 DNDS_assert_info(herr >= 0, "h5 error");
1455#ifdef H5_HAVE_FILTER_DEFLATE
1456 if (hdf5OutSetting.deflateLevel > 0)
1457 herr = H5Pset_deflate(dcpl_id, hdf5OutSetting.deflateLevel);
1458 if (hdf5OutSetting.deflateLevel > 0)
1459 herr = H5Pset_deflate(dcpl_id_3arr, hdf5OutSetting.deflateLevel);
1460#endif
1461 int scalarSize = mpi.rank == mRank ? 1 : 0;
1462 H5_WriteDataset(VTKHDF_group_id, "NumberOfCells", 1, 0, scalarSize,
1463 H5T_NATIVE_LLONG, H5T_NATIVE_LLONG, plist_id, dcpl_id,
1464 &numberOfCells);
1465 H5_WriteDataset(VTKHDF_group_id, "NumberOfPoints", 1, 0, scalarSize,
1466 H5T_NATIVE_LLONG, H5T_NATIVE_LLONG, plist_id, dcpl_id,
1467 &numberOfNodes);
1468 H5_WriteDataset(VTKHDF_group_id, "NumberOfConnectivityIds", 1, 0, scalarSize,
1469 H5T_NATIVE_LLONG, H5T_NATIVE_LLONG, plist_id, dcpl_id,
1470 &numberOfConnectivity);
1471
1472 /************************************************************/
1473 // log() << fmt::format("coordsSize {} cellSize {} ", coordOut->Size(), cell2node.father->Size());
1474 // log() << fmt::format("numberOfNodes {}, numberOfCells {} numberOfCon {}", numberOfNodes, numberOfCells, numberOfConnectivity) << std::endl;
1475 H5_WriteDataset(VTKHDF_group_id, "Points", numberOfNodes, vtkNodeOffset, coordOut->Size(),
1476 H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, plist_id, dcpl_id_3arr, coordOut->data(), 3);
1477 /************************************************************/
1478 DNDS_assert(vtkCell2nodeOffsets.size() == cell2node.father->Size() + 1);
1479 H5_WriteDataset(VTKHDF_group_id, "Offsets", numberOfCells + 1, vtkCellOffset, cell2node.father->Size() + (mpi.rank == mpi.size - 1 ? 1 : 0),
1480 H5T_NATIVE_LLONG, H5T_NATIVE_LLONG, plist_id, dcpl_id, vtkCell2nodeOffsets.data(), -1);
1481 H5_WriteDataset(VTKHDF_group_id, "Types", numberOfCells, vtkCellOffset, cell2node.father->Size(),
1482 H5T_NATIVE_UINT8, H5T_NATIVE_UINT8, plist_id, dcpl_id, vtkCellType.data(), -1);
1483 H5_WriteDataset(VTKHDF_group_id, "Connectivity", vtkCell2NodeGlobalSiz, vtkCell2nodeOffsets.front(), vtkCell2node.size(),
1484 H5T_NATIVE_LLONG, H5T_NATIVE_LLONG, plist_id, dcpl_id, vtkCell2node.data(), -1);
1485
1486 /************************************************************/
1487
1488 {
1489 std::vector<double> cellDataBuf(cell2node.father->Size() * 3);
1490 hid_t CellData_group_id = H5Gcreate(VTKHDF_group_id, "CellData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1491 DNDS_assert_info(H5I_INVALID_HID != CellData_group_id, "group create failed");
1492
1493 for (int iArr = 0; iArr < arraySiz; iArr++)
1494 {
1495 for (index iC = 0; iC < cell2node.father->Size(); iC++)
1496 cellDataBuf[iC] = data(iArr, iC);
1497
1498 auto arrName = names(iArr);
1499 H5_WriteDataset(CellData_group_id, arrName.c_str(), numberOfCells, vtkCellOffset, cell2node.father->Size(),
1500 H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, plist_id, dcpl_id, cellDataBuf.data(), -1);
1501 }
1502 for (int iArr = 0; iArr < vecArraySiz; iArr++)
1503 {
1504 for (index iC = 0; iC < cell2node.father->Size(); iC++)
1505 for (int iRow = 0; iRow < 3; iRow++)
1506 cellDataBuf[iC * 3 + iRow] = vectorData(iArr, iC, iRow);
1507
1508 auto arrName = vectorNames(iArr);
1509 H5_WriteDataset(CellData_group_id, arrName.c_str(), numberOfCells, vtkCellOffset, cell2node.father->Size(),
1510 H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, plist_id, dcpl_id_3arr, cellDataBuf.data(), 3);
1511 }
1512
1513 herr = H5Gclose(CellData_group_id), H5S_Close;
1514 }
1515 if (isPeriodic)
1516 {
1517 DNDS_assert(coordsPeriodicRecreated.father && coordsPeriodicRecreated.father->Size() == nodeRecreated2nodeLocal.size());
1518 }
1519
1520 {
1521 index iNMax = coordOut->Size();
1522 std::vector<double> nodeDataBuf(iNMax * 3);
1523 hid_t PointData_group_id = H5Gcreate(VTKHDF_group_id, "PointData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1524 DNDS_assert_info(H5I_INVALID_HID != PointData_group_id, "group create failed");
1525
1526 for (int iArr = 0; iArr < arraySizPoint; iArr++)
1527 {
1528 for (index iN = 0; iN < iNMax; iN++)
1529 nodeDataBuf[iN] = dataPoint(iArr, isPeriodic ? nodeRecreated2nodeLocal.at(iN) : iN);
1530
1531 auto arrName = namesPoint(iArr);
1532 H5_WriteDataset(PointData_group_id, arrName.c_str(), numberOfNodes, vtkNodeOffset, coordOut->Size(),
1533 H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, plist_id, dcpl_id, nodeDataBuf.data(), -1);
1534 }
1535
1536 for (int iArr = 0; iArr < vecArraySizPoint; iArr++)
1537 {
1538 for (index iN = 0; iN < iNMax; iN++)
1539 for (int iRow = 0; iRow < 3; iRow++)
1540 nodeDataBuf[iN * 3 + iRow] = vectorDataPoint(iArr, isPeriodic ? nodeRecreated2nodeLocal.at(iN) : iN, iRow);
1541
1542 auto arrName = vectorNamesPoint(iArr);
1543 H5_WriteDataset(PointData_group_id, arrName.c_str(), numberOfNodes, vtkNodeOffset, coordOut->Size(),
1544 H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, plist_id, dcpl_id_3arr, nodeDataBuf.data(), 3);
1545 }
1546
1547 herr = H5Gclose(PointData_group_id), H5S_Close;
1548 }
1549
1550 herr = H5Pclose(plist_id), H5S_Close;
1551 herr = H5Pclose(dcpl_id), H5S_Close;
1552 herr = H5Pclose(dcpl_id_3arr), H5S_Close;
1553
1554 herr = H5Gclose(VTKHDF_group_id), H5S_Close;
1555 herr = H5Fclose(file_id), H5S_Close;
1556 }
1557
1558 void UnstructuredMesh::PrintMeshCGNS(std::string fname, const t_FBCID_2_Name &fbcid2name, const std::vector<std::string> &allNames)
1559 {
1560 DNDS_assert_info(!isPeriodic, "PrintMeshCGNS does not handle periodic mesh!");
1561 /*****************************/
1562 /* Data preparation: */
1563 /*****************************/
1564 DNDS_assert(this->adjPrimaryState == Adj_PointToLocal);
1565 this->AdjLocal2GlobalPrimary();
1566
1567 const index nBnd = this->NumBnd();
1568 const index nBndGlobal = this->NumBndGlobal();
1569 const index nCell = this->NumCell();
1570 const index nCellGlobal = this->NumCellGlobal();
1571 const index nBndCell = nBnd + nCell;
1572 const index nNode = this->NumNode();
1573 const index nNodeGlobal = this->NumNodeGlobal();
1574 index nNodeProcOffset{0};
1575
1576 MPI::Scan(&nNode, &nNodeProcOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
1577 DNDS_assert(mpi.rank != mpi.size - 1 || nNodeProcOffset == nNodeGlobal);
1578 nNodeProcOffset -= nNode;
1579
1580 index nBndCellProcOffset{0};
1581 MPI::Scan(&nBndCell, &nBndCellProcOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
1582 nBndCellProcOffset -= nBndCell;
1583
1584 std::vector<cgsize_t> bndCellOffsets(nBnd + nCell + 1);
1585 bndCellOffsets[0] = 0;
1586 for (index iBnd = 0; iBnd < nBnd; iBnd++)
1587 bndCellOffsets[iBnd + 1] = bndCellOffsets[iBnd] + bnd2node[iBnd].size() + 1;
1588 for (index iCell = 0; iCell < nCell; iCell++)
1589 bndCellOffsets[iCell + nBnd + 1] = bndCellOffsets[iCell + nBnd] + cell2node[iCell].size() + 1;
1590 std::vector<cgsize_t> elementPolyData(bndCellOffsets.back());
1591 for (index iBnd = 0; iBnd < nBnd; iBnd++)
1592 {
1593 elementPolyData.at(bndCellOffsets[iBnd]) = __getCGNSTypeFromElemType(this->bndElemInfo(iBnd, 0).getElemType());
1594 for (rowsize ib2n = 0; ib2n < bnd2node[iBnd].size(); ib2n++)
1595 elementPolyData.at(bndCellOffsets[iBnd] + 1 + ib2n) = bnd2node[iBnd][ib2n] + 1; // global node index, 1 based
1596 }
1597 for (index iCell = 0; iCell < nCell; iCell++)
1598 {
1599 elementPolyData.at(bndCellOffsets[iCell + nBnd]) = __getCGNSTypeFromElemType(this->cellElemInfo(iCell, 0).getElemType());
1600 for (rowsize ic2n = 0; ic2n < cell2node[iCell].size(); ic2n++)
1601 elementPolyData.at(bndCellOffsets[iCell + nBnd] + 1 + ic2n) = cell2node[iCell][ic2n] + 1; // global node index, 1 based
1602 }
1603 index bndCellOffsetsProcOffset{0};
1604 MPI::Scan(&bndCellOffsets.back(), &bndCellOffsetsProcOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
1605 index bndCellOffsetsGlobalMax{bndCellOffsetsProcOffset};
1606 MPI::Bcast(&bndCellOffsetsGlobalMax, 1, DNDS_MPI_INDEX, mpi.size - 1, mpi.comm);
1607 bndCellOffsetsProcOffset -= bndCellOffsets.back();
1608 for (auto &v : bndCellOffsets)
1609 v += bndCellOffsetsProcOffset; // made into global offset (0-based)
1610
1611 std::vector<double> coordsX(nNode), coordsY(nNode), coordsZ(nNode);
1612 for (index iNode = 0; iNode < nNode; iNode++)
1613 coordsX[iNode] = coords[iNode](0),
1614 coordsY[iNode] = coords[iNode](1),
1615 coordsZ[iNode] = coords[iNode](2);
1616
1617 std::map<std::string, std::vector<cgsize_t>> bocoLists;
1618 for (auto &name : allNames)
1619 bocoLists.insert({name, std::vector<cgsize_t>()});
1620 for (index iBnd = 0; iBnd < nBnd; iBnd++)
1621 {
1622
1623 std::string bcName = fbcid2name(this->bndElemInfo(iBnd, 0).zone);
1624 bocoLists[bcName].push_back(iBnd + nBndCellProcOffset + 1);
1625 }
1626 std::vector<index> bocoListLengths(allNames.size());
1627 for (size_t i = 0; i < allNames.size(); i++)
1628 bocoListLengths[i] = bocoLists[allNames[i]].size();
1629 std::vector<index> bocoListLengthsGlobal{bocoListLengths};
1630 MPI::Allreduce(bocoListLengths.data(), bocoListLengthsGlobal.data(), bocoListLengths.size(), DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
1631
1632 this->AdjGlobal2LocalPrimary();
1633 /*****************************/
1634 /* CGNS operations: */
1635 /*****************************/
1636
1637 std::filesystem::path outFile{fname};
1638 std::filesystem::create_directories(outFile.parent_path() / ".");
1639
1640 int cgns_file{0};
1641 DNDS_CGNS_CALL_EXIT(cgp_open(fname.c_str(), CG_MODE_WRITE, &cgns_file));
1642
1643 int iBase{0};
1644 DNDS_CGNS_CALL_EXIT(cg_base_write(cgns_file, "Base_0", this->dim, 3, &iBase));
1645 int iZone{0};
1646 std::array<cgsize_t, 3> zone_sizes{nNodeGlobal, nCellGlobal, 0};
1647 DNDS_CGNS_CALL_EXIT(cg_zone_write(cgns_file, iBase, "Zone_0", zone_sizes.data(), Unstructured, &iZone));
1648 std::array<int, 3> iCoords;
1649 DNDS_CGNS_CALL_EXIT(cgp_coord_write(cgns_file, iBase, iZone, RealDouble, "CoordinateX", &iCoords[0]));
1650 DNDS_CGNS_CALL_EXIT(cgp_coord_write(cgns_file, iBase, iZone, RealDouble, "CoordinateY", &iCoords[1]));
1651 DNDS_CGNS_CALL_EXIT(cgp_coord_write(cgns_file, iBase, iZone, RealDouble, "CoordinateZ", &iCoords[2]));
1652 std::array<cgsize_t, 1> rmin{nNodeProcOffset + 1}; // note 1 based
1653 std::array<cgsize_t, 1> rmax{nNodeProcOffset + nNode}; // note inclusive end
1654
1655 DNDS_CGNS_CALL_EXIT(cgp_coord_write_data(cgns_file, iBase, iZone, iCoords[0], rmin.data(), rmax.data(), coordsX.data()));
1656 DNDS_CGNS_CALL_EXIT(cgp_coord_write_data(cgns_file, iBase, iZone, iCoords[1], rmin.data(), rmax.data(), coordsY.data()));
1657 DNDS_CGNS_CALL_EXIT(cgp_coord_write_data(cgns_file, iBase, iZone, iCoords[2], rmin.data(), rmax.data(), coordsZ.data()));
1658
1659 int iSec{0};
1660 DNDS_CGNS_CALL_EXIT(cgp_poly_section_write(cgns_file, iBase, iZone, "all_elements", ElementType_t::MIXED, 1, nCellGlobal + nBndGlobal,
1661 bndCellOffsetsGlobalMax, 0, &iSec));
1662
1663 DNDS_CGNS_CALL_EXIT(cgp_poly_elements_write_data(cgns_file, iBase, iZone, iSec,
1664 nBndCellProcOffset + 1, nBndCellProcOffset + nBndCell,
1665 elementPolyData.data(), bndCellOffsets.data()));
1666
1667 static_assert(CGNS_VERSION >= 4500, "Need at least CGNS 4.5.0");
1668 for (size_t i = 0; i < allNames.size(); i++)
1669 if (bocoListLengthsGlobal[i] > 0)
1670 {
1671 int iBC;
1672 DNDS_CGNS_CALL_EXIT(cg_boco_write(cgns_file, iBase, iZone, allNames[i].data(), BCType_t::BCTypeNull, PointList,
1673 bocoListLengthsGlobal[i], NULL, &iBC)); // use NULL to postpone the data write
1674 DNDS_CGNS_CALL_EXIT(cg_boco_gridlocation_write(cgns_file, iBase, iZone, iBC, this->dim == 3 ? GridLocation_t::FaceCenter : GridLocation_t::EdgeCenter));
1675 DNDS_CGNS_CALL_EXIT(cg_goto(cgns_file, iBase, "Zone_t", iZone, "ZoneBC_t", 1, "BC_t", iBC, "PointList", 0, ""));
1676 index nBCLocal = bocoListLengths[i];
1677 index nBCLocalProcOffSet{0};
1678 MPI::Scan(&nBCLocal, &nBCLocalProcOffSet, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
1679 nBCLocalProcOffSet -= nBCLocal;
1680 DNDS_CGNS_CALL_EXIT(cgp_ptlist_write_data(cgns_file, nBCLocalProcOffSet + 1, nBCLocalProcOffSet + nBCLocal, bocoLists.at(allNames[i]).data()));
1681 }
1682
1683 if (cgp_close(cgns_file))
1684 cg_error_exit();
1685 }
1686}
DNDS_CGNS_CALL_EXIT
#define DNDS_CGNS_CALL_EXIT(call)
Definition CGNS.hpp:9
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
H5SS
#define H5SS
connectivity
std::vector< int64_t > connectivity
Flat node-index list for all cells.
Definition Mesh_Plts.cpp:525
offsets
std::vector< int64_t > offsets
Per-cell cumulative offset into connectivity.
Definition Mesh_Plts.cpp:526
cellTypes
std::vector< uint8_t > cellTypes
VTK cell type per cell.
Definition Mesh_Plts.cpp:527
H5S_Close
#define H5S_Close
DNDS::Geom::Elem::ToVTKVertsAndData
std::pair< int, std::vector< index > > ToVTKVertsAndData(Element e, const TIn &vin)
Definition Elements.hpp:566
DNDS::Geom::tFGetName
std::function< std::string(int)> tFGetName
Definition Mesh_DeviceView.hpp:85
DNDS::Geom::tPbiPair
ArrayPair< ArrayNodePeriodicBits< DNDS::NonUniformSize > > tPbiPair
Definition Mesh_DeviceView.hpp:66
DNDS::Geom::SerialOutput
@ SerialOutput
Definition Mesh.hpp:822
DNDS::Geom::tAdj
decltype(tAdjPair::father) tAdj
Definition Mesh_DeviceView.hpp:65
DNDS::Geom::tFGetData
std::function< DNDS::real(int, DNDS::index)> tFGetData
Definition Mesh_DeviceView.hpp:86
DNDS::Geom::tFGetVecData
std::function< DNDS::real(int, DNDS::index, DNDS::rowsize)> tFGetVecData
Definition Mesh_DeviceView.hpp:87
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::Geom::tElemInfoArrayPair
DNDS::ArrayPair< DNDS::ParArray< ElemInfo > > tElemInfoArrayPair
Definition Mesh_DeviceView.hpp:80
DNDS::Geom::tElemInfoArray
DNDS::ssp< DNDS::ParArray< ElemInfo > > tElemInfoArray
Definition Mesh_DeviceView.hpp:81
DNDS::Geom::t_FBCID_2_Name
std::function< std::string(t_index)> t_FBCID_2_Name
Definition BoundaryCondition.hpp:47
DNDS::Geom::__getCGNSTypeFromElemType
constexpr ElementType_t __getCGNSTypeFromElemType(Elem::ElemType et)
Definition CGNS.hpp:48
DNDS::Geom::tPbi
decltype(tPbiPair::father) tPbi
Definition Mesh_DeviceView.hpp:67
DNDS::Geom::tAdjPair
DNDS::ArrayAdjacencyPair< DNDS::NonUniformSize > tAdjPair
Definition Mesh_DeviceView.hpp:64
DNDS::Geom::tCoord
decltype(tCoordPair::father) tCoord
Definition Mesh_DeviceView.hpp:79
DNDS::MPI::Bcast
MPI_int Bcast(void *buf, MPI_int num, MPI_Datatype type, MPI_int source_rank, MPI_Comm comm)
dumb wrapper
Definition MPI.cpp:150
DNDS::MPI::Scan
MPI_int Scan(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Scan (inclusive prefix reduction).
Definition MPI.cpp:220
DNDS::MPI::Allreduce
MPI_int Allreduce(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Allreduce.
Definition MPI.cpp:203
DNDS::DNDS_MPI_INDEX
const MPI_Datatype DNDS_MPI_INDEX
MPI datatype matching index (= MPI_INT64_T).
Definition MPI.hpp:90
DNDS::rowsize
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:109
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
DNDS::getStringForcePath
std::string getStringForcePath(const std::filesystem::path::string_type &v)
Portable conversion of a platform-native path string to std::string.
Definition Defines.hpp:630
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:49
DNDS::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:190
DNDS::MPI_int
int MPI_int
MPI counterpart type for MPI_int (= C int). Used for counts and ranks in MPI calls.
Definition MPI.hpp:43
_generate_explicit_inst.fname
str fname
Definition _generate_explicit_inst.py:52
test_basic_fv.mpi
mpi
Definition test_basic_fv.py:487
DNDS::ArrayPair::father
ssp< TArray > father
Owned-side array (must be resized before ghost setup).
Definition ArrayPair.hpp:170
DNDS::ArrayPair::Size
index Size() const
Combined row count (father->Size() + son->Size()).
Definition ArrayPair.hpp:292
DNDS::ArrayPair::son
ssp< TArray > son
Ghost-side array (sized automatically by createMPITypes / BorrowAndPull).
Definition ArrayPair.hpp:172
DNDS::ArrayPair::RowSize
auto RowSize() const
Uniform row width (delegates to father).
Definition ArrayPair.hpp:258
DNDS::Geom::ElemInfo::CommType
static MPI_Datatype CommType()
Definition Mesh_DeviceView.hpp:43
DNDS::Geom::NodePeriodicBits::CommType
static MPI_Datatype CommType()
Definition PeriodicInfo.hpp:41
DNDS::Geom::UnstructuredMeshSerialRW::GetCurrentOutputArrays
void GetCurrentOutputArrays(int flag, tCoordPair &coordOut, tAdjPair &cell2nodeOut, tPbiPair &cell2nodePbiOut, tElemInfoArrayPair &cellElemInfoOut, index &nCell, index &nNode)
Definition Mesh_Plts.cpp:20
DNDS::Geom::UnstructuredMeshSerialRW::cell2nodeSerialOutTrans
DNDS::ArrayTransformerType< tAdj::element_type >::Type cell2nodeSerialOutTrans
Definition Mesh.hpp:908
DNDS::Geom::UnstructuredMeshSerialRW::mesh
DNDS::ssp< UnstructuredMesh > mesh
Definition Mesh.hpp:856
DNDS::Geom::UnstructuredMeshSerialRW::PrintSerialPartVTKDataArray
void PrintSerialPartVTKDataArray(std::string fname, std::string seriesName, int arraySiz, int vecArraySiz, int arraySizPoint, int vecArraySizPoint, const tFGetName &names, const tFGetData &data, const tFGetName &vectorNames, const tFGetVecData &vectorData, const tFGetName &namesPoint, const tFGetData &dataPoint, const tFGetName &vectorNamesPoint, const tFGetVecData &vectorDataPoint, double t, int flag=0)
names(idata) data(idata, ivolume)
Definition Mesh_Plts.cpp:753
DNDS::Geom::UnstructuredMeshSerialRW::PrintSerialPartPltBinaryDataArray
void PrintSerialPartPltBinaryDataArray(std::string fname, int arraySiz, int arraySizPoint, const tFGetName &names, const tFGetData &data, const tFGetName &namesPoint, const tFGetData &dataPoint, double t, int flag)
names(idata) data(idata, ivolume) https://tecplot.azureedge.net/products/360/current/360_data_format_...
Definition Mesh_Plts.cpp:127
DNDS::Geom::UnstructuredMesh::t_deviceView
UnstructuredMeshDeviceView< B > t_deviceView
Definition Mesh.hpp:777
DNDS::Geom::UnstructuredMesh::coordsPeriodicRecreated
tCoordPair coordsPeriodicRecreated
Definition Mesh.hpp:141
DNDS::Geom::UnstructuredMesh::vtkCellType
std::vector< uint8_t > vtkCellType
Definition Mesh.hpp:133
DNDS::Geom::UnstructuredMesh::vtkCell2nodeOffsets
std::vector< index > vtkCell2nodeOffsets
for parallel out
Definition Mesh.hpp:132
DNDS::Geom::UnstructuredMesh::coords
tCoordPair coords
reader
Definition Mesh.hpp:54
DNDS::Geom::UnstructuredMesh::hdf5OutSetting
struct DNDS::Geom::UnstructuredMesh::HDF5OutSetting hdf5OutSetting
DNDS::Geom::UnstructuredMesh::PrintParallelVTKHDFDataArray
void PrintParallelVTKHDFDataArray(std::string fname, std::string seriesName, int arraySiz, int vecArraySiz, int arraySizPoint, int vecArraySizPoint, const tFGetName &names, const tFGetData &data, const tFGetName &vectorNames, const tFGetVecData &vectorData, const tFGetName &namesPoint, const tFGetData &dataPoint, const tFGetName &vectorNamesPoint, const tFGetVecData &vectorDataPoint, double t, MPI_Comm commDup)
Definition Mesh_Plts.cpp:1342
DNDS::Geom::UnstructuredMesh::adjPrimaryState
MeshAdjState adjPrimaryState
Definition Mesh.hpp:41
DNDS::Geom::UnstructuredMesh::cellElemInfo
tElemInfoArrayPair cellElemInfo
Definition Mesh.hpp:59
DNDS::Geom::UnstructuredMesh::PrintMeshCGNS
void PrintMeshCGNS(std::string fname, const t_FBCID_2_Name &fbcid2name, const std::vector< std::string > &allNames)
Definition Mesh_Plts.cpp:1558
DNDS::Geom::UnstructuredMesh::bndElemInfo
tElemInfoArrayPair bndElemInfo
Definition Mesh.hpp:60
DNDS::Geom::UnstructuredMesh::vtkCell2node
std::vector< index > vtkCell2node
Definition Mesh.hpp:134
DNDS::Geom::UnstructuredMesh::nodeRecreated2nodeLocal
std::vector< index > nodeRecreated2nodeLocal
Definition Mesh.hpp:142
DNDS::MPIInfo::size
int size
Number of ranks in comm (-1 until initialised).
Definition MPI.hpp:221
DNDS::MPIInfo::rank
int rank
This rank's 0-based index within comm (-1 until initialised).
Definition MPI.hpp:219
DNDS::MPIInfo::comm
MPI_Comm comm
The underlying MPI communicator handle.
Definition MPI.hpp:217
DNDS::ObjName
Tag type for naming objects created via make_ssp.
Definition Defines.hpp:249
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
nLocal
constexpr DNDS::index nLocal
Definition test_ArrayTransformer.cpp:151
r
tVec r(NCells)
result
auto result
Definition test_Limiters.cpp:110