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