DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Mesh_Serial_ReadFromCGNS.cpp
Go to the documentation of this file.
1#include "Mesh.hpp"
2#include "Geom/CGNS.hpp"
3#include "Geom/OpenFOAMMesh.hpp"
4
5#include <array>
6#include <fmt/core.h>
7
8#include <cgnslib.h>
9
10#include <set>
11
12namespace DNDS::Geom
13{
14 // =================================================================
15 // File-local helpers for ReadFromCGNSSerial
16 // =================================================================
17 namespace
18 {
19 /**
20 * \brief Deduplicate shared nodes across CGNS zones via DFS on the
21 * zone connectivity graph.
22 *
23 * \param[in] ZoneCoords Per-zone coordinate arrays.
24 * \param[in] ZoneConnect Per-zone connectivity point lists (1-based).
25 * \param[in] ZoneConnectDonor Per-zone donor point lists (1-based).
26 * \param[in] ZoneConnectTargetIZone Per-zone target zone indices.
27 * \param[out] coordSerial Assembled coordinate array (resized and filled).
28 * \return NodeOld2New mapping from concatenated zone-local indices to
29 * assembled (deduplicated) global indices, plus ZoneNodeStarts.
30 */
31 std::pair<std::vector<DNDS::index>, std::vector<DNDS::index>>
32 AssembleZoneNodes(
33 const std::vector<tCoord> &ZoneCoords,
34 const std::vector<std::vector<std::vector<cgsize_t>>> &ZoneConnect,
35 const std::vector<std::vector<std::vector<cgsize_t>>> &ZoneConnectDonor,
36 const std::vector<std::vector<int>> &ZoneConnectTargetIZone,
37 tCoord &coordSerial,
38 int dim)
39 {
40 std::vector<DNDS::index> ZoneNodeSizes(ZoneCoords.size());
41 std::vector<DNDS::index> ZoneNodeStarts(ZoneCoords.size() + 1);
42 ZoneNodeStarts[0] = 0;
43 for (size_t i = 0; i < ZoneNodeSizes.size(); i++)
44 ZoneNodeStarts[i + 1] = ZoneNodeStarts[i] + (ZoneNodeSizes[i] = ZoneCoords[i]->Size());
45
46 std::vector<DNDS::index> NodeOld2New(ZoneNodeStarts.back(), -1);
47 DNDS::index cTop = 0;
48
49 // DFS over zone connectivity graph to deduplicate shared nodes
50 std::set<int> zonesLeft;
51 for (size_t iGZ = 0; iGZ < ZoneNodeSizes.size(); iGZ++)
52 zonesLeft.insert(static_cast<int>(iGZ));
53 std::vector<int> zonesFront;
54 zonesFront.push_back(0);
55 DNDS_assert(!ZoneNodeSizes.empty());
56
57 while (!zonesFront.empty())
58 {
59 int iGZ = zonesFront.back();
60 zonesFront.pop_back();
61 zonesLeft.erase(iGZ);
62 for (size_t iOther = 0; iOther < ZoneConnect.at(iGZ).size(); iOther++)
63 {
64 int iGZOther = ZoneConnectTargetIZone.at(iGZ).at(iOther);
65 if (zonesLeft.count(iGZOther))
66 {
67 zonesFront.push_back(iGZOther);
68 }
69 else
70 {
71 for (size_t iNode = 0; iNode < ZoneConnect.at(iGZ).at(iOther).size(); iNode++)
72 {
73 NodeOld2New.at(ZoneNodeStarts.at(iGZ) + ZoneConnect.at(iGZ).at(iOther).at(iNode) - 1) //! note ZoneConnect is 1-based
74 = NodeOld2New.at(ZoneNodeStarts.at(iGZOther) + ZoneConnectDonor.at(iGZ).at(iOther).at(iNode) - 1);
75 DNDS_assert(NodeOld2New.at(ZoneNodeStarts.at(iGZ) + ZoneConnect.at(iGZ).at(iOther).at(iNode) - 1) >= 0);
76 }
77 }
78 }
79 for (DNDS::index iNode = ZoneNodeStarts.at(iGZ); iNode < ZoneNodeStarts.at(iGZ + 1); iNode++)
80 {
81 if (NodeOld2New.at(iNode) < 0)
82 {
83 NodeOld2New.at(iNode) = cTop;
84 cTop++;
85 }
86 }
87 }
88 DNDS::log() << "CGNS === Assembled Zones have NNodes " << cTop << std::endl;
89 DNDS_assert_info(zonesLeft.empty(), "Did not reach all the zones, might missing zone-connectivities");
90
91 coordSerial->Resize(cTop);
92 for (DNDS::index i = 0; i < coordSerial->Size(); i++)
93 coordSerial->operator[](i).setConstant(DNDS::UnInitReal);
94
95 for (size_t iGZ = 0; iGZ < ZoneNodeSizes.size(); iGZ++)
96 {
97 for (DNDS::index iNode = ZoneNodeStarts.at(iGZ); iNode < ZoneNodeStarts.at(static_cast<int>(iGZ) + 1); iNode++)
98 {
99 auto coordC = (*ZoneCoords[iGZ])[iNode - ZoneNodeStarts.at(iGZ)];
100 real dist = ((*coordSerial)[NodeOld2New.at(iNode)] - coordC).norm();
101 if (!DNDS::IsUnInitReal((*coordSerial)[NodeOld2New.at(iNode)][0]))
102 DNDS_assert_info(dist < 1e-10,
103 "Not same points on the connection, distance is " + std::to_string(dist));
104 (*coordSerial)[NodeOld2New.at(iNode)] = coordC;
105 }
106 }
107
108 return {std::move(NodeOld2New), std::move(ZoneNodeStarts)};
109 }
110
111 /**
112 * \brief Separate zone elements into volume cells and boundary faces.
113 *
114 * Converts element node indices from zone-local to assembled global
115 * (via NodeOld2New), then fills cell2nodeSerial/cellElemInfoSerial and
116 * bnd2nodeSerial/bndElemInfoSerial.
117 */
118 void SeparateVolumeAndBoundaryElements(
119 const std::vector<tAdj> &ZoneElems,
120 const std::vector<tElemInfoArray> &ZoneElemInfos,
121 const std::vector<DNDS::index> &NodeOld2New,
122 const std::vector<DNDS::index> &ZoneNodeStarts,
123 tAdj &cell2nodeSerial,
124 tElemInfoArray &cellElemInfoSerial,
125 tAdj1 &cell2cellOrigSerial,
126 tAdj &bnd2nodeSerial,
127 tElemInfoArray &bndElemInfoSerial,
128 tAdj1 &bnd2bndOrigSerial,
129 tAdj1 &node2nodeOrigSerial,
130 tCoord &coordSerial,
131 int dim)
132 {
133 // Convert element node indices to assembled global and count
134 DNDS::index nVolElem = 0;
135 DNDS::index nBndElem = 0;
136 for (size_t iGZ = 0; iGZ < ZoneElems.size(); iGZ++)
137 {
138 for (DNDS::index iElem = 0; iElem < ZoneElems[iGZ]->Size(); iElem++)
139 {
140 for (DNDS::rowsize j = 0; j < ZoneElems[iGZ]->RowSize(iElem); j++)
141 ZoneElems[iGZ]->operator()(iElem, j) = NodeOld2New[ZoneNodeStarts[iGZ] + ZoneElems[iGZ]->operator()(iElem, j)];
142 auto Elem = Elem::Element{ZoneElemInfos[iGZ]->operator()(iElem, 0).getElemType()};
143 if (Elem.GetDim() == dim)
144 nVolElem++;
145 if (Elem.GetDim() == dim - 1 && !FaceIDIsTrueInternal(ZoneElemInfos[iGZ]->operator()(iElem, 0).zone))
146 nBndElem++;
147 }
148 }
149
150 cell2nodeSerial->Resize(nVolElem);
151 bnd2nodeSerial->Resize(nBndElem);
152 cellElemInfoSerial->Resize(nVolElem);
153 cell2cellOrigSerial->Resize(nVolElem);
154 bndElemInfoSerial->Resize(nBndElem);
155 node2nodeOrigSerial->Resize(coordSerial->Size());
156 bnd2bndOrigSerial->Resize(nBndElem);
157 nVolElem = 0;
158 nBndElem = 0;
159 for (size_t iGZ = 0; iGZ < ZoneElems.size(); iGZ++)
160 {
161 for (DNDS::index iElem = 0; iElem < ZoneElems[iGZ]->Size(); iElem++)
162 {
163 auto Elem = Elem::Element{ZoneElemInfos[iGZ]->operator()(iElem, 0).getElemType()};
164 if (Elem.GetDim() == dim)
165 {
166 cell2nodeSerial->ResizeRow(nVolElem, ZoneElems[iGZ]->RowSize(iElem));
167 for (DNDS::rowsize j = 0; j < ZoneElems[iGZ]->RowSize(iElem); j++)
168 cell2nodeSerial->operator()(nVolElem, j) = ZoneElems[iGZ]->operator()(iElem, j);
169 cellElemInfoSerial->operator()(nVolElem, 0) = ZoneElemInfos[iGZ]->operator()(iElem, 0);
170 nVolElem++;
171 }
172 if (Elem.GetDim() == dim - 1 && !FaceIDIsTrueInternal(ZoneElemInfos[iGZ]->operator()(iElem, 0).zone))
173 {
174 bnd2nodeSerial->ResizeRow(nBndElem, ZoneElems[iGZ]->RowSize(iElem));
175 for (DNDS::rowsize j = 0; j < ZoneElems[iGZ]->RowSize(iElem); j++)
176 bnd2nodeSerial->operator()(nBndElem, j) = ZoneElems[iGZ]->operator()(iElem, j);
177 bndElemInfoSerial->operator()(nBndElem, 0) = ZoneElemInfos[iGZ]->operator()(iElem, 0);
178 nBndElem++;
179 }
180 }
181 }
182 log() << "CGNS === Vol Elem [ " << nVolElem << " ]"
183 << ", "
184 << " Bnd Elem [ " << nBndElem << " ]" << std::endl;
185
186 bnd2nodeSerial->Compress();
187 cell2nodeSerial->Compress();
188 }
189
190 /**
191 * \brief Build bnd2cellSerial: for each boundary face, find its parent
192 * volume cell by vertex-set inclusion.
193 *
194 * Uses node-to-boundary index for linear complexity.
195 */
196 void BuildBnd2CellSerial(
197 tAdj2 &bnd2cellSerial,
198 const tAdj &bnd2nodeSerial,
199 const tAdj &cell2nodeSerial,
200 const tElemInfoArray &cellElemInfoSerial,
201 const tElemInfoArray &bndElemInfoSerial,
202 const tCoord &coordSerial)
203 {
204 bnd2cellSerial->Resize(bnd2nodeSerial->Size());
205 for (DNDS::index iB = 0; iB < bnd2cellSerial->Size(); iB++)
206 (*bnd2cellSerial)[iB][0] = (*bnd2cellSerial)[iB][1] = DNDS::UnInitIndex;
207
208 // Build node-to-boundary index for linear complexity
209 std::vector<std::vector<DNDS::index>> node2bnd(coordSerial->Size());
210 std::vector<DNDS::rowsize> node2bndSiz(coordSerial->Size(), 0);
211 for (DNDS::index iBFace = 0; iBFace < bnd2nodeSerial->Size(); iBFace++)
212 for (DNDS::rowsize iN = 0; iN < Elem::Element{(*bndElemInfoSerial)(iBFace, 0).getElemType()}.GetNumVertices(); iN++)
213 node2bndSiz[(*bnd2nodeSerial)(iBFace, iN)]++;
214 for (DNDS::index iNode = 0; iNode < static_cast<DNDS::index>(node2bnd.size()); iNode++)
215 node2bnd[iNode].reserve(node2bndSiz[iNode]);
216 for (DNDS::index iBFace = 0; iBFace < bnd2nodeSerial->Size(); iBFace++)
217 for (DNDS::rowsize iN = 0; iN < Elem::Element{(*bndElemInfoSerial)(iBFace, 0).getElemType()}.GetNumVertices(); iN++)
218 node2bnd[(*bnd2nodeSerial)(iBFace, iN)].push_back(iBFace);
219
220 // Search each cell for matching boundary faces
221 for (DNDS::index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
222 {
223 auto CellElem = Elem::Element{(*cellElemInfoSerial)[iCell]->getElemType()};
224 std::vector<DNDS::index> cell2nodeRow{(*cell2nodeSerial)[iCell].begin(), (*cell2nodeSerial)[iCell].begin() + CellElem.GetNumVertices()};
225 std::sort(cell2nodeRow.begin(), cell2nodeRow.end());
226 for (auto iNode : cell2nodeRow)
227 {
228 for (auto iB : node2bnd[iNode])
229 {
230 auto BndElem = Elem::Element{(*bndElemInfoSerial)(iB, 0).getElemType()};
231 std::vector<DNDS::index> bnd2nodeRow{(*bnd2nodeSerial)[iB].begin(), (*bnd2nodeSerial)[iB].begin() + BndElem.GetNumVertices()};
232 std::sort(bnd2nodeRow.begin(), bnd2nodeRow.end());
233 if (std::includes(cell2nodeRow.begin(), cell2nodeRow.end(), bnd2nodeRow.begin(), bnd2nodeRow.end()))
234 {
235 DNDS_assert_info(
236 (*bnd2cellSerial)[iB][0] == DNDS::UnInitIndex || (*bnd2cellSerial)[iB][0] == iCell, "bnd2cell not untouched!");
237 (*bnd2cellSerial)[iB][0] = iCell;
238 }
239 }
240 }
241 }
242 for (DNDS::index iB = 0; iB < bnd2cellSerial->Size(); iB++)
243 DNDS_assert((*bnd2cellSerial)[iB][0] != DNDS::UnInitIndex);
244 }
245 } // anonymous namespace
246 void UnstructuredMeshSerialRW::
247 ReadFromCGNSSerial(const std::string &fName, const t_FBCName_2_ID &FBCName_2_ID)
248 {
249 mode = SerialReadAndDistribute;
250 this->dataIsSerialIn = true;
251
252 int cgErr = CG_OK;
253
254 cell2nodeSerial = make_ssp<decltype(cell2nodeSerial)::element_type>(ObjName{"ReadFromCGNSSerial::cell2nodeSerial"}, mesh->getMPI());
255 bnd2nodeSerial = make_ssp<decltype(bnd2nodeSerial)::element_type>(ObjName{"ReadFromCGNSSerial::bnd2nodeSerial"}, mesh->getMPI());
256 coordSerial = make_ssp<decltype(coordSerial)::element_type>(ObjName{"ReadFromCGNSSerial::coordSerial"}, mesh->getMPI());
257 cellElemInfoSerial = make_ssp<decltype(cellElemInfoSerial)::element_type>(ObjName{"ReadFromCGNSSerial::cellElemInfoSerial"}, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
258 bndElemInfoSerial = make_ssp<decltype(bndElemInfoSerial)::element_type>(ObjName{"ReadFromCGNSSerial::bndElemInfoSerial"}, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
259 bnd2cellSerial = make_ssp<decltype(bnd2cellSerial)::element_type>(ObjName{"ReadFromCGNSSerial::bnd2cellSerial"}, mesh->getMPI());
260 cell2cellOrigSerial = make_ssp<decltype(cell2cellOrigSerial)::element_type>(ObjName{"ReadFromCGNSSerial::cell2cellOrigSerial"}, mesh->getMPI());
261 node2nodeOrigSerial = make_ssp<decltype(node2nodeOrigSerial)::element_type>(ObjName{"ReadFromCGNSSerial::node2nodeOrigSerial"}, mesh->getMPI());
262 bnd2bndOrigSerial = make_ssp<decltype(bnd2bndOrigSerial)::element_type>(ObjName{"ReadFromCGNSSerial::bnd2bndOrigSerial"}, mesh->getMPI());
263
264 if (mRank != mesh->getMPI().rank) //! parallel done!!! now serial!!!
265 return;
266
267 int cgns_file = -1;
268 if (cg_open(fName.c_str(), CG_MODE_READ, &cgns_file))
269 {
270 std::cerr << fmt::format("cgns file cannot open: [{}]", fName) << std::endl;
271 cg_error_exit();
272 }
273 int n_bases = -1;
274 if (cg_nbases(cgns_file, &n_bases))
275 cg_error_exit();
276 DNDS::log() << "CGNS === N bases: " << n_bases << std::endl;
277 DNDS_assert(n_bases > 0);
278 DNDS_assert(n_bases == 1);
279 std::vector<std::pair<int, int>> Base_Zone;
280 std::vector<std::string> BaseNames;
281 std::vector<std::string> ZoneNames;
282 // std::vector<std::string> ZoneFamilyNames;
283 std::vector<std::array<cgsize_t, 9>> ZoneSizes;
284 std::vector<tCoord> ZoneCoords; //! purely serial
285 std::vector<tAdj> ZoneElems;
286 std::vector<tElemInfoArray> ZoneElemInfos;
287 std::vector<std::vector<std::vector<cgsize_t>>> ZoneConnect;
288 std::vector<std::vector<std::vector<cgsize_t>>> ZoneConnectDonor;
289 std::vector<std::vector<int>> ZoneConnectTargetIZone;
290 /***************************************************************************/
291 // TRAVERSE 1
292 for (int iBase = 1; iBase <= n_bases; iBase++)
293 {
294
295 int cgns_base = -1;
296 std::array<char, 48> basename{};
297 int celldim = -1;
298 int physdim = -1;
299 int nzones = -1;
300
301 if (cg_base_read(cgns_file, iBase, basename.data(), &celldim, &physdim))
302 cg_error_exit();
303 DNDS_assert_info(celldim == mesh->dim, "CGNS file need to be with correct dim");
304 if (cg_nzones(cgns_file, iBase, &nzones))
305 cg_error_exit();
306 for (int iZone = 1; iZone <= nzones; iZone++)
307 {
308 std::array<char, 48> zonename{};
309 std::array<cgsize_t, 9> size{0, 0, 0, 0, 0, 0, 0, 0, 0};
310 if (cg_zone_read(cgns_file, iBase, iZone, zonename.data(), size.data()))
311 cg_error_exit();
312 ZoneType_t zoneType;
313 if (cg_zone_type(cgns_file, iBase, iZone, &zoneType))
314 cg_error_exit();
315 DNDS_assert(zoneType == Unstructured); //! only supports unstructured
316 Base_Zone.emplace_back(iBase, iZone);
317 BaseNames.emplace_back(basename.data());
318 ZoneNames.emplace_back(zonename.data());
319
320 // ! family name used for Volume condition not used now
321 // if (cg_goto(cgns_file, iBase, "Zone_t", iZone, ""))
322 // cg_error_exit();
323 // char famname[48]{0};
324 // if (cg_famname_read(famname))
325 // cg_error_exit();
326 // ZoneFamilyNames.push_back(famname); //* family name used for Volume condition
327 //
328
329 DNDS::index nNodes = size[0];
330 DNDS::index nVols = size[1];
331 DNDS::index nBVertex = size[2];
332
333 //*** read coords
334 {
335 int ncoords;
336 DNDS_CGNS_CALL_EXIT(cg_ncoords(cgns_file, iBase, iZone, &ncoords));
337 // DNDS_assert(ncoords == mesh->dim); //!Even if has Z coord, we only use X,Y in 2d
338
339 ZoneCoords.emplace_back(std::make_shared<decltype(ZoneCoords)::value_type::element_type>());
340 ZoneCoords.back()->Resize(nNodes);
341 std::vector<double> coordsRead(nNodes);
342 std::array<cgsize_t, 3> RMin{1, 0, 0};
343 std::array<cgsize_t, 3> RMax{nNodes, 0, 0};
344 //* READ X
345 DNDS_CGNS_CALL_EXIT(cg_coord_read(cgns_file, iBase, iZone, "CoordinateX", RealDouble, RMin.data(), RMax.data(), coordsRead.data()));
346 for (DNDS::index i = 0; i < nNodes; i++)
347 ZoneCoords.back()->operator[](i)[0] = coordsRead.at(i);
348 //* READ Y
349 DNDS_CGNS_CALL_EXIT(cg_coord_read(cgns_file, iBase, iZone, "CoordinateY", RealDouble, RMin.data(), RMax.data(), coordsRead.data()));
350 for (DNDS::index i = 0; i < nNodes; i++)
351 ZoneCoords.back()->operator[](i)[1] = coordsRead.at(i);
352 //* READ Z
353 if (mesh->dim == 3)
354 {
355 DNDS_CGNS_CALL_EXIT(cg_coord_read(cgns_file, iBase, iZone, "CoordinateZ", RealDouble, RMin.data(), RMax.data(), coordsRead.data()));
356 }
357 else
358 for (auto &i : coordsRead)
359 i = 0;
360 for (DNDS::index i = 0; i < nNodes; i++)
361 ZoneCoords.back()->operator[](i)[2] = coordsRead.at(i);
362
363 DNDS::log() << "CGNS === Zone " << iZone << " Coords Reading Done" << std::endl;
364 }
365 // for (DNDS::index i = 0; i < ZoneCoords.back()->Size(); i++)
366 // std::cout << (*ZoneCoords.back())[i].transpose() << std::endl;
367
368 //*** read Elems
369 {
370 int nSections;
371 DNDS_CGNS_CALL_EXIT(cg_nsections(cgns_file, iBase, iZone, &nSections));
372 DNDS_assert(nSections >= 1);
373 // [[maybe_unused]]: These track the current section's start/end for potential
374 // validation of section contiguity. The original assertion (cend == start - 1) was
375 // removed because CGNS sections may have gaps. Preserved for documentation and
376 // in case stricter validation is needed for specific file formats.
377 [[maybe_unused]] cgsize_t cstart = 0;
378 [[maybe_unused]] cgsize_t cend = 0;
379 cgsize_t maxend = 0;
380 //*Total size
381 for (int iSection = 1; iSection <= nSections; iSection++)
382 {
383 std::array<char, 48> sectionName{};
384 ElementType_t etype;
385 cgsize_t start;
386 cgsize_t end;
387 int nBnd{0}, parentFlag{0};
388 DNDS_CGNS_CALL_EXIT(cg_section_read(cgns_file, iBase, iZone, iSection, sectionName.data(), &etype, &start, &end, &nBnd, &parentFlag));
389 // cgsize_t elemDataSize{0};
390 // DNDS_CGNS_CALL_EXIT(cg_ElementDataSize(cgns_file, iBase, iZone, iSection, &elemDataSize));
391 // DNDS_assert(cend == start - 1); //? testing//!not valid!
392 cstart = start;
393 cend = end;
394 maxend = std::max(end, maxend);
395 }
396 ZoneElems.emplace_back(std::make_shared<decltype(ZoneElems)::value_type::element_type>());
397 ZoneElems.back()->Resize(maxend);
398 ZoneElemInfos.emplace_back(std::make_shared<decltype(ZoneElemInfos)::value_type::element_type>());
399 ZoneElemInfos.back()->Resize(maxend);
400 //*Resize row
401 for (int iSection = 1; iSection <= nSections; iSection++)
402 {
403 std::array<char, 48> sectionName{};
404 ElementType_t etype;
405 cgsize_t start;
406 cgsize_t end;
407 int nBnd{0}, parentFlag{0};
408 DNDS_CGNS_CALL_EXIT(cg_section_read(cgns_file, iBase, iZone, iSection, sectionName.data(), &etype, &start, &end, &nBnd, &parentFlag));
409 cgsize_t elemDataSize{0};
410 DNDS_CGNS_CALL_EXIT(cg_ElementDataSize(cgns_file, iBase, iZone, iSection, &elemDataSize));
411 std::vector<cgsize_t> elemsRead(elemDataSize);
412
413 int nElemSec = end - start + 1;
414 if (etype == MIXED)
415 {
416 std::vector<cgsize_t> elemStarts(nElemSec + 1); //* note size
417 DNDS_CGNS_CALL_EXIT(cg_poly_elements_read(cgns_file, iBase, iZone, iSection, elemsRead.data(), elemStarts.data(), nullptr));
418 for (cgsize_t i = 0; i < nElemSec; i++)
419 {
420 auto c_etype = static_cast<ElementType_t>(elemsRead.at(elemStarts[i]));
421 if (_getElemTypeFromCGNSType(c_etype) == Elem::UnknownElem)
422 {
423 DNDS::log() << "Error ETYPE " << std::to_string(c_etype) << std::endl;
424 DNDS_assert_info(false, "Unsupported Element! ");
425 }
426 Elem::ElemType ct = _getElemTypeFromCGNSType(c_etype);
427 DNDS_assert_info(Elem::Element{ct}.GetNumNodes() + 1 == elemStarts[i + 1] - elemStarts[i],
428 "Element Node Number Mismatch!");
429 ZoneElems.back()->ResizeRow(start - 1 + i, Elem::Element{ct}.GetNumNodes());
430 for (t_index iNode = 0; iNode < Elem::Element{ct}.GetNumNodes(); iNode++)
431 ZoneElems.back()->operator()(start - 1 + i, iNode) =
432 elemsRead.at(elemStarts[i] + 1 + iNode) - 1; //! convert to 0 based; pointing to zonal index
433 ZoneElemInfos.back()->operator()(start - 1 + i, 0).setElemType(ct);
434 ZoneElemInfos.back()->operator()(start - 1 + i, 0).zone = Geom::BC_ID_INTERNAL; //! initialized as inner,need new way of doing vol condition
435 }
436 /// @todo //TODO: TEST with actual data (MIXED TYPE) !!!!!!
437 }
438 else if (_getElemTypeFromCGNSType(etype) != Elem::UnknownElem)
439 {
440 Elem::ElemType ct = _getElemTypeFromCGNSType(etype);
441 DNDS_CGNS_CALL_EXIT(cg_elements_read(cgns_file, iBase, iZone, iSection, elemsRead.data(), nullptr));
442 DNDS_assert(elemDataSize / Elem::Element{ct}.GetNumNodes() == nElemSec);
443 for (cgsize_t i = 0; i < nElemSec; i++)
444 {
445
446 ZoneElems.back()->ResizeRow(start - 1 + i, Elem::Element{ct}.GetNumNodes());
447 for (t_index iNode = 0; iNode < Elem::Element{ct}.GetNumNodes(); iNode++)
448 ZoneElems.back()->operator()(start - 1 + i, iNode) =
449 elemsRead.at(Elem::Element{ct}.GetNumNodes() * i + iNode) - 1; //! convert to 0 based; pointing to zonal index
450 ZoneElemInfos.back()->operator()(start - 1 + i, 0).setElemType(ct);
451 ZoneElemInfos.back()->operator()(start - 1 + i, 0).zone = Geom::BC_ID_INTERNAL; //! initialized as inner,need new way of doing vol condition
452 }
453 }
454 else
455 {
456 DNDS::log() << "Error ETYPE " << std::to_string(etype) << std::endl;
457 DNDS_assert_info(false, "Unsupported Element! ");
458 }
459 }
460 DNDS::log() << "CGNS === Zone " << iZone << " Elems Reading Done" << std::endl;
461 }
462
463 //* read BCs
464 {
465 int nBC;
466 DNDS_CGNS_CALL_EXIT(cg_nbocos(cgns_file, iBase, iZone, &nBC));
467 for (int iBC = 1; iBC <= nBC; iBC++)
468 {
469 std::array<char, 48> boconame{};
470 PointSetType_t pType;
471 cgsize_t nPts, normalListSize;
472 BCType_t bcType;
473 std::array<int, 3> NormalIndex{};
474 DataType_t normalDataType;
475 int nDataset;
476 GridLocation_t gloc;
477 DNDS_CGNS_CALL_EXIT(cg_boco_info(cgns_file, iBase, iZone, iBC, boconame.data(), &bcType, &pType, &nPts, NormalIndex.data(), &normalListSize, &normalDataType, &nDataset));
478 DNDS_CGNS_CALL_EXIT(cg_boco_gridlocation_read(cgns_file, iBase, iZone, iBC, &gloc));
479 DNDS_assert_info(pType == PointRange || pType == PointList, "Only PointRange / PointList supported in BC!");
480 if (mesh->dim == 2)
481 DNDS_assert(gloc == EdgeCenter);
482 if (mesh->dim == 3)
483 DNDS_assert(gloc == FaceCenter);
484 std::vector<cgsize_t> pts(nPts);
485 std::vector<double> normalBuf(normalListSize); // should have checked normalDataType, but it is not used here so not checked
486 DNDS_CGNS_CALL_EXIT(cg_boco_read(cgns_file, iBase, iZone, iBC, pts.data(), normalBuf.data()));
487 DNDS_assert(pts[0] >= 1 && pts[1] <= ZoneElems.back()->Size());
488
489 t_index BCCode = FBCName_2_ID(std::string(boconame.data()));
490 if (BCCode == BC_ID_NULL)
491 {
492 DNDS_assert_info(false, fmt::format("BC NAME [{}] NOT FOUND IN DATABASE", boconame.data()));
493 }
494 if (pType == PointRange)
495 {
496 for (DNDS::index i = pts[0] - 1; i < pts[1]; i++)
497 ZoneElemInfos.back()->operator()(i, 0).zone = BCCode; //! setting BC code
498 }
499 else if (pType == PointList)
500 {
501 for (auto i : pts)
502 ZoneElemInfos.back()->operator()(i - 1, 0).zone = BCCode; //* note that pts is 1-based
503 }
504 else
505 DNDS_assert(false);
506 }
507
508 DNDS::log() << "CGNS === Zone " << iZone << " BCs Reading Done" << std::endl;
509 }
510 }
511 }
512 /***************************************************************************/
513
514 /***************************************************************************/
515 // TRAVERSE 2
516 for (int iBase = 1; iBase <= n_bases; iBase++)
517 {
518
519 int cgns_base = -1;
520 std::array<char, 48> basename{};
521 int celldim = -1;
522 int physdim = -1;
523 int nzones = -1;
524
525 if (cg_base_read(cgns_file, iBase, basename.data(), &celldim, &physdim))
526 cg_error_exit();
527 DNDS_assert(celldim == mesh->dim);
528 if (cg_nzones(cgns_file, iBase, &nzones))
529 cg_error_exit();
530 for (int iZone = 1; iZone <= nzones; iZone++)
531 {
532 std::array<char, 48> zonename{};
533 std::array<cgsize_t, 9> size{0, 0, 0, 0, 0, 0, 0, 0, 0};
534 if (cg_zone_read(cgns_file, iBase, iZone, zonename.data(), size.data()))
535 cg_error_exit();
536 ZoneType_t zoneType;
537 if (cg_zone_type(cgns_file, iBase, iZone, &zoneType))
538 cg_error_exit();
539 DNDS_assert(zoneType == Unstructured); //! only supports unstructured
540
541 DNDS::index nNodes = size[0];
542 DNDS::index nVols = size[1];
543 DNDS::index nBVertex = size[2];
544
545 //*** read Connectivity
546 {
547 ZoneConnect.emplace_back();
548 ZoneConnectDonor.emplace_back();
549 ZoneConnectTargetIZone.emplace_back();
550 int nConns;
551 DNDS_CGNS_CALL_EXIT(cg_nconns(cgns_file, iBase, iZone, &nConns));
552 for (int iConn = 1; iConn <= nConns; iConn++)
553 {
554 std::array<char, 48> connName{}, donorName{};
555 GridLocation_t gLoc;
556 GridConnectivityType_t connType;
557 PointSetType_t ptType, ptType_donor;
558 cgsize_t npts, npts_donor;
559 ZoneType_t donorZT;
560 DataType_t donorDT;
561
562 DNDS_CGNS_CALL_EXIT(cg_conn_info(cgns_file, iBase, iZone, iConn, connName.data(), &gLoc, &connType, &ptType, &npts,
563 donorName.data(), &donorZT, &ptType_donor, &donorDT, &npts_donor));
564
565 DNDS_assert_info(connType == Abutting1to1, "Only support Abutting1to1 in connection!");
566 DNDS_assert_info(ptType == PointList, "Only Supports PointList in connection!");
567 DNDS_assert_info(ptType_donor == PointListDonor, "Only Supports PointListDonor in connection!");
568 DNDS_assert_info(donorZT == Unstructured, "Only Supports Unstructured in connection!");
569 DNDS_assert_info(gLoc == Vertex, "Only Supports Vertex in connection!");
570 DNDS_assert_info(npts_donor == npts, "Only Supports npts_donor == npts in connection!");
571 // std::vector<cgsize_t> ptSet(npts);
572 // std::vector<cgsize_t> ptSet_donor(npts);
573 ZoneConnectDonor.back().emplace_back(npts);
574 ZoneConnect.back().emplace_back(npts);
575 DNDS_CGNS_CALL_EXIT(cg_conn_read(cgns_file, iBase, iZone, iConn, ZoneConnect.back().back().data(), donorDT, ZoneConnectDonor.back().back().data()));
576 int iGZFound = -1;
577 for (size_t iGZ = 0; iGZ < ZoneNames.size(); iGZ++) // find the donor
578 {
579 if (Base_Zone.at(iGZ).first == iBase)
580 {
581 if (ZoneNames.at(iGZ) == donorName.data())
582 {
583 iGZFound = iGZ;
584 }
585 }
586 }
587 DNDS_assert(iGZFound >= 0);
588 ZoneConnectTargetIZone.back().push_back(iGZFound);
589 DNDS::log() << "CGNS === Connection at Zone-Zone: " << iZone << " - " << iGZFound + 1 << std::endl;
590 }
591 }
592 }
593 }
594
595 /***************************************************************************/
596
597 /***************************************************************************/
598 // ASSEMBLE: deduplicate shared nodes across zones
599 auto [NodeOld2New, ZoneNodeStarts] = AssembleZoneNodes(
600 ZoneCoords, ZoneConnect, ZoneConnectDonor, ZoneConnectTargetIZone,
601 coordSerial, mesh->dim);
602
603 // Separate volume cells and boundary faces
604 SeparateVolumeAndBoundaryElements(
605 ZoneElems, ZoneElemInfos, NodeOld2New, ZoneNodeStarts,
606 cell2nodeSerial, cellElemInfoSerial, cell2cellOrigSerial,
607 bnd2nodeSerial, bndElemInfoSerial, bnd2bndOrigSerial,
608 node2nodeOrigSerial, coordSerial, mesh->dim);
609 /***************************************************************************/
610
611 /***************************************************************************/
612 // Build bnd2cell inverse mapping
613 BuildBnd2CellSerial(
614 bnd2cellSerial, bnd2nodeSerial, cell2nodeSerial,
615 cellElemInfoSerial, bndElemInfoSerial, coordSerial);
616
617 // fill in original indices
618 for (DNDS::index iCell = 0; iCell < cell2cellOrigSerial->Size(); iCell++)
619 cell2cellOrigSerial->operator()(iCell, 0) = iCell;
620 for (DNDS::index iNode = 0; iNode < node2nodeOrigSerial->Size(); iNode++)
621 node2nodeOrigSerial->operator()(iNode, 0) = iNode;
622 for (DNDS::index iBnd = 0; iBnd < bnd2bndOrigSerial->Size(); iBnd++)
623 bnd2bndOrigSerial->operator()(iBnd, 0) = iBnd;
624
625 cg_close(cgns_file);
626
627 log() << "CGNS === Serial Read Done" << std::endl;
628 // Memory with DM240-120 here: 18G ; after deconstruction done: 7.5G
629 }
630
631 void UnstructuredMeshSerialRW::
632 ReadFromOpenFOAMAndConvertSerial(const std::string &fName, const std::map<std::string, std::string> &nameMapping, const t_FBCName_2_ID &FBCName_2_ID)
633 {
634 mode = SerialReadAndDistribute;
635 this->dataIsSerialIn = true;
636
637 cell2nodeSerial = make_ssp<decltype(cell2nodeSerial)::element_type>(ObjName{"ReadFromOpenFOAMAndConvertSerial::cell2nodeSerial"}, mesh->getMPI());
638 bnd2nodeSerial = make_ssp<decltype(bnd2nodeSerial)::element_type>(ObjName{"ReadFromOpenFOAMAndConvertSerial::bnd2nodeSerial"}, mesh->getMPI());
639 coordSerial = make_ssp<decltype(coordSerial)::element_type>(ObjName{"ReadFromOpenFOAMAndConvertSerial::coordSerial"}, mesh->getMPI());
640 cellElemInfoSerial = make_ssp<decltype(cellElemInfoSerial)::element_type>(ObjName{"ReadFromOpenFOAMAndConvertSerial::cellElemInfoSerial"}, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
641 bndElemInfoSerial = make_ssp<decltype(bndElemInfoSerial)::element_type>(ObjName{"ReadFromOpenFOAMAndConvertSerial::bndElemInfoSerial"}, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
642 bnd2cellSerial = make_ssp<decltype(bnd2cellSerial)::element_type>(ObjName{"ReadFromOpenFOAMAndConvertSerial::bnd2cellSerial"}, mesh->getMPI());
643
644 if (mRank != mesh->getMPI().rank) //! parallel done!!! now serial!!!
645 return;
646
647 std::filesystem::path ofFilePath(fName);
648 std::ifstream pointsIFS(ofFilePath / "points");
649 std::ifstream facesIFS(ofFilePath / "faces");
650 std::ifstream ownerIFS(ofFilePath / "owner");
651 std::ifstream neighbourIFS(ofFilePath / "neighbour");
652 std::ifstream boundaryIFS(ofFilePath / "boundary");
653
654 DNDS_assert_info(pointsIFS.is_open(), "points file not found!");
655 DNDS_assert_info(facesIFS.is_open(), "faces file not found!");
656 DNDS_assert_info(ownerIFS.is_open(), "owner file not found!");
657 DNDS_assert_info(neighbourIFS.is_open(), "neighbour file not found!");
658 DNDS_assert_info(boundaryIFS.is_open(), "boundary file not found!");
659
660 OpenFOAM::OpenFOAMReader ofReader;
661 ofReader.ReadPoints(pointsIFS);
662 ofReader.ReadFaces(facesIFS);
663 ofReader.ReadOwner(ownerIFS);
664 ofReader.ReadNeighbour(neighbourIFS);
665 ofReader.ReadBoundary(boundaryIFS);
666
667 OpenFOAM::OpenFOAMConverter ofConverter;
668 ofConverter.BuildFaceElemInfo(ofReader);
669 ofConverter.BuildCell2Face(ofReader);
670 ofConverter.BuildCell2Node(ofReader);
671
672 coordSerial->Resize(ofReader.points.size());
673 for (index iN = 0; iN < coordSerial->Size(); iN++)
674 (*coordSerial)[iN] = ofReader.points.at(iN);
675 log() << "OpenFOAM === got num node: " << coordSerial->Size() << std::endl;
676
677 cell2nodeSerial->Resize(ofConverter.cell2node.size());
678 cellElemInfoSerial->Resize(ofConverter.cell2node.size());
679 for (index iC = 0; iC < cell2nodeSerial->Size(); iC++)
680 {
681 cell2nodeSerial->ResizeRow(iC, ofConverter.cell2node[iC].size());
682 for (size_t iN = 0; iN < ofConverter.cell2node[iC].size(); iN++)
683 (*cell2nodeSerial)[iC][iN] = ofConverter.cell2node[iC][iN];
684 cellElemInfoSerial->operator()(iC, 0) = ofConverter.cellElemInfo.at(iC);
685 }
686
687 log() << "OpenFOAM === got num cell: " << cell2nodeSerial->Size() << std::endl;
688
689 index nBnd = ofReader.owner.size() - ofReader.neighbour.size();
690 bnd2nodeSerial->Resize(nBnd);
691 bndElemInfoSerial->Resize(nBnd);
692 bnd2cellSerial->Resize(nBnd);
693 for (index iBnd = 0; iBnd < nBnd; iBnd++)
694 {
695 index iFaceOF = iBnd + ofReader.neighbour.size();
696 bnd2nodeSerial->ResizeRow(iBnd, ofReader.faces[iFaceOF].size());
697 for (size_t ib2c = 0; ib2c < ofReader.faces[iFaceOF].size(); ib2c++)
698 (*bnd2nodeSerial)[iBnd][ib2c] = ofReader.faces[iFaceOF][ib2c];
699 bnd2cellSerial->operator()(iBnd, 0) = ofReader.owner.at(iFaceOF);
700 bnd2cellSerial->operator()(iBnd, 1) = UnInitIndex;
701 bndElemInfoSerial->operator()(iBnd, 0) = ofConverter.faceElemInfo.at(iFaceOF);
702 }
703
704 log() << "nameMapping size: " << nameMapping.size() << std::endl;
705
706 for (auto &bc : ofReader.boundaryConditions)
707 {
708 auto boconame = bc.first;
709 if (nameMapping.count(boconame))
710 boconame = nameMapping.at(boconame);
711 t_index BCCode = FBCName_2_ID(std::string(boconame));
712 if (BCCode == BC_ID_NULL)
713 {
714 DNDS_assert_info(false, fmt::format("BC NAME [{}] NOT FOUND IN DATABASE", boconame));
715 }
716 for (index iBndC = 0; iBndC < bc.second.nFaces; iBndC++)
717 bndElemInfoSerial->operator()(iBndC + bc.second.startFace - ofReader.neighbour.size(), 0).zone = BCCode;
718 }
719 for (index iBnd = 0; iBnd < nBnd; iBnd++)
720 {
721 DNDS_assert(bndElemInfoSerial->operator()(iBnd, 0).zone != BC_ID_NULL);
722 }
723
724 log() << "OpenFOAM === got num bnd: " << nBnd << std::endl;
725
726 log() << "OpenFOAM === Serial Read Done" << std::endl;
727 }
728}
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
DNDS::Geom::tAdj1
decltype(tAdj1Pair::father) tAdj1
Definition Mesh_DeviceView.hpp:69
DNDS::Geom::FaceIDIsTrueInternal
DNDS_DEVICE_CALLABLE bool FaceIDIsTrueInternal(t_index id)
Definition BoundaryCondition.hpp:90
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::SerialReadAndDistribute
@ SerialReadAndDistribute
Definition Mesh.hpp:1080
DNDS::Geom::tAdj2
decltype(tAdj2Pair::father) tAdj2
Definition Mesh_DeviceView.hpp:71
DNDS::Geom::tAdj
decltype(tAdjPair::father) tAdj
Definition Mesh_DeviceView.hpp:65
DNDS::Geom::tElemInfoArray
DNDS::ssp< DNDS::ParArray< ElemInfo > > tElemInfoArray
Definition Mesh_DeviceView.hpp:81
DNDS::Geom::t_FBCName_2_ID
std::function< t_index(const std::string &)> t_FBCName_2_ID
Definition BoundaryCondition.hpp:46
DNDS::Geom::_getElemTypeFromCGNSType
constexpr Elem::ElemType _getElemTypeFromCGNSType(ElementType_t cgns_et)
Definition CGNS.hpp:16
DNDS::Geom::tCoord
decltype(tCoordPair::father) tCoord
Definition Mesh_DeviceView.hpp:79
DNDS::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:181
DNDS::rowsize
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:114
DNDS::UnInitReal
DNDS_CONSTANT const real UnInitReal
Sentinel "not initialised" real value (NaN). Cheap to detect with std::isnan or IsUnInitReal; survive...
Definition Defines.hpp:179
DNDS::IsUnInitReal
bool IsUnInitReal(real v)
Whether v equals the NaN sentinel UnInitReal (tested via isnan).
Definition Defines.hpp:188
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:112
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::Geom::ElemInfo::CommType
static MPI_Datatype CommType()
Definition Mesh_DeviceView.hpp:43
DNDS::Geom::Elem::Element::GetNumNodes
DNDS_DEVICE_CALLABLE constexpr t_index GetNumNodes() const
Definition Elements.hpp:235
DNDS::Geom::OpenFOAM::OpenFOAMConverter::cell2node
std::vector< std::vector< index > > cell2node
Definition OpenFOAMMesh.hpp:391
DNDS::Geom::OpenFOAM::OpenFOAMConverter::BuildCell2Face
void BuildCell2Face(OpenFOAMReader &reader)
Definition OpenFOAMMesh.hpp:412
DNDS::Geom::OpenFOAM::OpenFOAMConverter::BuildFaceElemInfo
void BuildFaceElemInfo(OpenFOAMReader &reader)
Definition OpenFOAMMesh.hpp:395
DNDS::Geom::OpenFOAM::OpenFOAMConverter::BuildCell2Node
void BuildCell2Node(OpenFOAMReader &reader)
Definition OpenFOAMMesh.hpp:447
DNDS::Geom::OpenFOAM::OpenFOAMReader::boundaryConditions
std::map< std::string, OpenFOAMBoundaryCondition > boundaryConditions
Definition OpenFOAMMesh.hpp:157
DNDS::Geom::OpenFOAM::OpenFOAMReader::ReadBoundary
void ReadBoundary(std::istream &IFS)
Definition OpenFOAMMesh.hpp:268
DNDS::Geom::OpenFOAM::OpenFOAMReader::ReadNeighbour
void ReadNeighbour(std::istream &IFS)
Definition OpenFOAMMesh.hpp:245
DNDS::Geom::OpenFOAM::OpenFOAMReader::faces
std::vector< std::vector< index > > faces
Definition OpenFOAMMesh.hpp:154
DNDS::Geom::UnstructuredMeshSerialRW::ReadFromCGNSSerial
void ReadFromCGNSSerial(const std::string &fName, const t_FBCName_2_ID &FBCName_2_ID)
reads a cgns file as serial input
Definition Mesh_Serial_ReadFromCGNS.cpp:247
DNDS::Geom::UnstructuredMeshSerialRW::mesh
DNDS::ssp< UnstructuredMesh > mesh
Definition Mesh.hpp:1117
DNDS::ObjName
Tag type for naming objects created via make_ssp.
Definition Defines.hpp:254
RowSize
for(DNDS::index iCell=0;iCell< 8;iCell++) for(DNDS j< res.parent2entity.father-> RowSize(iCell)
nBnd
DNDS::index nBnd
Definition test_MeshConnectivity_Interpolate.cpp:1456
i
res i
Definition test_MeshConnectivity_Interpolate.cpp:205
mesh
DistributedHex3D mesh
Definition test_MeshConnectivity_Interpolate.cpp:929
push_back
input explicitMaps push_back(EntityReorderMap{EntityKind::Cell, cellPartition})