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