DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
test_MeshPipeline.cpp
Go to the documentation of this file.
1/**
2 * @file test_MeshPipeline.cpp
3 * @brief Parameterized doctest tests for the full mesh pipeline.
4 *
5 * Each mesh configuration (file, dimension, periodic settings, expected counts)
6 * is tested through the complete pipeline:
7 * - Pipeline state tracking: every state variable checked
8 * - InterpolateFace: face topology, face2cell, bnd2face
9 * - ReorderLocalCells: cell permutation, partition starts
10 * - N2CB, C2F, Facial, C2CFace adjacency validation
11 * - ConstructBndMesh, BuildVTKConnectivity
12 * - Global count conservation
13 * - All 5 Adj round-trip inverses (on one mesh config)
14 *
15 * Mesh configurations:
16 * [0] UniformSquare_10 -- 2D, 100 quad cells, non-periodic
17 * [1] IV10_10 -- 2D, 100 quad cells, periodic (isentropic vortex)
18 * [2] NACA0012_H2 -- 2D, 20816 unstructured quad cells, non-periodic (airfoil)
19 * [3] IV10U_10 -- 2D, 322 unstructured tri cells, periodic (isentropic vortex)
20 */
21
22#define DOCTEST_CONFIG_IMPLEMENT
23#include "doctest.h"
24
25#include "Geom/Mesh/Mesh.hpp"
27#include <string>
28#include <vector>
29#include <array>
30#include <cmath>
31#include <set>
32#include <fmt/core.h>
33
34using namespace DNDS;
35using namespace DNDS::Geom;
36
37// NOTE: DNDS::index, DNDS::real, DNDS::rowsize clash with POSIX symbols.
38
39// ---------------------------------------------------------------------------
40// Mesh configuration
41// ---------------------------------------------------------------------------
42struct MeshConfig
43{
44 const char *name; ///< Human-readable label
45 const char *file; ///< Filename relative to data/mesh/
46 int dim; ///< Spatial dimension
47 bool periodic; ///< Has periodic boundaries
48 tPoint translation1; ///< Periodic translation axis 1
49 tPoint translation2; ///< Periodic translation axis 2
50 tPoint translation3; ///< Periodic translation axis 3
51 DNDS::index expectedCells; ///< Expected global cell count (-1 = skip check)
52 DNDS::index expectedBnds; ///< Expected global bnd count (-1 = skip check)
53};
54
55static const MeshConfig g_configs[] = {
56 {"UniformSquare_10", "UniformSquare_10.cgns", 2, false, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, 100, 40},
57 {"IV10_10", "IV10_10.cgns", 2, true, {10, 0, 0}, {0, 10, 0}, {0, 0, 10}, 100, -1},
58 {"NACA0012_H2", "NACA0012_H2.cgns", 2, false, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, 20816, 484},
59 {"IV10U_10", "IV10U_10.cgns", 2, true, {10, 0, 0}, {0, 10, 0}, {0, 0, 10}, 322, -1},
60 {"Ball2", "Ball2.cgns", 3, false, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, 958994, 16402}, // 3D mixed-element mesh (Ball2 is np=8 only)
61};
62static constexpr int N_CONFIGS = sizeof(g_configs) / sizeof(g_configs[0]);
63
64// ---------------------------------------------------------------------------
65// Globals
66// ---------------------------------------------------------------------------
67static MPIInfo g_mpi;
68
69/// One full-pipeline mesh per config.
70static ssp<UnstructuredMesh> g_full[N_CONFIGS];
71
72// ---------------------------------------------------------------------------
73// Helpers
74// ---------------------------------------------------------------------------
75
76static std::string meshPath(const std::string &name)
77{
78 std::string f(__FILE__);
79 for (int i = 0; i < 4; i++)
80 {
81 auto pos = f.rfind('/');
82 if (pos == std::string::npos)
83 pos = f.rfind('\\');
84 if (pos != std::string::npos)
85 f = f.substr(0, pos);
86 }
87 return f + "/data/mesh/" + name;
88}
89
90static ssp<UnstructuredMesh> buildFullMesh(
91 int meshId, const MeshConfig &cfg, int nReorderParts = 1)
92{
93 if (g_mpi.rank == 0)
94 std::cout << "[buildMesh " << meshId << " " << cfg.name << "] START" << std::endl;
95
96 auto mesh = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
98
99 if (cfg.periodic)
100 {
101 tPoint zero{0, 0, 0};
102 mesh->SetPeriodicGeometry(
103 cfg.translation1, zero, zero,
104 cfg.translation2, zero, zero,
105 cfg.translation3, zero, zero);
106 }
107
108 reader.ReadFromCGNSSerial(meshPath(cfg.file));
109 reader.Deduplicate1to1Periodic(1e-8);
110 reader.BuildCell2Cell();
111
113 pOpt.metisType = "KWAY";
114 pOpt.metisUfactor = 30;
115 pOpt.metisSeed = 42;
116 pOpt.metisNcuts = 1;
117 reader.MeshPartitionCell2Cell(pOpt);
118 reader.PartitionReorderToMeshCell2Cell();
119
120 mesh->RecoverNode2CellAndNode2Bnd();
121 mesh->RecoverCell2CellAndBnd2Cell();
122 mesh->BuildGhostPrimary();
123
124 mesh->AdjGlobal2LocalPrimary();
125 mesh->AdjGlobal2LocalN2CB();
126
127 if (nReorderParts > 1)
128 mesh->ReorderLocalCells(nReorderParts);
129
130 mesh->InterpolateFace();
131 mesh->AssertOnFaces();
132
133 mesh->AdjLocal2GlobalN2CB();
134 mesh->BuildGhostN2CB();
135 mesh->AdjGlobal2LocalN2CB();
136 mesh->AssertOnN2CB();
137
138 mesh->BuildCell2CellFace();
139 mesh->AdjGlobal2LocalC2CFace();
140
141 if (cfg.periodic)
142 mesh->RecreatePeriodicNodes();
143 mesh->BuildVTKConnectivity();
144
145 if (g_mpi.rank == 0)
146 std::cout << "[buildMesh " << meshId << " " << cfg.name << "] DONE" << std::endl;
147
148 return mesh;
149}
150
151static std::vector<std::vector<DNDS::index>> snapshotAdj(
152 const tAdjPair &adj, DNDS::index nRows)
153{
154 std::vector<std::vector<DNDS::index>> out(nRows);
155 for (DNDS::index i = 0; i < nRows; i++)
156 {
157 out[i].resize(adj.RowSize(i));
158 for (DNDS::rowsize j = 0; j < adj.RowSize(i); j++)
159 out[i][j] = adj(i, j);
160 }
161 return out;
162}
163
164static std::vector<std::pair<DNDS::index, DNDS::index>> snapshotAdj2(
165 const tAdj2Pair &adj, DNDS::index nRows)
166{
167 std::vector<std::pair<DNDS::index, DNDS::index>> out(nRows);
168 for (DNDS::index i = 0; i < nRows; i++)
169 out[i] = {adj(i, 0), adj(i, 1)};
170 return out;
171}
172
173static std::vector<DNDS::index> snapshotAdj1(
174 const tAdj1Pair &adj, DNDS::index nRows)
175{
176 std::vector<DNDS::index> out(nRows);
177 for (DNDS::index i = 0; i < nRows; i++)
178 out[i] = adj(i, 0);
179 return out;
180}
181
182static void checkAdjMatchesSnapshot(
183 const tAdjPair &adj, DNDS::index nRows,
184 const std::vector<std::vector<DNDS::index>> &snap)
185{
186 for (DNDS::index i = 0; i < nRows; i++)
187 {
188 REQUIRE(adj.RowSize(i) == static_cast<DNDS::rowsize>(snap[i].size()));
189 for (DNDS::rowsize j = 0; j < adj.RowSize(i); j++)
190 CHECK(adj(i, j) == snap[i][j]);
191 }
192}
193
194// ---------------------------------------------------------------------------
195int main(int argc, char **argv)
196{
197 MPI_Init(&argc, &argv);
198 g_mpi.setWorld();
199
200 // Build one full-pipeline mesh per config.
201 for (int ic = 0; ic < N_CONFIGS; ic++)
202 {
203 // Ball2 (index 4) is only tested with np=8 — skip the expensive
204 // build at other rank counts to avoid wasting ~60s on serial
205 // read+partition of a 958K-cell 3D mesh that no test will use.
206 if (ic == 4 && g_mpi.size != 8)
207 continue;
208 g_full[ic] = buildFullMesh(ic, g_configs[ic]);
209 }
210
211 doctest::Context ctx;
212 ctx.applyCommandLine(argc, argv);
213 int res = ctx.run();
214
215 for (auto &m : g_full)
216 m.reset();
217 MPI_Finalize();
218 return res;
219}
220
221// ===========================================================================
222// Parameterized tests over all mesh configs (read-only)
223// ===========================================================================
224
225#define FOR_EACH_MESH_CONFIG(body) \
226 for (int _ci = 0; _ci < N_CONFIGS; _ci++) \
227 { \
228 /* Ball2 (index 4) is only tested with np=8 */ \
229 if (_ci == 4 && g_mpi.size != 8) \
230 continue; \
231 CAPTURE(_ci); \
232 const auto &cfg = g_configs[_ci]; \
233 auto m = g_full[_ci]; \
234 SUBCASE(cfg.name) { body } \
235 }
236
237// --- Pipeline state ---
238
239TEST_CASE("Pipeline: all five states are Local"){
241 CHECK(m->adjPrimaryState == Adj_PointToLocal);
242 CHECK(m->adjFacialState == Adj_PointToLocal);
243 CHECK(m->adjC2FState == Adj_PointToLocal);
244 CHECK(m->adjN2CBState == Adj_PointToLocal);
245 CHECK(m->adjC2CFaceState == Adj_PointToLocal);
246 })}
247
248// --- Global count conservation ---
249
250TEST_CASE("Pipeline: global cell count matches expected"){
252 DNDS::index localCells = m->NumCell();
253 DNDS::index globalCells = 0;
254 MPI_Allreduce(&localCells, &globalCells, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
255 CHECK(globalCells > 0);
256 if (cfg.expectedCells > 0)
257 CHECK(globalCells == cfg.expectedCells);
258 })}
259
260TEST_CASE("Pipeline: global bnd count matches expected"){
262 DNDS::index localBnds = m->NumBnd();
263 DNDS::index globalBnds = 0;
264 MPI_Allreduce(&localBnds, &globalBnds, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
265 if (cfg.expectedBnds > 0)
266 CHECK(globalBnds == cfg.expectedBnds);
267 else if (!cfg.periodic)
268 CHECK(globalBnds > 0); // non-periodic meshes must have boundaries
269 })}
270
271// --- InterpolateFace topology ---
272
273TEST_CASE("InterpolateFace: face count is positive"){
275 CHECK(m->NumFace() > 0);
276 })}
277
278TEST_CASE("InterpolateFace: face2cell has exactly 2 entries per face"){
280 for (DNDS::index iF = 0; iF < m->NumFace(); iF++)
281 CHECK(m->face2cell.RowSize(iF) == 2);
282 })}
283
284TEST_CASE("InterpolateFace: face2cell owner is a valid local cell"){
286 DNDS::index totalCells = m->NumCell() + m->NumCellGhost();
287 for (DNDS::index iF = 0; iF < m->NumFace(); iF++)
288 {
289 DNDS::index owner = m->face2cell(iF, 0);
290 CHECK(owner >= 0);
291 CHECK(owner < totalCells);
292 }
293 })}
294
295TEST_CASE("InterpolateFace: face2node indices in valid range"){
297 DNDS::index totalNodes = m->NumNode() + m->NumNodeGhost();
298 for (DNDS::index iF = 0; iF < m->NumFace(); iF++)
299 for (DNDS::rowsize j = 0; j < m->face2node.RowSize(iF); j++)
300 {
301 DNDS::index iN = m->face2node(iF, j);
302 CHECK(iN >= 0);
303 CHECK(iN < totalNodes);
304 }
305 })}
306
307TEST_CASE("InterpolateFace: cell2face row sizes match expected face count"){
309 for (DNDS::index iC = 0; iC < m->NumCell(); iC++)
310 {
311 auto elem = m->GetCellElement(iC);
312 int nFaceExpected = elem.GetNumFaces();
313 CHECK(m->cell2face.RowSize(iC) == nFaceExpected);
314 }
315 })}
316
317TEST_CASE("InterpolateFace: cell2face entries are valid face indices"){
319 DNDS::index totalFaces = m->NumFace() + m->NumFaceGhost();
320 for (DNDS::index iC = 0; iC < m->NumCell(); iC++)
321 for (DNDS::rowsize j = 0; j < m->cell2face.RowSize(iC); j++)
322 {
323 DNDS::index iF = m->cell2face(iC, j);
324 CHECK(iF >= 0);
325 CHECK(iF < totalFaces);
326 }
327 })}
328
329TEST_CASE("InterpolateFace: bnd2face maps to valid faces"){
331 for (DNDS::index ib = 0; ib < m->NumBnd(); ib++)
332 {
333 DNDS::index iFace = m->bnd2face(ib, 0);
334 // Periodic boundaries that became internal faces may have bnd2face == -1
335 if (iFace >= 0)
336 CHECK(iFace < m->NumFace() + m->NumFaceGhost());
337 }
338 })}
339
340TEST_CASE("InterpolateFace: face element types are valid")
341{
343 for (DNDS::index iF = 0; iF < m->NumFace(); iF++)
344 CHECK(m->GetFaceElement(iF).type != Elem::ElemType::UnknownElem);
345 })
346}
347
348// --- InterpolateFace: DSL vs Legacy comparison ---
349
350/// Build a mesh through the pipeline up to (and including) InterpolateFace,
351/// using the specified face interpolation method.
352static ssp<UnstructuredMesh> buildMeshUpToFace(
353 const MeshConfig &cfg, bool useLegacy)
354{
355 auto mesh = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
357
358 if (cfg.periodic)
359 {
360 tPoint zero{0, 0, 0};
361 mesh->SetPeriodicGeometry(
362 cfg.translation1, zero, zero,
363 cfg.translation2, zero, zero,
364 cfg.translation3, zero, zero);
365 }
366
367 reader.ReadFromCGNSSerial(meshPath(cfg.file));
368 reader.Deduplicate1to1Periodic(1e-8);
369 reader.BuildCell2Cell();
370
372 pOpt.metisType = "KWAY";
373 pOpt.metisUfactor = 30;
374 pOpt.metisSeed = 42;
375 pOpt.metisNcuts = 1;
376 reader.MeshPartitionCell2Cell(pOpt);
377 reader.PartitionReorderToMeshCell2Cell();
378
379 mesh->RecoverNode2CellAndNode2Bnd();
380 mesh->RecoverCell2CellAndBnd2Cell();
381 mesh->BuildGhostPrimary();
382 mesh->AdjGlobal2LocalPrimary();
383 mesh->AdjGlobal2LocalN2CB();
384
385 if (useLegacy)
386 mesh->InterpolateFaceLegacy();
387 else
388 mesh->InterpolateFace();
389
390 return mesh;
391}
392
393/// Collect sorted face vertex sets for all local cells.
394static std::vector<std::set<std::set<DNDS::index>>>
395collectCellFaceVertexSets(const ssp<UnstructuredMesh> &m)
396{
397 DNDS::index nCells = m->cell2node.father->Size();
398 std::vector<std::set<std::set<DNDS::index>>> result(nCells);
399 for (DNDS::index iCell = 0; iCell < nCells; iCell++)
400 {
401 auto eCell = Elem::Element{m->cellElemInfo[iCell]->getElemType()};
402 for (DNDS::rowsize j = 0; j < m->cell2face.RowSize(iCell); j++)
403 {
404 DNDS::index iFace = m->cell2face(iCell, j);
405 auto eFace = eCell.ObtainFace(j);
406 int nVerts = eFace.GetNumVertices();
407 std::set<DNDS::index> vs;
408 for (int k = 0; k < nVerts; k++)
409 vs.insert(m->face2node[iFace][k]);
410 result[iCell].insert(vs);
411 }
412 }
413 return result;
414}
415
416TEST_CASE("InterpolateFace: DSL matches Legacy on all mesh configs")
417{
418 for (int ci = 0; ci < N_CONFIGS; ci++)
419 {
420 if (ci == 4 && g_mpi.size != 8)
421 continue;
422 CAPTURE(ci);
423 const auto &cfg = g_configs[ci];
424 SUBCASE(cfg.name)
425 {
426 auto mDSL = buildMeshUpToFace(cfg, false);
427 auto mLeg = buildMeshUpToFace(cfg, true);
428
429 // Both should have the same number of local faces
430 CHECK(mDSL->NumFace() == mLeg->NumFace());
431 CHECK(mDSL->NumFaceGhost() == mLeg->NumFaceGhost());
432
433 // For each local cell, the face vertex sets must match
434 auto dslSets = collectCellFaceVertexSets(mDSL);
435 auto legSets = collectCellFaceVertexSets(mLeg);
436 DNDS::index nCells = mDSL->cell2node.father->Size();
437 CHECK(nCells == static_cast<DNDS::index>(legSets.size()));
438
439 for (DNDS::index iCell = 0; iCell < nCells; iCell++)
440 {
441 CAPTURE(iCell);
442 CHECK(dslSets[iCell] == legSets[iCell]);
443 }
444
445 // face2cell: for each local face, the set of (face vertex set, cell pair)
446 // should match. We compare by face vertex set → cell global indices.
447 // (Face ordering may differ, so we compare via vertex sets.)
448 std::map<std::set<DNDS::index>, std::pair<DNDS::index, DNDS::index>> dslF2C, legF2C;
449 for (DNDS::index iF = 0; iF < mDSL->NumFace(); iF++)
450 {
451 auto eFace = Elem::Element{mDSL->faceElemInfo(iF, 0).getElemType()};
452 std::set<DNDS::index> vs;
453 for (int k = 0; k < eFace.GetNumVertices(); k++)
454 vs.insert(mDSL->face2node[iF][k]);
455 dslF2C[vs] = {mDSL->face2cell(iF, 0), mDSL->face2cell(iF, 1)};
456 }
457 for (DNDS::index iF = 0; iF < mLeg->NumFace(); iF++)
458 {
459 auto eFace = Elem::Element{mLeg->faceElemInfo(iF, 0).getElemType()};
460 std::set<DNDS::index> vs;
461 for (int k = 0; k < eFace.GetNumVertices(); k++)
462 vs.insert(mLeg->face2node[iF][k]);
463 legF2C[vs] = {mLeg->face2cell(iF, 0), mLeg->face2cell(iF, 1)};
464 }
465 CHECK(dslF2C.size() == legF2C.size());
466
467 for (auto &[vs, dslPair] : dslF2C)
468 {
469 auto it = legF2C.find(vs);
470 REQUIRE(it != legF2C.end());
471 // Cell pairs should be same set (order of L/R might differ)
472 std::set<DNDS::index> dslCells{dslPair.first, dslPair.second};
473 std::set<DNDS::index> legCells{it->second.first, it->second.second};
474 CHECK(dslCells == legCells);
475 }
476
477 // Boundary zone assignment: same face vertex set should have same zone
478 std::map<std::set<DNDS::index>, t_index> dslZones, legZones;
479 for (DNDS::index iF = 0; iF < mDSL->NumFace(); iF++)
480 {
481 auto eFace = Elem::Element{mDSL->faceElemInfo(iF, 0).getElemType()};
482 std::set<DNDS::index> vs;
483 for (int k = 0; k < eFace.GetNumVertices(); k++)
484 vs.insert(mDSL->face2node[iF][k]);
485 dslZones[vs] = mDSL->faceElemInfo(iF, 0).zone;
486 }
487 for (DNDS::index iF = 0; iF < mLeg->NumFace(); iF++)
488 {
489 auto eFace = Elem::Element{mLeg->faceElemInfo(iF, 0).getElemType()};
490 std::set<DNDS::index> vs;
491 for (int k = 0; k < eFace.GetNumVertices(); k++)
492 vs.insert(mLeg->face2node[iF][k]);
493 legZones[vs] = mLeg->faceElemInfo(iF, 0).zone;
494 }
495 for (auto &[vs, dslZone] : dslZones)
496 {
497 auto it = legZones.find(vs);
498 REQUIRE(it != legZones.end());
499 CHECK(dslZone == it->second);
500 }
501 }
502 }
503}
504
505// --- Periodic face invariants ---
506
507using idx_pbi_t = std::pair<DNDS::index, uint8_t>;
508
509TEST_CASE("Periodic: face pbi has uniform XOR between owner and neighbor (collaborating check)")
510{
511 for (int _ci = 0; _ci < N_CONFIGS; _ci++)
512 {
513 if (_ci == 4 && g_mpi.size != 8)
514 continue;
515 const auto &cfg = g_configs[_ci];
516 if (!cfg.periodic)
517 continue;
518 auto m = g_full[_ci];
519 CAPTURE(_ci);
520 SUBCASE(cfg.name)
521 {
522 for (DNDS::index iFace = 0; iFace < m->NumFace(); iFace++)
523 {
524 DNDS::index iCellL = m->face2cell(iFace, 0);
525 DNDS::index iCellR = m->face2cell(iFace, 1);
526 if (iCellR == DNDS::UnInitIndex)
527 continue;
528
529 auto eCellL = Elem::Element{m->cellElemInfo[iCellL]->getElemType()};
530 auto eCellR = Elem::Element{m->cellElemInfo[iCellR]->getElemType()};
531
532 int slotL = -1, slotR = -1;
533 for (DNDS::rowsize j = 0; j < m->cell2face.RowSize(iCellL); j++)
534 if (m->cell2face(iCellL, j) == iFace)
535 {
536 slotL = j;
537 break;
538 }
539 for (DNDS::rowsize j = 0; j < m->cell2face.RowSize(iCellR); j++)
540 if (m->cell2face(iCellR, j) == iFace)
541 {
542 slotR = j;
543 break;
544 }
545 REQUIRE(slotL >= 0);
546 REQUIRE(slotR >= 0);
547
548 auto eFaceL = eCellL.ObtainFace(slotL);
549 int nFN = eFaceL.GetNumNodes();
550
551 std::vector<DNDS::index> nodesL(nFN), nodesR(nFN);
552 std::vector<NodePeriodicBits> pbiL(nFN), pbiR(nFN);
553 eCellL.ExtractFaceNodes(slotL, m->cell2node[iCellL], nodesL);
554 eCellL.ExtractFaceNodes(slotL, m->cell2nodePbi[iCellL], pbiL);
555 eCellR.ExtractFaceNodes(slotR, m->cell2node[iCellR], nodesR);
556 eCellR.ExtractFaceNodes(slotR, m->cell2nodePbi[iCellR], pbiR);
557
558 std::vector<idx_pbi_t> pairsL(nFN), pairsR(nFN);
559 for (int k = 0; k < nFN; k++)
560 {
561 pairsL[k] = {nodesL[k], uint8_t(pbiL[k])};
562 pairsR[k] = {nodesR[k], uint8_t(pbiR[k])};
563 }
564 std::sort(pairsL.begin(), pairsL.end());
565 std::sort(pairsR.begin(), pairsR.end());
566
567 for (int k = 0; k < nFN; k++)
568 {
569 CAPTURE(iFace);
570 CAPTURE(k);
571 CHECK(pairsL[k].first == pairsR[k].first);
572 }
573
574 uint8_t v0 = pairsL[0].second ^ pairsR[0].second;
575 for (int k = 1; k < nFN; k++)
576 {
577 uint8_t vk = pairsL[k].second ^ pairsR[k].second;
578 CAPTURE(iFace);
579 CAPTURE(k);
580 CAPTURE(v0);
581 CAPTURE(vk);
582 CHECK(vk == v0);
583 }
584 }
585 }
586 }
587}
588
589TEST_CASE("Periodic: face2nodePbi is consistent with owner cell's cell2nodePbi")
590{
591 for (int _ci = 0; _ci < N_CONFIGS; _ci++)
592 {
593 if (_ci == 4 && g_mpi.size != 8)
594 continue;
595 const auto &cfg = g_configs[_ci];
596 if (!cfg.periodic)
597 continue;
598 auto m = g_full[_ci];
599 CAPTURE(_ci);
600 SUBCASE(cfg.name)
601 {
602 for (DNDS::index iFace = 0; iFace < m->NumFace(); iFace++)
603 {
604 DNDS::index iCellOwner = m->face2cell(iFace, 0);
605 auto eCell = Elem::Element{m->cellElemInfo[iCellOwner]->getElemType()};
606
607 int iSlot = -1;
608 for (DNDS::rowsize j = 0; j < m->cell2face.RowSize(iCellOwner); j++)
609 if (m->cell2face(iCellOwner, j) == iFace)
610 {
611 iSlot = j;
612 break;
613 }
614 REQUIRE(iSlot >= 0);
615
616 auto eFace = eCell.ObtainFace(iSlot);
617 int nFN = eFace.GetNumNodes();
618
619 std::vector<DNDS::index> cellFaceNodes(nFN);
620 std::vector<NodePeriodicBits> cellFacePbi(nFN);
621 eCell.ExtractFaceNodes(iSlot, m->cell2node[iCellOwner], cellFaceNodes);
622 eCell.ExtractFaceNodes(iSlot, m->cell2nodePbi[iCellOwner], cellFacePbi);
623
624 std::set<idx_pbi_t> fromCell, fromFace;
625 for (int k = 0; k < nFN; k++)
626 {
627 fromCell.insert({cellFaceNodes[k], uint8_t(cellFacePbi[k])});
628 fromFace.insert({m->face2node(iFace, k), uint8_t(m->face2nodePbi(iFace, k))});
629 }
630 CAPTURE(iFace);
631 CHECK(fromCell == fromFace);
632 }
633 }
634 }
635}
636
637TEST_CASE("Periodic: RecreatePeriodicNodes produces consistent counts DSL vs Legacy")
638{
639 for (int ci = 0; ci < N_CONFIGS; ci++)
640 {
641 if (ci == 4 && g_mpi.size != 8)
642 continue;
643 const auto &cfg = g_configs[ci];
644 if (!cfg.periodic)
645 continue;
646 CAPTURE(ci);
647 SUBCASE(cfg.name)
648 {
649 auto mDSL = buildMeshUpToFace(cfg, false);
650 auto mLeg = buildMeshUpToFace(cfg, true);
651
652 mDSL->AdjLocal2GlobalN2CB();
653 mDSL->BuildGhostN2CB();
654 mDSL->AdjGlobal2LocalN2CB();
655 mDSL->RecreatePeriodicNodes();
656
657 mLeg->AdjLocal2GlobalN2CB();
658 mLeg->BuildGhostN2CB();
659 mLeg->AdjGlobal2LocalN2CB();
660 mLeg->RecreatePeriodicNodes();
661
662 CHECK(mDSL->coordsPeriodicRecreated.father->Size() ==
663 mLeg->coordsPeriodicRecreated.father->Size());
664
665 auto dslGlobal = mDSL->coordsPeriodicRecreated.father->globalSize();
666 auto legGlobal = mLeg->coordsPeriodicRecreated.father->globalSize();
667 CHECK(dslGlobal == legGlobal);
668 }
669 }
670}
671
672TEST_CASE("Periodic: face2nodePbi matches DSL vs Legacy")
673{
674 for (int ci = 0; ci < N_CONFIGS; ci++)
675 {
676 if (ci == 4 && g_mpi.size != 8)
677 continue;
678 const auto &cfg = g_configs[ci];
679 if (!cfg.periodic)
680 continue;
681 CAPTURE(ci);
682 SUBCASE(cfg.name)
683 {
684 auto mDSL = buildMeshUpToFace(cfg, false);
685 auto mLeg = buildMeshUpToFace(cfg, true);
686
687 REQUIRE(mDSL->NumFace() == mLeg->NumFace());
688
689 // Map face vertex set → sorted (node, pbi) pairs for each path
690 using FaceKey = std::set<DNDS::index>;
691 using NodePbiSet = std::set<std::pair<DNDS::index, uint8_t>>;
692 std::map<FaceKey, NodePbiSet> dslMap, legMap;
693
694 for (DNDS::index iF = 0; iF < mDSL->NumFace(); iF++)
695 {
696 auto eF = Elem::Element{mDSL->faceElemInfo(iF, 0).getElemType()};
697 FaceKey vs;
698 for (int k = 0; k < eF.GetNumVertices(); k++)
699 vs.insert(mDSL->face2node(iF, k));
700 NodePbiSet nps;
701 for (int k = 0; k < eF.GetNumNodes(); k++)
702 nps.insert({mDSL->face2node(iF, k), uint8_t(mDSL->face2nodePbi(iF, k))});
703 dslMap[vs] = nps;
704 }
705 for (DNDS::index iF = 0; iF < mLeg->NumFace(); iF++)
706 {
707 auto eF = Elem::Element{mLeg->faceElemInfo(iF, 0).getElemType()};
708 FaceKey vs;
709 for (int k = 0; k < eF.GetNumVertices(); k++)
710 vs.insert(mLeg->face2node(iF, k));
711 NodePbiSet nps;
712 for (int k = 0; k < eF.GetNumNodes(); k++)
713 nps.insert({mLeg->face2node(iF, k), uint8_t(mLeg->face2nodePbi(iF, k))});
714 legMap[vs] = nps;
715 }
716
717 CHECK(dslMap.size() == legMap.size());
718 for (auto &[vs, dslNps] : dslMap)
719 {
720 auto it = legMap.find(vs);
721 REQUIRE(it != legMap.end());
722 CHECK(dslNps == it->second);
723 }
724 }
725 }
726}
727
728// --- bnd2cell DSL vs Legacy ---
729
730/// Build mesh through RecoverCell2CellAndBnd2Cell only (not InterpolateFace).
731static ssp<UnstructuredMesh> buildMeshUpToBnd2Cell(
732 const MeshConfig &cfg, bool useLegacy)
733{
734 auto mesh = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
736
737 if (cfg.periodic)
738 {
739 tPoint zero{0, 0, 0};
740 mesh->SetPeriodicGeometry(
741 cfg.translation1, zero, zero,
742 cfg.translation2, zero, zero,
743 cfg.translation3, zero, zero);
744 }
745
746 reader.ReadFromCGNSSerial(meshPath(cfg.file));
747 reader.Deduplicate1to1Periodic(1e-8);
748 reader.BuildCell2Cell();
749
751 pOpt.metisType = "KWAY";
752 pOpt.metisUfactor = 30;
753 pOpt.metisSeed = 42;
754 pOpt.metisNcuts = 1;
755 reader.MeshPartitionCell2Cell(pOpt);
756 reader.PartitionReorderToMeshCell2Cell();
757
758 if (useLegacy)
759 {
760 mesh->RecoverNode2CellAndNode2BndLegacy();
761 mesh->RecoverCell2CellAndBnd2CellLegacy();
762 }
763 else
764 {
765 mesh->RecoverNode2CellAndNode2Bnd();
766 mesh->RecoverCell2CellAndBnd2Cell();
767 }
768
769 return mesh;
770}
771
772TEST_CASE("bnd2cell: DSL matches Legacy on all mesh configs")
773{
774 for (int ci = 0; ci < N_CONFIGS; ci++)
775 {
776 if (ci == 4 && g_mpi.size != 8)
777 continue;
778 const auto &cfg = g_configs[ci];
779 CAPTURE(ci);
780 SUBCASE(cfg.name)
781 {
782 auto mDSL = buildMeshUpToBnd2Cell(cfg, false);
783 auto mLeg = buildMeshUpToBnd2Cell(cfg, true);
784
785 REQUIRE(mDSL->NumBnd() == mLeg->NumBnd());
786
787 // Compare bnd2cell: for each bnd, the cell pair (as a set) must match.
788 for (DNDS::index iBnd = 0; iBnd < mDSL->NumBnd(); iBnd++)
789 {
790 CAPTURE(iBnd);
791 std::set<DNDS::index> dslCells{mDSL->bnd2cell(iBnd, 0), mDSL->bnd2cell(iBnd, 1)};
792 std::set<DNDS::index> legCells{mLeg->bnd2cell(iBnd, 0), mLeg->bnd2cell(iBnd, 1)};
793 CHECK(dslCells == legCells);
794
795 // For periodic boundaries with 2 cells, the ordering matters
796 // (bnd2cell(i,0) is donor-side). Check exact match.
797 if (mDSL->bnd2cell(iBnd, 1) != DNDS::UnInitIndex)
798 {
799 CHECK(mDSL->bnd2cell(iBnd, 0) == mLeg->bnd2cell(iBnd, 0));
800 CHECK(mDSL->bnd2cell(iBnd, 1) == mLeg->bnd2cell(iBnd, 1));
801 }
802 }
803
804 // Also compare cell2cell: for each cell, the neighbor set must match.
805 REQUIRE(mDSL->NumCell() == mLeg->NumCell());
806 for (DNDS::index iCell = 0; iCell < mDSL->NumCell(); iCell++)
807 {
808 CAPTURE(iCell);
809 std::set<DNDS::index> dslNeighbors, legNeighbors;
810 for (DNDS::rowsize j = 0; j < mDSL->cell2cell.father->RowSize(iCell); j++)
811 dslNeighbors.insert(mDSL->cell2cell.father->operator()(iCell, j));
812 for (DNDS::rowsize j = 0; j < mLeg->cell2cell.father->RowSize(iCell); j++)
813 legNeighbors.insert(mLeg->cell2cell.father->operator()(iCell, j));
814 CHECK(dslNeighbors == legNeighbors);
815 }
816 }
817 }
818}
819
820// --- BuildGhostPrimary: DSL vs Legacy ---
821
822/// Build mesh through BuildGhostPrimary, using either DSL or Legacy path.
823static ssp<UnstructuredMesh> buildMeshUpToGhost(
824 const MeshConfig &cfg, bool useLegacy)
825{
826 auto mesh = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
828
829 if (cfg.periodic)
830 {
831 tPoint zero{0, 0, 0};
832 mesh->SetPeriodicGeometry(
833 cfg.translation1, zero, zero,
834 cfg.translation2, zero, zero,
835 cfg.translation3, zero, zero);
836 }
837
838 reader.ReadFromCGNSSerial(meshPath(cfg.file));
839 reader.Deduplicate1to1Periodic(1e-8);
840 reader.BuildCell2Cell();
841
843 pOpt.metisType = "KWAY";
844 pOpt.metisUfactor = 30;
845 pOpt.metisSeed = 42;
846 pOpt.metisNcuts = 1;
847 reader.MeshPartitionCell2Cell(pOpt);
848 reader.PartitionReorderToMeshCell2Cell();
849
850 mesh->RecoverNode2CellAndNode2Bnd();
851 mesh->RecoverCell2CellAndBnd2Cell();
852
853 if (useLegacy)
854 mesh->BuildGhostPrimaryLegacy();
855 else
856 mesh->BuildGhostPrimary();
857
858 return mesh;
859}
860
861TEST_CASE("BuildGhostPrimary: DSL matches Legacy ghost sets")
862{
863 for (int ci = 0; ci < N_CONFIGS; ci++)
864 {
865 if (ci == 4 && g_mpi.size != 8)
866 continue;
867 const auto &cfg = g_configs[ci];
868 CAPTURE(ci);
869 SUBCASE(cfg.name)
870 {
871 auto mDSL = buildMeshUpToGhost(cfg, false);
872 auto mLeg = buildMeshUpToGhost(cfg, true);
873
874 // Cell ghost sets must be identical (same cell2cell traversal).
875 DNDS::index nCellDSL = mDSL->NumCellGhost();
876 DNDS::index nCellLeg = mLeg->NumCellGhost();
877 CHECK(nCellDSL == nCellLeg);
878
879 // Ghost cell indices must match (from the ghost mapping).
880 {
881 auto &dslGhost = mDSL->cell2cell.trans.pLGhostMapping->ghostIndex;
882 auto &legGhost = mLeg->cell2cell.trans.pLGhostMapping->ghostIndex;
883 std::set<DNDS::index> dslSet(dslGhost.begin(), dslGhost.end());
884 std::set<DNDS::index> legSet(legGhost.begin(), legGhost.end());
885 CHECK(dslSet == legSet);
886 }
887
888 // Node ghost: DSL must be a superset of Legacy.
889 // DSL may find additional ghost nodes via its better chain
890 // traversal; Legacy is the baseline.
891 DNDS::index nNodeDSL = mDSL->NumNodeGhost();
892 DNDS::index nNodeLeg = mLeg->NumNodeGhost();
893 CHECK(nNodeDSL >= nNodeLeg);
894 // CHECK(nNodeDSL == nNodeLeg);
895
896 {
897 auto &dslGhost = mDSL->coords.trans.pLGhostMapping->ghostIndex;
898 auto &legGhost = mLeg->coords.trans.pLGhostMapping->ghostIndex;
899 std::set<DNDS::index> dslSet(dslGhost.begin(), dslGhost.end());
900 std::set<DNDS::index> legSet(legGhost.begin(), legGhost.end());
901 CHECK(std::includes(dslSet.begin(), dslSet.end(),
902 legSet.begin(), legSet.end()));
903 // CHECK(dslSet == legSet);
904 }
905
906 // Bnd ghost: DSL must be a superset of Legacy.
907 // DSL finds cross-rank bnd entries that Legacy misses because
908 // Legacy iterates node2bnd.father only (built from bnd2node.father).
909 DNDS::index nBndDSL = mDSL->NumBndGhost();
910 DNDS::index nBndLeg = mLeg->NumBndGhost();
911 CHECK(nBndDSL >= nBndLeg);
912 // CHECK(nBndDSL == nBndLeg);
913
914 {
915 auto &dslGhost = mDSL->bnd2cell.trans.pLGhostMapping->ghostIndex;
916 auto &legGhost = mLeg->bnd2cell.trans.pLGhostMapping->ghostIndex;
917 std::set<DNDS::index> dslSet(dslGhost.begin(), dslGhost.end());
918 std::set<DNDS::index> legSet(legGhost.begin(), legGhost.end());
919 CHECK(std::includes(dslSet.begin(), dslSet.end(),
920 legSet.begin(), legSet.end()));
921 // CHECK(dslSet == legSet);
922 }
923
924 if (g_mpi.rank == 0)
925 fmt::print(" [{}] ghost: cells={}, nodes={}, bnds={}\n",
926 cfg.name, nCellDSL, nNodeDSL, nBndDSL);
927
928 if (ci == 0 || ci == 2)
929 {
930 MPI_Barrier(g_mpi.comm);
931 fmt::print(" [{}] r{} DSL nBndGhost={} LEG nBndGhost={}\n",
932 cfg.name, g_mpi.rank, nBndDSL, nBndLeg);
933 // Ghost bnds: global bnd index and their bnd2node entries (global node indices)
934 for (DNDS::index iBnd = mDSL->bnd2node.father->Size(); iBnd < mDSL->bnd2node.Size(); iBnd++)
935 {
936 DNDS::index iSon = iBnd - mDSL->bnd2node.father->Size();
937 DNDS::index bndGlobal = mDSL->bnd2node.trans.pLGhostMapping->ghostIndex[iSon];
938 fmt::print(" r{} gBnd[{}]:", g_mpi.rank, bndGlobal);
939 for (DNDS::rowsize j = 0; j < mDSL->bnd2node.RowSize(iBnd); j++)
940 fmt::print(" {}", mDSL->bnd2node(iBnd, j));
941 fmt::print("\n");
942 }
943 }
944 }
945 }
946}
947
948// --- Benchmark: DSL vs Legacy for all migrated APIs ---
949
950/// Build a mesh up to the point where we can benchmark each API independently.
951/// Returns mesh in Adj_PointToGlobal state after serial read + partition.
952static ssp<UnstructuredMesh> buildMeshForBenchmark(const MeshConfig &cfg)
953{
954 auto mesh = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
956
957 if (cfg.periodic)
958 {
959 tPoint zero{0, 0, 0};
960 mesh->SetPeriodicGeometry(
961 cfg.translation1, zero, zero,
962 cfg.translation2, zero, zero,
963 cfg.translation3, zero, zero);
964 }
965
966 reader.ReadFromCGNSSerial(meshPath(cfg.file));
967 reader.Deduplicate1to1Periodic(1e-8);
968 reader.BuildCell2Cell();
969
971 pOpt.metisType = "KWAY";
972 pOpt.metisUfactor = 30;
973 pOpt.metisSeed = 42;
974 pOpt.metisNcuts = 1;
975 reader.MeshPartitionCell2Cell(pOpt);
976 reader.PartitionReorderToMeshCell2Cell();
977 return mesh;
978}
979
980TEST_CASE("Benchmark: DSL vs Legacy pipeline steps")
981{
982 // Use NACA0012_H2 (20816 cells) — large enough to measure, small enough to
983 // run quickly. Ball2 is too expensive and the small meshes have too much
984 // noise. Index 2 in g_configs.
985 const auto &cfg = g_configs[2];
986
987 if (g_mpi.size < 2)
988 return; // benchmarking single-rank is not interesting
989
990 const int nRuns = 3; // repeat for stability
991
992 // --- RecoverNode2CellAndNode2Bnd (Inverse) ---
993 {
994 double tDSL = 0, tLeg = 0;
995 for (int r = 0; r < nRuns; r++)
996 {
997 auto m1 = buildMeshForBenchmark(cfg);
998 MPI_Barrier(g_mpi.comm);
999 double t0 = MPI_Wtime();
1000 m1->RecoverNode2CellAndNode2Bnd();
1001 MPI_Barrier(g_mpi.comm);
1002 double t1 = MPI_Wtime();
1003 tDSL += t1 - t0;
1004
1005 auto m2 = buildMeshForBenchmark(cfg);
1006 MPI_Barrier(g_mpi.comm);
1007 double t2 = MPI_Wtime();
1008 m2->RecoverNode2CellAndNode2BndLegacy();
1009 MPI_Barrier(g_mpi.comm);
1010 double t3 = MPI_Wtime();
1011 tLeg += t3 - t2;
1012 }
1013 if (g_mpi.rank == 0)
1014 fmt::print(" Inverse (RecoverN2CB): DSL {:.4f}s Legacy {:.4f}s ({:.2f}x, {} runs)\n",
1015 tDSL / nRuns, tLeg / nRuns, tLeg / tDSL, nRuns);
1016 }
1017
1018 // --- RecoverCell2CellAndBnd2Cell (ComposeFiltered) ---
1019 {
1020 double tDSL = 0, tLeg = 0;
1021 for (int r = 0; r < nRuns; r++)
1022 {
1023 auto m1 = buildMeshForBenchmark(cfg);
1024 m1->RecoverNode2CellAndNode2Bnd();
1025 MPI_Barrier(g_mpi.comm);
1026 double t0 = MPI_Wtime();
1027 m1->RecoverCell2CellAndBnd2Cell();
1028 MPI_Barrier(g_mpi.comm);
1029 double t1 = MPI_Wtime();
1030 tDSL += t1 - t0;
1031
1032 auto m2 = buildMeshForBenchmark(cfg);
1033 m2->RecoverNode2CellAndNode2BndLegacy();
1034 MPI_Barrier(g_mpi.comm);
1035 double t2 = MPI_Wtime();
1036 m2->RecoverCell2CellAndBnd2CellLegacy();
1037 MPI_Barrier(g_mpi.comm);
1038 double t3 = MPI_Wtime();
1039 tLeg += t3 - t2;
1040 }
1041 if (g_mpi.rank == 0)
1042 fmt::print(" Compose (RecoverC2C): DSL {:.4f}s Legacy {:.4f}s ({:.2f}x, {} runs)\n",
1043 tDSL / nRuns, tLeg / nRuns, tLeg / tDSL, nRuns);
1044 }
1045
1046 // --- BuildGhostPrimary (evaluateGhostTree for cell ghost) ---
1047 {
1048 double tDSL = 0, tLeg = 0;
1049 for (int r = 0; r < nRuns; r++)
1050 {
1051 auto m1 = buildMeshForBenchmark(cfg);
1052 m1->RecoverNode2CellAndNode2Bnd();
1053 m1->RecoverCell2CellAndBnd2Cell();
1054 MPI_Barrier(g_mpi.comm);
1055 double t0 = MPI_Wtime();
1056 m1->BuildGhostPrimary();
1057 MPI_Barrier(g_mpi.comm);
1058 double t1 = MPI_Wtime();
1059 tDSL += t1 - t0;
1060
1061 auto m2 = buildMeshForBenchmark(cfg);
1062 m2->RecoverNode2CellAndNode2Bnd();
1063 m2->RecoverCell2CellAndBnd2Cell();
1064 MPI_Barrier(g_mpi.comm);
1065 double t2 = MPI_Wtime();
1066 m2->BuildGhostPrimaryLegacy();
1067 MPI_Barrier(g_mpi.comm);
1068 double t3 = MPI_Wtime();
1069 tLeg += t3 - t2;
1070 }
1071 if (g_mpi.rank == 0)
1072 fmt::print(" Ghost (BuildGhostPri): DSL {:.4f}s Legacy {:.4f}s ({:.2f}x, {} runs)\n",
1073 tDSL / nRuns, tLeg / nRuns, tLeg / tDSL, nRuns);
1074 }
1075
1076 // --- InterpolateFace (Interpolate) ---
1077 {
1078 double tDSL = 0, tLeg = 0;
1079 for (int r = 0; r < nRuns; r++)
1080 {
1081 auto m1 = buildMeshForBenchmark(cfg);
1082 m1->RecoverNode2CellAndNode2Bnd();
1083 m1->RecoverCell2CellAndBnd2Cell();
1084 m1->BuildGhostPrimary();
1085 m1->AdjGlobal2LocalPrimary();
1086 m1->AdjGlobal2LocalN2CB();
1087 MPI_Barrier(g_mpi.comm);
1088 double t0 = MPI_Wtime();
1089 m1->InterpolateFace();
1090 MPI_Barrier(g_mpi.comm);
1091 double t1 = MPI_Wtime();
1092 tDSL += t1 - t0;
1093
1094 auto m2 = buildMeshForBenchmark(cfg);
1095 m2->RecoverNode2CellAndNode2Bnd();
1096 m2->RecoverCell2CellAndBnd2Cell();
1097 m2->BuildGhostPrimary();
1098 m2->AdjGlobal2LocalPrimary();
1099 m2->AdjGlobal2LocalN2CB();
1100 MPI_Barrier(g_mpi.comm);
1101 double t2 = MPI_Wtime();
1102 m2->InterpolateFaceLegacy();
1103 MPI_Barrier(g_mpi.comm);
1104 double t3 = MPI_Wtime();
1105 tLeg += t3 - t2;
1106 }
1107 if (g_mpi.rank == 0)
1108 fmt::print(" Interpolate (Face): DSL {:.4f}s Legacy {:.4f}s ({:.2f}x, {} runs)\n",
1109 tDSL / nRuns, tLeg / nRuns, tLeg / tDSL, nRuns);
1110 }
1111}
1112
1113// --- N2CB ---
1114
1115TEST_CASE("N2CB: every local node has at least one adjacent cell"){
1117 for (DNDS::index iNode = 0; iNode < m->NumNode(); iNode++)
1118 {
1119 bool hasCell = false;
1120 for (DNDS::rowsize j = 0; j < m->node2cell.RowSize(iNode); j++)
1121 if (m->node2cell(iNode, j) >= 0)
1122 hasCell = true;
1123 CHECK(hasCell);
1124 }
1125 })}
1126
1127TEST_CASE("N2CB: node2cell entries are valid local indices"){
1129 DNDS::index totalCells = m->NumCell() + m->NumCellGhost();
1130 for (DNDS::index iNode = 0; iNode < m->NumNodeProc(); iNode++)
1131 for (DNDS::rowsize j = 0; j < m->node2cell.RowSize(iNode); j++)
1132 {
1133 DNDS::index iC = m->node2cell(iNode, j);
1134 if (iC >= 0)
1135 CHECK(iC < totalCells);
1136 }
1137 })}
1138
1139// --- BuildCell2CellFace ---
1140
1141TEST_CASE("BuildCell2CellFace: row sizes match cell2face"){
1143 for (DNDS::index iC = 0; iC < m->NumCell(); iC++)
1144 CHECK(m->cell2cellFace.RowSize(iC) == m->cell2face.RowSize(iC));
1145 })}
1146
1147TEST_CASE("BuildCell2CellFace: entries are valid cells or negative"){
1149 DNDS::index totalCells = m->NumCell() + m->NumCellGhost();
1150 for (DNDS::index iC = 0; iC < m->NumCell(); iC++)
1151 for (DNDS::rowsize j = 0; j < m->cell2cellFace.RowSize(iC); j++)
1152 {
1153 DNDS::index nb = m->cell2cellFace(iC, j);
1154 if (nb >= 0)
1155 CHECK(nb < totalCells);
1156 }
1157 })}
1158
1159// --- ConstructBndMesh ---
1160
1161TEST_CASE("ConstructBndMesh: correct dimensions and state"){
1163 if (cfg.periodic)
1164 continue; // periodic meshes may have no real boundaries after deduplication
1165 UnstructuredMesh bMesh(g_mpi, m->dim - 1);
1166 m->ConstructBndMesh(bMesh);
1167 CHECK(bMesh.dim == m->dim - 1);
1168 CHECK(bMesh.NumCell() == m->NumBnd());
1169 CHECK(bMesh.adjPrimaryState == Adj_PointToLocal);
1170 })}
1171
1172TEST_CASE("ConstructBndMesh: node2parentNode valid"){
1174 if (cfg.periodic)
1175 continue;
1176 UnstructuredMesh bMesh(g_mpi, m->dim - 1);
1177 m->ConstructBndMesh(bMesh);
1178 CHECK(bMesh.node2parentNode.size() ==
1179 static_cast<size_t>(bMesh.NumNode() + bMesh.NumNodeGhost()));
1180 for (DNDS::index i = 0; i < bMesh.NumNode(); i++)
1181 {
1182 DNDS::index pn = bMesh.node2parentNode[i];
1183 CHECK(pn >= 0);
1184 CHECK(pn < m->NumNode() + m->NumNodeGhost());
1185 }
1186 })}
1187
1188// --- BuildVTKConnectivity ---
1189
1190TEST_CASE("BuildVTKConnectivity: arrays populated correctly"){
1192 CHECK(m->vtkCell2node.size() > 0);
1193 CHECK(m->vtkCell2nodeOffsets.size() == static_cast<size_t>(m->NumCell() + 1));
1194 CHECK(m->vtkCellType.size() == static_cast<size_t>(m->NumCell()));
1195 for (size_t i = 1; i < m->vtkCell2nodeOffsets.size(); i++)
1196 CHECK(m->vtkCell2nodeOffsets[i] >= m->vtkCell2nodeOffsets[i - 1]);
1197 for (auto ct : m->vtkCellType)
1198 CHECK(ct != 0);
1199 })}
1200
1201// ===========================================================================
1202// ReorderLocalCells (config 0 only)
1203// ===========================================================================
1204
1205TEST_CASE("ReorderLocalCells: all states preserved"){
1207 m->ReorderLocalCells(2);
1208 CHECK(m->adjPrimaryState == Adj_PointToLocal);
1209 CHECK(m->adjFacialState == Adj_PointToLocal);
1210 CHECK(m->adjC2FState == Adj_PointToLocal);
1211 CHECK(m->adjN2CBState == Adj_PointToLocal);
1212 CHECK(m->adjC2CFaceState == Adj_PointToLocal);
1213 })}
1214
1215TEST_CASE("ReorderLocalCells: cell count preserved")
1216{
1218 DNDS::index g1 = 0;
1219 DNDS::index g2 = 0;
1220 DNDS::index loc1 = m->NumCell();
1221 MPI_Allreduce(&loc1, &g1, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
1222 m->ReorderLocalCells(2);
1223 DNDS::index loc2 = m->NumCell();
1224 MPI_Allreduce(&loc2, &g2, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
1225 CHECK(g1 == g2);)
1226}
1227
1228TEST_CASE("ReorderLocalCells: partition starts valid"){
1230 m->ReorderLocalCells(2);
1231 int nParts = m->NLocalParts();
1232 CHECK(nParts == 2);
1233 CHECK(m->LocalPartStart(0) == 0);
1234 CHECK(m->LocalPartEnd(nParts - 1) == m->NumCell());
1235 for (int ip = 0; ip < nParts; ip++)
1236 CHECK(m->LocalPartStart(ip) <= m->LocalPartEnd(ip));
1237 })}
1238
1239TEST_CASE("ReorderLocalCells: cell2node still valid"){
1241 m->ReorderLocalCells(2);
1242 DNDS::index totalNodes = m->NumNode() + m->NumNodeGhost();
1243 for (DNDS::index iC = 0; iC < m->NumCell(); iC++)
1244 for (DNDS::rowsize j = 0; j < m->cell2node.RowSize(iC); j++)
1245 {
1246 DNDS::index iN = m->cell2node(iC, j);
1247 CHECK(iN >= 0);
1248 CHECK(iN < totalNodes);
1249 }
1250 })}
1251
1252TEST_CASE("ReorderLocalCells: face count preserved")
1253{
1255 DNDS::index g1 = 0;
1256 DNDS::index g2 = 0;
1257 DNDS::index loc1 = m->NumFace();
1258 MPI_Allreduce(&loc1, &g1, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
1259 m->ReorderLocalCells(2);
1260 DNDS::index loc2 = m->NumFace();
1261 MPI_Allreduce(&loc2, &g2, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
1262 CHECK(g1 == g2);)
1263}
1264
1265// ===========================================================================
1266// Adj round-trip tests (all mesh configs)
1267// ===========================================================================
1268
1269TEST_CASE("AdjFacial round-trip: face2node and face2cell identity")
1270{
1272 REQUIRE(m->adjFacialState == Adj_PointToLocal);
1273
1274 auto f2nSnap = snapshotAdj(m->face2node, m->NumFace());
1275 auto f2cSnap = snapshotAdj2(m->face2cell, m->NumFace());
1276
1277 m->AdjLocal2GlobalFacial();
1278 CHECK(m->adjFacialState == Adj_PointToGlobal);
1279
1280 for (DNDS::index iF = 0; iF < m->NumFace(); iF++) for (DNDS::rowsize j = 0; j < m->face2node.RowSize(iF); j++)
1281 CHECK(m->face2node(iF, j) >= 0);
1282
1283 m->AdjGlobal2LocalFacial();
1284 CHECK(m->adjFacialState == Adj_PointToLocal);
1285
1286 checkAdjMatchesSnapshot(m->face2node, m->NumFace(), f2nSnap);
1287 for (DNDS::index iF = 0; iF < m->NumFace(); iF++) {
1288 CHECK(m->face2cell(iF, 0) == f2cSnap[iF].first);
1289 CHECK(m->face2cell(iF, 1) == f2cSnap[iF].second);
1290 }
1291 // NOTE: face2bnd intentionally not round-tripped (asymmetric by design)
1292 )
1293}
1294
1295TEST_CASE("AdjC2F round-trip: cell2face and bnd2face identity")
1296{
1298 REQUIRE(m->adjC2FState == Adj_PointToLocal);
1299
1300 auto c2fSnap = snapshotAdj(m->cell2face, m->NumCell());
1301 auto b2fSnap = snapshotAdj1(m->bnd2face, m->NumBnd());
1302
1303 m->AdjLocal2GlobalC2F();
1304 CHECK(m->adjC2FState == Adj_PointToGlobal);
1305
1306 m->AdjGlobal2LocalC2F();
1307 CHECK(m->adjC2FState == Adj_PointToLocal);
1308
1309 checkAdjMatchesSnapshot(m->cell2face, m->NumCell(), c2fSnap);
1310 for (DNDS::index ib = 0; ib < m->NumBnd(); ib++)
1311 CHECK(m->bnd2face(ib, 0) == b2fSnap[ib]);)
1312}
1313
1314TEST_CASE("AdjN2CB round-trip: node2cell and node2bnd identity")
1315{
1317 REQUIRE(m->adjN2CBState == Adj_PointToLocal);
1318
1319 auto n2cSnap = snapshotAdj(m->node2cell, m->NumNodeProc());
1320 auto n2bSnap = snapshotAdj(m->node2bnd, m->NumNodeProc());
1321
1322 m->AdjLocal2GlobalN2CB();
1323 CHECK(m->adjN2CBState == Adj_PointToGlobal);
1324
1325 m->AdjGlobal2LocalN2CB();
1326 CHECK(m->adjN2CBState == Adj_PointToLocal);
1327
1328 checkAdjMatchesSnapshot(m->node2cell, m->NumNodeProc(), n2cSnap);
1329 checkAdjMatchesSnapshot(m->node2bnd, m->NumNodeProc(), n2bSnap);)
1330}
1331
1332TEST_CASE("AdjC2CFace round-trip: cell2cellFace identity")
1333{
1335 REQUIRE(m->adjC2CFaceState == Adj_PointToLocal);
1336
1337 auto c2cfSnap = snapshotAdj(m->cell2cellFace, m->NumCell());
1338
1339 m->AdjLocal2GlobalC2CFace();
1340 CHECK(m->adjC2CFaceState == Adj_PointToGlobal);
1341
1342 m->AdjGlobal2LocalC2CFace();
1343 CHECK(m->adjC2CFaceState == Adj_PointToLocal);
1344
1345 checkAdjMatchesSnapshot(m->cell2cellFace, m->NumCell(), c2cfSnap);)
1346}
1347
1348TEST_CASE("AdjPrimary round-trip: serial-out pattern")
1349{
1351 REQUIRE(m->adjPrimaryState == Adj_PointToLocal);
1352
1353 auto c2nSnap = snapshotAdj(m->cell2node, m->NumCell());
1354 auto c2cSnap = snapshotAdj(m->cell2cell, m->NumCell());
1355 auto b2nSnap = snapshotAdj(m->bnd2node, m->NumBnd());
1356 auto b2cSnap = snapshotAdj2(m->bnd2cell, m->NumBnd());
1357
1358 m->AdjLocal2GlobalPrimary();
1359 CHECK(m->adjPrimaryState == Adj_PointToGlobal);
1360
1361 DNDS::index globalNodeSize = m->coords.father->globalSize();
1362 for (DNDS::index iC = 0; iC < m->NumCell(); iC++) for (DNDS::rowsize j = 0; j < m->cell2node.RowSize(iC); j++) {
1363 DNDS::index gN = m->cell2node(iC, j);
1364 CHECK(gN >= 0);
1365 CHECK(gN < globalNodeSize);
1366 }
1367
1368 m->AdjGlobal2LocalPrimary();
1369 CHECK(m->adjPrimaryState == Adj_PointToLocal);
1370
1371 checkAdjMatchesSnapshot(m->cell2node, m->NumCell(), c2nSnap);
1372 checkAdjMatchesSnapshot(m->cell2cell, m->NumCell(), c2cSnap);
1373 checkAdjMatchesSnapshot(m->bnd2node, m->NumBnd(), b2nSnap);
1374 for (DNDS::index ib = 0; ib < m->NumBnd(); ib++) {
1375 CHECK(m->bnd2cell(ib, 0) == b2cSnap[ib].first);
1376 CHECK(m->bnd2cell(ib, 1) == b2cSnap[ib].second);
1377 })
1378}
1379
1380TEST_CASE("AdjForBnd round-trip: ConstructBndMesh + ForBnd")
1381{
1383 UnstructuredMesh bMesh(g_mpi, m->dim - 1);
1384 m->ConstructBndMesh(bMesh);
1385 REQUIRE(bMesh.adjPrimaryState == Adj_PointToLocal);
1386
1387 auto c2nSnap = snapshotAdj(bMesh.cell2node, bMesh.NumCell());
1388
1389 bMesh.AdjLocal2GlobalPrimaryForBnd();
1390 CHECK(bMesh.adjPrimaryState == Adj_PointToGlobal);
1391
1392 DNDS::index globalNodeSize = bMesh.coords.father->globalSize();
1393 for (DNDS::index iC = 0; iC < bMesh.NumCell(); iC++) for (DNDS::rowsize j = 0; j < bMesh.cell2node.RowSize(iC); j++) {
1394 DNDS::index gN = bMesh.cell2node(iC, j);
1395 CHECK(gN >= 0);
1396 CHECK(gN < globalNodeSize);
1397 }
1398
1399 bMesh.AdjGlobal2LocalPrimaryForBnd();
1400 CHECK(bMesh.adjPrimaryState == Adj_PointToLocal);
1401
1402 checkAdjMatchesSnapshot(bMesh.cell2node, bMesh.NumCell(), c2nSnap);)
1403}
1404
1405// ===========================================================================
1406// Ball2: 3D mixed-element mesh (np=8 only)
1407// ===========================================================================
1408
1409TEST_CASE("Ball2: mixed 3D element types")
1410{
1411 // Ball2 (config 4) is only tested with np=8
1412 if (g_mpi.size != 8)
1413 return;
1414
1415 auto m = g_full[4]; // Ball2 is index 4
1416 DNDS_assert(m != nullptr);
1417
1418 // Count element types locally
1419 std::map<Elem::ElemType, int> typeCounts;
1420 for (DNDS::index iC = 0; iC < m->NumCell(); iC++)
1421 {
1422 auto elem = m->GetCellElement(iC);
1423 typeCounts[elem.type]++;
1424 }
1425
1426 // Gather all unique element types from all ranks
1427 std::vector<int> localTypes;
1428 for (const auto &[type, count] : typeCounts)
1429 localTypes.push_back(static_cast<int>(type));
1430
1431 int nLocalTypes = localTypes.size();
1432 std::vector<int> allRankSizes(g_mpi.size);
1433 MPI_Gather(&nLocalTypes, 1, MPI_INT, allRankSizes.data(), 1, MPI_INT, 0, g_mpi.comm);
1434
1435 std::vector<int> allTypes;
1436 std::vector<int> displs;
1437 int totalTypes = 0;
1438
1439 if (g_mpi.rank == 0)
1440 {
1441 displs.resize(g_mpi.size);
1442 for (int i = 0; i < g_mpi.size; i++)
1443 {
1444 displs[i] = totalTypes;
1445 totalTypes += allRankSizes[i];
1446 }
1447 allTypes.resize(totalTypes);
1448 }
1449
1450 MPI_Gatherv(localTypes.data(), nLocalTypes, MPI_INT,
1451 allTypes.data(), allRankSizes.data(), displs.data(), MPI_INT, 0, g_mpi.comm);
1452
1453 // Build set of unique types (on rank 0)
1454 std::set<Elem::ElemType> uniqueTypes;
1455 if (g_mpi.rank == 0)
1456 {
1457 for (int typeInt : allTypes)
1458 uniqueTypes.insert(static_cast<Elem::ElemType>(typeInt));
1459 }
1460
1461 // Broadcast unique types to all ranks
1462 int nUnique = uniqueTypes.size();
1463 MPI_Bcast(&nUnique, 1, MPI_INT, 0, g_mpi.comm);
1464 std::vector<int> uniqueTypesVec;
1465 if (g_mpi.rank == 0)
1466 {
1467 for (auto type : uniqueTypes)
1468 uniqueTypesVec.push_back(static_cast<int>(type));
1469 }
1470 uniqueTypesVec.resize(nUnique);
1471 MPI_Bcast(uniqueTypesVec.data(), nUnique, MPI_INT, 0, g_mpi.comm);
1472
1473 // All ranks iterate over the same types in the same order
1474 std::map<Elem::ElemType, int> globalCounts;
1475 for (int typeInt : uniqueTypesVec)
1476 {
1477 Elem::ElemType type = static_cast<Elem::ElemType>(typeInt);
1478 int localCount = typeCounts.count(type) ? typeCounts[type] : 0;
1479 int globalCount;
1480 MPI_Reduce(&localCount, &globalCount, 1, MPI_INT, MPI_SUM, 0, g_mpi.comm);
1481 if (g_mpi.rank == 0)
1482 globalCounts[type] = globalCount;
1483 }
1484
1485 if (g_mpi.rank == 0)
1486 {
1487 std::cout << "\nBall2 element type counts:" << std::endl;
1488 for (const auto &[type, count] : globalCounts)
1489 {
1490 std::cout << " Type " << static_cast<int>(type) << ": " << count << std::endl;
1491 }
1492
1493 // Verify expected total
1494 int total = 0;
1495 for (const auto &[type, count] : globalCounts)
1496 total += count;
1497 CHECK(total == 958994);
1498 }
1499}
1500
1501// ===========================================================================
1502// Elevation/Bisection tests
1503// ===========================================================================
1504
1505TEST_CASE("Elevation: O2 mesh cell count equals O1 cell count")
1506{
1507 // Only test with config 0 (UniformSquare_10) for now
1508 const auto &cfg = g_configs[0];
1509
1510 // Build fresh O1 mesh
1511 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1512 UnstructuredMeshSerialRW reader(mO1, 0);
1513
1514 if (cfg.periodic)
1515 {
1516 tPoint zero{0, 0, 0};
1517 mO1->SetPeriodicGeometry(cfg.translation1, zero, zero, cfg.translation2, zero, zero, cfg.translation3, zero, zero);
1518 }
1519
1520 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1521 reader.Deduplicate1to1Periodic(1e-8);
1522 reader.BuildCell2Cell();
1523
1525 pOpt.metisType = "KWAY";
1526 pOpt.metisUfactor = 30;
1527 pOpt.metisSeed = 42;
1528 pOpt.metisNcuts = 1;
1529 reader.MeshPartitionCell2Cell(pOpt);
1530 reader.PartitionReorderToMeshCell2Cell();
1531
1532 mO1->RecoverNode2CellAndNode2Bnd();
1533 mO1->RecoverCell2CellAndBnd2Cell();
1534 mO1->BuildGhostPrimary();
1535 mO1->AdjGlobal2LocalPrimary();
1536
1537 // Build O2 mesh via elevation
1538 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1539 mO2->BuildO2FromO1Elevation(*mO1);
1540
1541 // Set up MPI state for mO2
1542 mO2->RecoverNode2CellAndNode2Bnd();
1543 mO2->RecoverCell2CellAndBnd2Cell();
1544 mO2->BuildGhostPrimary();
1545 mO2->AdjGlobal2LocalPrimary();
1546
1547 DNDS::index nCellO1 = mO1->NumCellGlobal();
1548 DNDS::index nCellO2 = mO2->NumCellGlobal();
1549
1550 // Elevation preserves cell count
1551 CHECK(nCellO2 == nCellO1);
1552
1553 // All cells should be O2 type
1554 for (DNDS::index iC = 0; iC < mO2->NumCell(); iC++)
1555 {
1556 auto elem = mO2->GetCellElement(iC);
1557 CHECK(elem.GetOrder() == 2);
1558 }
1559}
1560
1561TEST_CASE("Elevation: O2 mesh has more nodes than O1")
1562{
1563 const auto &cfg = g_configs[0];
1564
1565 // Build fresh O1 mesh
1566 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1567 UnstructuredMeshSerialRW reader(mO1, 0);
1568
1569 if (cfg.periodic)
1570 {
1571 tPoint zero{0, 0, 0};
1572 mO1->SetPeriodicGeometry(cfg.translation1, zero, zero, cfg.translation2, zero, zero, cfg.translation3, zero, zero);
1573 }
1574
1575 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1576 reader.Deduplicate1to1Periodic(1e-8);
1577 reader.BuildCell2Cell();
1578
1580 pOpt.metisType = "KWAY";
1581 pOpt.metisUfactor = 30;
1582 pOpt.metisSeed = 42;
1583 pOpt.metisNcuts = 1;
1584 reader.MeshPartitionCell2Cell(pOpt);
1585 reader.PartitionReorderToMeshCell2Cell();
1586
1587 mO1->RecoverNode2CellAndNode2Bnd();
1588 mO1->RecoverCell2CellAndBnd2Cell();
1589 mO1->BuildGhostPrimary();
1590 mO1->AdjGlobal2LocalPrimary();
1591
1592 // Build O2 mesh via elevation
1593 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1594 mO2->BuildO2FromO1Elevation(*mO1);
1595
1596 // Set up MPI state for mO2
1597 mO2->RecoverNode2CellAndNode2Bnd();
1598 mO2->RecoverCell2CellAndBnd2Cell();
1599 mO2->BuildGhostPrimary();
1600 mO2->AdjGlobal2LocalPrimary();
1601
1602 DNDS::index nNodeO1 = mO1->coords.father->globalSize();
1603 DNDS::index nNodeO2 = mO2->coords.father->globalSize();
1604
1605 // Elevation adds nodes
1606 CHECK(nNodeO2 > nNodeO1);
1607}
1608
1609TEST_CASE("Bisection: O1 mesh has more cells than O2")
1610{
1611 const auto &cfg = g_configs[0];
1612
1613 // Build fresh O1 mesh
1614 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1615 UnstructuredMeshSerialRW reader(mO1, 0);
1616
1617 if (cfg.periodic)
1618 {
1619 tPoint zero{0, 0, 0};
1620 mO1->SetPeriodicGeometry(cfg.translation1, zero, zero, cfg.translation2, zero, zero, cfg.translation3, zero, zero);
1621 }
1622
1623 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1624 reader.Deduplicate1to1Periodic(1e-8);
1625 reader.BuildCell2Cell();
1626
1628 pOpt.metisType = "KWAY";
1629 pOpt.metisUfactor = 30;
1630 pOpt.metisSeed = 42;
1631 pOpt.metisNcuts = 1;
1632 reader.MeshPartitionCell2Cell(pOpt);
1633 reader.PartitionReorderToMeshCell2Cell();
1634
1635 mO1->RecoverNode2CellAndNode2Bnd();
1636 mO1->RecoverCell2CellAndBnd2Cell();
1637 mO1->BuildGhostPrimary();
1638 mO1->AdjGlobal2LocalPrimary();
1639
1640 // Build O2 mesh
1641 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1642 mO2->BuildO2FromO1Elevation(*mO1);
1643 mO2->RecoverNode2CellAndNode2Bnd();
1644 mO2->RecoverCell2CellAndBnd2Cell();
1645 mO2->BuildGhostPrimary();
1646 mO2->AdjGlobal2LocalPrimary();
1647
1648 // Build O1 mesh via bisection
1649 mO2->AdjLocal2GlobalPrimary();
1650 auto mO1B = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1651 mO1B->BuildBisectO1FormO2(*mO2);
1652
1653 DNDS::index nCellO2 = mO2->NumCellGlobal();
1654 DNDS::index nCellO1B = mO1B->NumCellGlobal();
1655
1656 // Bisection increases cell count
1657 CHECK(nCellO1B > nCellO2);
1658
1659 // All cells should be O1 (linear) element types
1660 for (DNDS::index iC = 0; iC < mO1B->NumCell(); iC++)
1661 {
1662 auto elem = mO1B->GetCellElement(iC);
1663 CHECK(elem.GetOrder() == 1);
1664 }
1665}
1666
1667TEST_CASE("Bisection: global node count is preserved from O2 mesh")
1668{
1669 const auto &cfg = g_configs[0];
1670
1671 // Build fresh O1 mesh
1672 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1673 UnstructuredMeshSerialRW reader(mO1, 0);
1674
1675 if (cfg.periodic)
1676 {
1677 tPoint zero{0, 0, 0};
1678 mO1->SetPeriodicGeometry(cfg.translation1, zero, zero, cfg.translation2, zero, zero, cfg.translation3, zero, zero);
1679 }
1680
1681 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1682 reader.Deduplicate1to1Periodic(1e-8);
1683 reader.BuildCell2Cell();
1684
1686 pOpt.metisType = "KWAY";
1687 pOpt.metisUfactor = 30;
1688 pOpt.metisSeed = 42;
1689 pOpt.metisNcuts = 1;
1690 reader.MeshPartitionCell2Cell(pOpt);
1691 reader.PartitionReorderToMeshCell2Cell();
1692
1693 mO1->RecoverNode2CellAndNode2Bnd();
1694 mO1->RecoverCell2CellAndBnd2Cell();
1695 mO1->BuildGhostPrimary();
1696 mO1->AdjGlobal2LocalPrimary();
1697
1698 // Build O2 mesh
1699 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1700 mO2->BuildO2FromO1Elevation(*mO1);
1701 mO2->RecoverNode2CellAndNode2Bnd();
1702 mO2->RecoverCell2CellAndBnd2Cell();
1703 mO2->BuildGhostPrimary();
1704 mO2->AdjGlobal2LocalPrimary();
1705
1706 // Build O1 mesh via bisection
1707 mO2->AdjLocal2GlobalPrimary();
1708 auto mO1B = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1709 mO1B->BuildBisectO1FormO2(*mO2);
1710
1711 // No new nodes created during bisection
1712 DNDS::index nNodeO2 = mO2->coords.father->globalSize();
1713 DNDS::index nNodeO1B = mO1B->coords.father->globalSize();
1714 CHECK(nNodeO1B == nNodeO2);
1715}
1716
1717TEST_CASE("Elevation+Bisection: exact cell count progression")
1718{
1719 // UniformSquare_10: 100 Quad4 cells
1720 // After elevation: 100 Quad9 cells (same count, higher order)
1721 // After bisection: 400 Quad4 cells (4 sub-cells per Quad9)
1722 const auto &cfg = g_configs[0];
1723
1724 // Build fresh O1 mesh
1725 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1726 UnstructuredMeshSerialRW reader(mO1, 0);
1727
1728 if (cfg.periodic)
1729 {
1730 tPoint zero{0, 0, 0};
1731 mO1->SetPeriodicGeometry(cfg.translation1, zero, zero, cfg.translation2, zero, zero, cfg.translation3, zero, zero);
1732 }
1733
1734 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1735 reader.Deduplicate1to1Periodic(1e-8);
1736 reader.BuildCell2Cell();
1737
1739 pOpt.metisType = "KWAY";
1740 pOpt.metisUfactor = 30;
1741 pOpt.metisSeed = 42;
1742 pOpt.metisNcuts = 1;
1743 reader.MeshPartitionCell2Cell(pOpt);
1744 reader.PartitionReorderToMeshCell2Cell();
1745
1746 mO1->RecoverNode2CellAndNode2Bnd();
1747 mO1->RecoverCell2CellAndBnd2Cell();
1748 mO1->BuildGhostPrimary();
1749 mO1->AdjGlobal2LocalPrimary();
1750
1751 // Build O2 mesh
1752 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1753 mO2->BuildO2FromO1Elevation(*mO1);
1754 mO2->RecoverNode2CellAndNode2Bnd();
1755 mO2->RecoverCell2CellAndBnd2Cell();
1756 mO2->BuildGhostPrimary();
1757 mO2->AdjGlobal2LocalPrimary();
1758
1759 // Build O1 mesh via bisection
1760 mO2->AdjLocal2GlobalPrimary();
1761 auto mO1B = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1762 mO1B->BuildBisectO1FormO2(*mO2);
1763
1764 DNDS::index nCellOrig = mO1->NumCellGlobal();
1765 DNDS::index nCellElevated = mO2->NumCellGlobal();
1766 DNDS::index nCellBisected = mO1B->NumCellGlobal();
1767
1768 CHECK(nCellOrig == 100);
1769 CHECK(nCellElevated == 100);
1770 CHECK(nCellBisected == 400);
1771}
1772
1773// ===========================================================================
1774// Elevation + Boundary/Internal Smooth tests (NACA0012_H2, config 2)
1775// ===========================================================================
1776
1777/**
1778 * \brief Build an O2 mesh with full face interpolation, ready for smoothing.
1779 *
1780 * Sequence: read O1 -> partition -> ghost -> elevate -> rebuild topology
1781 * -> face interpolation -> N2CB ghost.
1782 */
1783static ssp<UnstructuredMesh> buildO2MeshWithFaces(
1784 const MeshConfig &cfg, ssp<UnstructuredMesh> &mO1Out)
1785{
1786 // Build O1 mesh
1787 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1788 UnstructuredMeshSerialRW reader(mO1, 0);
1789
1790 if (cfg.periodic)
1791 {
1792 tPoint zero{0, 0, 0};
1793 mO1->SetPeriodicGeometry(
1794 cfg.translation1, zero, zero,
1795 cfg.translation2, zero, zero,
1796 cfg.translation3, zero, zero);
1797 }
1798
1799 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1800 reader.Deduplicate1to1Periodic(1e-8);
1801 reader.BuildCell2Cell();
1802
1804 pOpt.metisType = "KWAY";
1805 pOpt.metisUfactor = 30;
1806 pOpt.metisSeed = 42;
1807 pOpt.metisNcuts = 1;
1808 reader.MeshPartitionCell2Cell(pOpt);
1809 reader.PartitionReorderToMeshCell2Cell();
1810
1811 mO1->RecoverNode2CellAndNode2Bnd();
1812 mO1->RecoverCell2CellAndBnd2Cell();
1813 mO1->BuildGhostPrimary();
1814 mO1->AdjGlobal2LocalPrimary();
1815
1816 // Build O2 mesh via elevation
1817 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1818 mO2->BuildO2FromO1Elevation(*mO1);
1819
1820 // Rebuild full topology for O2
1821 mO2->RecoverNode2CellAndNode2Bnd();
1822 mO2->RecoverCell2CellAndBnd2Cell();
1823 mO2->BuildGhostPrimary();
1824 mO2->AdjGlobal2LocalPrimary();
1825 mO2->AdjGlobal2LocalN2CB();
1826
1827 mO2->InterpolateFace();
1828 mO2->AssertOnFaces();
1829
1830 mO2->AdjLocal2GlobalN2CB();
1831 mO2->BuildGhostN2CB();
1832 mO2->AdjGlobal2LocalN2CB();
1833
1834 mO2->BuildCell2CellFace();
1835 mO2->AdjGlobal2LocalC2CFace();
1836
1837 mO1Out = mO1;
1838 return mO2;
1839}
1840
1841TEST_CASE("BoundarySmooth: NACA0012 moves boundary O2 nodes")
1842{
1843 // NACA0012_H2 is config 2 — non-periodic, curved airfoil boundary
1844 const auto &cfg = g_configs[2];
1845
1847 auto mO2 = buildO2MeshWithFaces(cfg, mO1);
1848
1849 // Snapshot pre-smooth coordinates
1850 std::vector<tPoint> preSmooth(mO2->NumNode());
1851 for (DNDS::index i = 0; i < mO2->NumNode(); i++)
1852 preSmooth[i] = mO2->coords[i];
1853
1854 // Run boundary smooth — all external BCs are smoothed
1855 mO2->ElevatedNodesGetBoundarySmooth(
1856 [](Geom::t_index bndId)
1857 { return Geom::FaceIDIsExternalBC(bndId); });
1858
1859 // Check that nTotalMoved > 0 (some O2 boundary nodes moved)
1860 CHECK(mO2->nTotalMoved > 0);
1861
1862 // Check that coordsElevDisp was initialized
1863 CHECK(mO2->coordsElevDisp.father != nullptr);
1864
1865 // All coordinates must be finite (no NaN/Inf)
1866 for (DNDS::index i = 0; i < mO2->NumNode(); i++)
1867 {
1868 for (int d = 0; d < 3; d++)
1869 CHECK(std::isfinite(mO2->coords[i](d)));
1870 }
1871}
1872
1873TEST_CASE("BoundarySmooth: no NaN in displacement field")
1874{
1875 const auto &cfg = g_configs[2]; // NACA0012_H2
1876
1878 auto mO2 = buildO2MeshWithFaces(cfg, mO1);
1879
1880 mO2->ElevatedNodesGetBoundarySmooth(
1881 [](Geom::t_index bndId)
1882 { return Geom::FaceIDIsExternalBC(bndId); });
1883
1884 // coordsElevDisp should have no NaN for moved nodes
1885 DNDS::index nMovedLocal = 0;
1886 for (DNDS::index i = 0; i < mO2->coordsElevDisp.father->Size(); i++)
1887 {
1888 auto disp = mO2->coordsElevDisp[i];
1889 // Entries with disp(0) != largeReal were set by the smooth
1890 if (disp(0) != UnInitReal && std::abs(disp(0)) < 1e30)
1891 {
1892 for (int d = 0; d < 3; d++)
1893 CHECK(std::isfinite(disp(d)));
1894 nMovedLocal++;
1895 }
1896 }
1897
1898 DNDS::index nMovedGlobal = 0;
1899 MPI_Allreduce(&nMovedLocal, &nMovedGlobal, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
1900 CHECK(nMovedGlobal > 0);
1901}
1902
1903TEST_CASE("InternalSmooth V2: coordinates remain finite after solve")
1904{
1905 const auto &cfg = g_configs[2]; // NACA0012_H2
1906
1908 auto mO2 = buildO2MeshWithFaces(cfg, mO1);
1909
1910 // Set up elevation info for fast convergence (small problem)
1911 mO2->elevationInfo.nIter = 10;
1912 mO2->elevationInfo.nSearch = 20;
1913 mO2->elevationInfo.RBFRadius = 5.0;
1914
1915 // Run boundary smooth
1916 mO2->ElevatedNodesGetBoundarySmooth(
1917 [](Geom::t_index bndId)
1918 { return Geom::FaceIDIsExternalBC(bndId); });
1919
1920 if (mO2->nTotalMoved == 0)
1921 return; // nothing to smooth
1922
1923 // Save pre-smooth coords
1924 std::vector<tPoint> preSmoothCoords(mO2->NumNode());
1925 for (DNDS::index i = 0; i < mO2->NumNode(); i++)
1926 preSmoothCoords[i] = mO2->coords[i];
1927
1928 // Run internal smooth V2
1929 mO2->ElevatedNodesSolveInternalSmoothV2();
1930
1931 // All coordinates must be finite
1932 for (DNDS::index i = 0; i < mO2->NumNode(); i++)
1933 for (int d = 0; d < 3; d++)
1934 CHECK(std::isfinite(mO2->coords[i](d)));
1935
1936 // At least some coordinates should have changed
1937 DNDS::index nChangedLocal = 0;
1938 for (DNDS::index i = 0; i < mO2->NumNode(); i++)
1939 if ((mO2->coords[i] - preSmoothCoords[i]).norm() > 1e-15)
1940 nChangedLocal++;
1941
1942 DNDS::index nChangedGlobal = 0;
1943 MPI_Allreduce(&nChangedLocal, &nChangedGlobal, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
1944 CHECK(nChangedGlobal > 0);
1945}
1946
1947TEST_CASE("Elevation: O2 coordinates have no NaN")
1948{
1949 // Test on all non-periodic configs (0, 2)
1950 for (int ci : {0, 2})
1951 {
1952 const auto &cfg = g_configs[ci];
1953 SUBCASE(cfg.name)
1954 {
1955 auto mO1 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1956 UnstructuredMeshSerialRW reader(mO1, 0);
1957
1958 reader.ReadFromCGNSSerial(meshPath(cfg.file));
1959 reader.Deduplicate1to1Periodic(1e-8);
1960 reader.BuildCell2Cell();
1961
1963 pOpt.metisType = "KWAY";
1964 pOpt.metisUfactor = 30;
1965 pOpt.metisSeed = 42;
1966 pOpt.metisNcuts = 1;
1967 reader.MeshPartitionCell2Cell(pOpt);
1968 reader.PartitionReorderToMeshCell2Cell();
1969
1970 mO1->RecoverNode2CellAndNode2Bnd();
1971 mO1->RecoverCell2CellAndBnd2Cell();
1972 mO1->BuildGhostPrimary();
1973 mO1->AdjGlobal2LocalPrimary();
1974
1975 auto mO2 = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1976 mO2->BuildO2FromO1Elevation(*mO1);
1977
1978 // All O2 coordinates must be finite
1979 for (DNDS::index i = 0; i < mO2->coords.father->Size(); i++)
1980 for (int d = 0; d < 3; d++)
1981 CHECK(std::isfinite(mO2->coords[i](d)));
1982 }
1983 }
1984}
1985
1986// ===========================================================================
1987// Multi-layer ghost tests (UniformSquare_10, config 0)
1988// ===========================================================================
1989
1990/// Build mesh through BuildGhostPrimary with a specific ghost layer count,
1991/// then run the full prepare pipeline.
1992static ssp<UnstructuredMesh> buildMeshWithGhostLayers(
1993 const MeshConfig &cfg, int nGhostLayers)
1994{
1995 auto mesh = std::make_shared<UnstructuredMesh>(g_mpi, cfg.dim);
1996 UnstructuredMeshSerialRW reader(mesh, 0);
1997
1998 if (cfg.periodic)
1999 {
2000 tPoint zero{0, 0, 0};
2001 mesh->SetPeriodicGeometry(
2002 cfg.translation1, zero, zero,
2003 cfg.translation2, zero, zero,
2004 cfg.translation3, zero, zero);
2005 }
2006
2007 reader.ReadFromCGNSSerial(meshPath(cfg.file));
2008 reader.Deduplicate1to1Periodic(1e-8);
2009 reader.BuildCell2Cell();
2010
2012 pOpt.metisType = "KWAY";
2013 pOpt.metisUfactor = 30;
2014 pOpt.metisSeed = 42;
2015 pOpt.metisNcuts = 1;
2016 reader.MeshPartitionCell2Cell(pOpt);
2017 reader.PartitionReorderToMeshCell2Cell();
2018
2019 mesh->RecoverNode2CellAndNode2Bnd();
2020 mesh->RecoverCell2CellAndBnd2Cell();
2021 mesh->BuildGhostPrimary(nGhostLayers);
2022 mesh->AdjGlobal2LocalPrimary();
2023 mesh->AdjGlobal2LocalN2CB();
2024
2025 mesh->InterpolateFace();
2026 mesh->AssertOnFaces();
2027
2028 mesh->AdjLocal2GlobalN2CB();
2029 mesh->BuildGhostN2CB();
2030 mesh->AdjGlobal2LocalN2CB();
2031
2032 return mesh;
2033}
2034
2035TEST_CASE("MultiLayerGhost: ghost counts increase with layer count")
2036{
2037 const auto &cfg = g_configs[0]; // UniformSquare_10
2038
2039 DNDS::index cellGhost[3], nodeGhost[3], bndGhost[3];
2040
2041 for (int nL = 1; nL <= 3; nL++)
2042 {
2043 auto m = buildMeshWithGhostLayers(cfg, nL);
2044 cellGhost[nL - 1] = m->NumCellGhost();
2045 nodeGhost[nL - 1] = m->NumNodeGhost();
2046 bndGhost[nL - 1] = m->NumBndGhost();
2047
2048 if (g_mpi.rank == 0)
2049 fmt::print(" nLayers={}: cellGhost={}, nodeGhost={}, bndGhost={}\n",
2050 nL, cellGhost[nL - 1], nodeGhost[nL - 1], bndGhost[nL - 1]);
2051 }
2052
2053 // Monotonically non-decreasing.
2054 for (int i = 0; i < 2; i++)
2055 {
2056 CHECK(cellGhost[i + 1] >= cellGhost[i]);
2057 CHECK(nodeGhost[i + 1] >= nodeGhost[i]);
2058 CHECK(bndGhost[i + 1] >= bndGhost[i]);
2059 }
2060
2061 // With np >= 2, cell and node ghosts must strictly increase.
2062 if (g_mpi.size >= 2)
2063 {
2064 CHECK(cellGhost[1] > cellGhost[0]);
2065 CHECK(cellGhost[2] > cellGhost[1]);
2066 CHECK(nodeGhost[1] > nodeGhost[0]);
2067 CHECK(nodeGhost[2] > nodeGhost[1]);
2068 }
2069}
2070
2071TEST_CASE("MultiLayerGhost: 2-layer cell2cell neighbors are all local")
2072{
2073 if (g_mpi.size < 2)
2074 return; // single-rank has no ghosts
2075
2076 const auto &cfg = g_configs[0]; // UniformSquare_10
2077 auto m = buildMeshWithGhostLayers(cfg, 2);
2078
2079 DNDS::index nCellOwned = m->NumCell();
2080 DNDS::index nCellProc = m->NumCellProc();
2081
2082 // Every owned cell's cell2cell neighbor's cell2cell neighbors
2083 // must be within [0, nCellProc).
2084 for (DNDS::index iCell = 0; iCell < nCellOwned; iCell++)
2085 {
2086 for (DNDS::rowsize j = 0; j < m->cell2cell.RowSize(iCell); j++)
2087 {
2088 DNDS::index iNeighbor = m->cell2cell(iCell, j);
2089 if (iNeighbor < 0)
2090 continue;
2091 REQUIRE(iNeighbor < nCellProc);
2092
2093 for (DNDS::rowsize k = 0; k < m->cell2cell.RowSize(iNeighbor); k++)
2094 {
2095 DNDS::index iNeighbor2 = m->cell2cell(iNeighbor, k);
2096 if (iNeighbor2 < 0)
2097 continue;
2098 CHECK(iNeighbor2 < nCellProc);
2099 }
2100 }
2101 }
2102}
2103
2104TEST_CASE("MultiLayerGhost: global cell count preserved across layers")
2105{
2106 const auto &cfg = g_configs[0]; // UniformSquare_10
2107
2108 for (int nL = 1; nL <= 3; nL++)
2109 {
2110 CAPTURE(nL);
2111 auto m = buildMeshWithGhostLayers(cfg, nL);
2112
2113 DNDS::index localCells = m->NumCell();
2114 DNDS::index globalCells = 0;
2115 MPI_Allreduce(&localCells, &globalCells, 1, DNDS_MPI_INDEX, MPI_SUM, g_mpi.comm);
2116 CHECK(globalCells == cfg.expectedCells);
2117 }
2118}
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:112
Eigen::Matrix< real, 3, 3 > m
int main()
Definition dummy.cpp:3
int32_t t_index
Definition Geometric.hpp:6
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
the host side operators are provided as implemented
const MPI_Datatype DNDS_MPI_INDEX
MPI datatype matching index (= MPI_INT64_T).
Definition MPI.hpp:106
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:181
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:114
DNDS_CONSTANT const real UnInitReal
Sentinel "not initialised" real value (NaN). Cheap to detect with std::isnan or IsUnInitReal; survive...
Definition Defines.hpp:179
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:112
std::shared_ptr< T > ssp
Shortened alias for std::shared_ptr used pervasively in DNDSR.
Definition Defines.hpp:143
Convenience bundle of a father, son, and attached ArrayTransformer.
DNDS_DEVICE_CALLABLE constexpr Element ObtainFace(t_index iFace) const
Definition Elements.hpp:257
DNDS_DEVICE_CALLABLE constexpr t_index GetNumVertices() const
Definition Elements.hpp:229
Lightweight bundle of an MPI communicator and the calling rank's coordinates.
Definition MPI.hpp:231
int size
Number of ranks in comm (-1 until initialised).
Definition MPI.hpp:237
int rank
This rank's 0-based index within comm (-1 until initialised).
Definition MPI.hpp:235
MPI_Comm comm
The underlying MPI communicator handle.
Definition MPI.hpp:233
void setWorld()
Initialise the object to MPI_COMM_WORLD. Requires MPI_Init to have run.
Definition MPI.hpp:258
DNDS::index expectedCells
Expected global cell count (-1 = skip check)
tPoint translation1
Periodic translation axis 1.
bool periodic
Has periodic boundaries.
int dim
Spatial dimension.
tPoint translation3
Periodic translation axis 3.
const char * file
Filename relative to data/mesh/.
tPoint translation2
Periodic translation axis 2.
const char * name
Human-readable label.
DNDS::index expectedBnds
Expected global bnd count (-1 = skip check)
tVec r(NCells)
auto adj
CHECK(result.size()==3)
auto result
REQUIRE(bool(result.parent2entityPbi.father))
DistributedHex3D mesh
#define FOR_EACH_MESH_CONFIG(body)
std::pair< DNDS::index, uint8_t > idx_pbi_t
input explicitMaps push_back(EntityReorderMap{EntityKind::Cell, cellPartition})
TEST_CASE("3D: VFV P2 HQM error < P1 on sinCos3D")