DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Mesh_Serial_BuildCell2Cell.cpp
Go to the documentation of this file.
1#include "Mesh.hpp"
2
3#include <set>
4
5#include <nanoflann.hpp>
6
7#include "PointCloud.hpp"
8#include <unordered_map>
9#include <omp.h>
10
11namespace DNDS::Geom
12{
13
14 void UnstructuredMeshSerialRW::
15 Deduplicate1to1Periodic(real search_eps) // currently does not handle parallel input mode
16 {
17 DNDS_assert(this->dataIsSerialIn);
18 // TODO: build periodic donor oct-trees
19 // TODO: search in donors if periodic main and connect cell2cell
20
21 PointCloudKDTree periodicDonorCenter1;
22 PointCloudKDTree periodicDonorCenter2;
23 PointCloudKDTree periodicDonorCenter3;
24 std::vector<index> periodicDonorIBnd1;
25 std::vector<index> periodicDonorIBnd2;
26 std::vector<index> periodicDonorIBnd3;
27
28 for (DNDS::index iBnd = 0; iBnd < bnd2cellSerial->Size(); iBnd++)
29 {
30 auto faceID = bndElemInfoSerial->operator()(iBnd, 0).zone;
31 Eigen::Matrix<real, 3, Eigen::Dynamic> coordBnd;
32 {
33 coordBnd.resize(3, bnd2nodeSerial->RowSize(iBnd));
34 for (rowsize ib2n = 0; ib2n < bnd2nodeSerial->RowSize(iBnd); ib2n++)
35 coordBnd(EigenAll, ib2n) = coordSerial->operator[]((*bnd2nodeSerial)(iBnd, ib2n));
36 }
37 tPoint faceCent = coordBnd.rowwise().mean();
38
39 if (FaceIDIsPeriodicMain(faceID))
40 {
41 }
42 if (FaceIDIsPeriodicDonor(faceID))
43 {
44 if (faceID == BC_ID_PERIODIC_1_DONOR)
45 periodicDonorCenter1.pts.push_back(faceCent), periodicDonorIBnd1.push_back(iBnd);
46 if (faceID == BC_ID_PERIODIC_2_DONOR)
47 periodicDonorCenter2.pts.push_back(faceCent), periodicDonorIBnd2.push_back(iBnd);
48 if (faceID == BC_ID_PERIODIC_3_DONOR)
49 periodicDonorCenter3.pts.push_back(faceCent), periodicDonorIBnd3.push_back(iBnd);
50 }
51 }
52 index nDonor = periodicDonorCenter1.pts.size() +
53 periodicDonorCenter2.pts.size() +
54 periodicDonorCenter3.pts.size();
55 index nDonorAll{0};
56 MPI::Allreduce(&nDonor, &nDonorAll, 1, DNDS_MPI_INDEX, MPI_SUM, mesh->getMPI().comm);
57 if (nDonorAll)
58 mesh->isPeriodic = true; //! below are all periodic code
59 else
60 {
61 mesh->isPeriodic = false;
62 return;
63 }
64
65 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
66 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTree>,
67 PointCloudKDTree,
68 3,
69 index>;
70 kdtree_t periodicDonorTree1(3, periodicDonorCenter1);
71 kdtree_t periodicDonorTree2(3, periodicDonorCenter2);
72 kdtree_t periodicDonorTree3(3, periodicDonorCenter3);
73
74 std::unordered_map<index, index> periodicMainToDonor1;
75 std::unordered_map<index, index> periodicMainToDonor2;
76 std::unordered_map<index, index> periodicMainToDonor3;
77
78 std::unordered_map<index, index> iNodeDonorToMain1;
79 std::unordered_map<index, index> iNodeDonorToMain2;
80 std::unordered_map<index, index> iNodeDonorToMain3;
81
82 for (DNDS::index iBnd = 0; iBnd < bnd2cellSerial->Size(); iBnd++)
83 {
84 auto faceID = bndElemInfoSerial->operator()(iBnd, 0).zone;
85 Eigen::Matrix<real, 3, Eigen::Dynamic> coordBnd;
86 {
87 coordBnd.resize(3, bnd2nodeSerial->RowSize(iBnd));
88 for (rowsize ib2n = 0; ib2n < bnd2nodeSerial->RowSize(iBnd); ib2n++)
89 coordBnd(EigenAll, ib2n) = coordSerial->operator[]((*bnd2nodeSerial)(iBnd, ib2n));
90 }
91 tPoint faceCent = coordBnd.rowwise().mean();
92
93 if (FaceIDIsPeriodicMain(faceID))
94 {
95 tPoint faceCentTrans = mesh->periodicInfo.TransCoord(faceCent, faceID);
96 std::vector<index> outIndices(2);
97 Eigen::Vector<real, Eigen::Dynamic> outDistancesSqr;
98 outDistancesSqr.resize(2);
99 size_t nResults{0};
100 if (faceID == BC_ID_PERIODIC_1)
101 nResults = periodicDonorTree1.knnSearch(faceCentTrans.data(), 2, outIndices.data(), outDistancesSqr.data());
102 if (faceID == BC_ID_PERIODIC_2)
103 nResults = periodicDonorTree2.knnSearch(faceCentTrans.data(), 2, outIndices.data(), outDistancesSqr.data());
104 if (faceID == BC_ID_PERIODIC_3)
105 nResults = periodicDonorTree3.knnSearch(faceCentTrans.data(), 2, outIndices.data(), outDistancesSqr.data());
106 DNDS_assert_info(nResults >= 1, "Tree Search of periodic result number not enough");
107
108 // std::cout << outDistancesSqr[0] << std::endl;
109 DNDS_assert_info(outDistancesSqr[0] <= sqr(search_eps),
110 fmt::format("nearest neighbour dist{} not matching under the threshold",
111 outDistancesSqr[0]));
112 if (nResults > 1)
113 {
114 // DNDS_assert_info(outDistancesSqr[1] > sqr(search_eps),
115 // fmt::format("second nearest neighbour dist{} not matching over the threshold",
116 // outDistancesSqr[1]));
117 if (!(outDistancesSqr[1] > sqr(search_eps)))
118 std::cerr << TermColor::Red
119 << fmt::format("Warning: second nearest neighbour dist{} not matching over the threshold",
120 outDistancesSqr[1])
121 << "\n"
122 << fmt::format("at {:.8e}", Eigen::RowVectorFMTSafe<real, Eigen::Dynamic>(faceCentTrans.transpose()))
123 << TermColor::Reset
124 << std::endl;
125 }
126 index donorIBnd{-1};
127 if (faceID == BC_ID_PERIODIC_1)
128 donorIBnd = periodicDonorIBnd1.at(outIndices[0]), periodicMainToDonor1[iBnd] = donorIBnd;
129 if (faceID == BC_ID_PERIODIC_2)
130 donorIBnd = periodicDonorIBnd2.at(outIndices[0]), periodicMainToDonor2[iBnd] = donorIBnd;
131 if (faceID == BC_ID_PERIODIC_3)
132 donorIBnd = periodicDonorIBnd3.at(outIndices[0]), periodicMainToDonor3[iBnd] = donorIBnd;
133 Eigen::Matrix<real, 3, Eigen::Dynamic> coordBndOther;
134 {
135 coordBndOther.resize(3, bnd2nodeSerial->RowSize(donorIBnd));
136 for (rowsize ib2n = 0; ib2n < bnd2nodeSerial->RowSize(donorIBnd); ib2n++)
137 coordBndOther(EigenAll, ib2n) = // put onto main's data
138 coordSerial->operator[]((*bnd2nodeSerial)(donorIBnd, ib2n));
139 }
140 DNDS_assert(coordBndOther.cols() == coordBnd.cols());
141 for (rowsize ib2n = 0; ib2n < coordBnd.cols(); ib2n++)
142 {
143 bool found = false;
144 real minDist = veryLargeReal;
145 rowsize jb2nMin = 0;
146 for (rowsize jb2n = 0; jb2n < coordBnd.cols(); jb2n++)
147 {
148 real dist = (coordBndOther(EigenAll, jb2n) -
149 mesh->periodicInfo.TransCoord(coordBnd(EigenAll, ib2n), faceID))
150 .squaredNorm();
151 if (dist < minDist)
152 {
153 jb2nMin = jb2n;
154 minDist = dist;
155 }
156 if (dist < sqr(search_eps))
157 found = true;
158 }
159 if (faceID == BC_ID_PERIODIC_1)
160 iNodeDonorToMain1[(*bnd2nodeSerial)(donorIBnd, jb2nMin)] = (*bnd2nodeSerial)(iBnd, ib2n);
161 if (faceID == BC_ID_PERIODIC_2)
162 iNodeDonorToMain2[(*bnd2nodeSerial)(donorIBnd, jb2nMin)] = (*bnd2nodeSerial)(iBnd, ib2n);
163 if (faceID == BC_ID_PERIODIC_3)
164 iNodeDonorToMain3[(*bnd2nodeSerial)(donorIBnd, jb2nMin)] = (*bnd2nodeSerial)(iBnd, ib2n);
165 // DNDS_assert_info(found, fmt::format("min dist was: [{}]", minDist));
166 if (!found)
167 {
168 std::cerr << TermColor::Red << fmt::format("face-face node search min dist was: [{}]", minDist) << "\n"
169 << fmt::format("at {:.8e}", Eigen::RowVectorFMTSafe<real, Eigen::Dynamic>(faceCentTrans.transpose()))
170 << TermColor::Reset
171 << std::endl;
172 }
173 }
174 }
175 if (FaceIDIsPeriodicDonor(faceID))
176 {
177 }
178 }
179 /******************************************************************/
180 // need to be done before cell2node points to de-duplicated
181 /**********************************************************************************************************************/
182 // getting node2cell // only primary vertices
183 std::unordered_map<index, std::vector<index>> donorNode2Cell;
184 for (DNDS::index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
185 for (DNDS::rowsize iN = 0; iN < Elem::Element{(*cellElemInfoSerial)(iCell, 0).getElemType()}.GetNumVertices(); iN++)
186 {
187 auto iNode = (*cell2nodeSerial)(iCell, iN);
188 if (iNodeDonorToMain1.count(iNode) || iNodeDonorToMain2.count(iNode) || iNodeDonorToMain3.count(iNode))
189 {
190 donorNode2Cell[iNode].push_back(iCell);
191 }
192 }
193 std::unordered_map<index, std::vector<index>> donorNode2Bnd;
194 for (DNDS::index iBnd = 0; iBnd < bnd2nodeSerial->Size(); iBnd++)
195 for (DNDS::rowsize iN = 0; iN < Elem::Element{(*bndElemInfoSerial)(iBnd, 0).getElemType()}.GetNumVertices(); iN++)
196 {
197 auto iNode = (*bnd2nodeSerial)(iBnd, iN);
198 if (iNodeDonorToMain1.count(iNode) || iNodeDonorToMain2.count(iNode) || iNodeDonorToMain3.count(iNode))
199 {
200 donorNode2Bnd[iNode].push_back(iBnd);
201 }
202 }
203 /**********************************************************************************************************************/
204
205 cell2nodePbiSerial = make_ssp<decltype(cell2nodePbiSerial)::element_type>(ObjName{"Deduplicate1to1Periodic::cell2nodePbiSerial"}, mesh->getMPI());
206 cell2nodePbiSerial->Resize(cell2nodeSerial->Size());
207 for (index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
208 cell2nodePbiSerial->ResizeRow(iCell, cell2nodeSerial->RowSize(iCell));
209 // scan the periodic bnds
210 for (index iBnd = 0; iBnd < bnd2nodeSerial->Size(); iBnd++)
211 {
212 auto faceID = bndElemInfoSerial->operator()(iBnd, 0).zone;
213 if (FaceIDIsPeriodicDonor(faceID))
214 {
215 for (auto iNodeFace : (*bnd2nodeSerial)[iBnd])
216 {
217 DNDS_assert(donorNode2Cell.count(iNodeFace)); // node 2 cell means all the cells in bound are calculated
218 for (auto iCell : donorNode2Cell[iNodeFace])
219 for (auto iNode : (*bnd2nodeSerial)[iBnd])
220 for (rowsize ic2n = 0; ic2n < (*cell2nodeSerial).RowSize(iCell); ic2n++)
221 {
222 if ((*cell2nodeSerial)(iCell, ic2n) == iNode)
223 {
224 if (faceID == BC_ID_PERIODIC_1_DONOR)
225 (*cell2nodePbiSerial)(iCell, ic2n).setP1True();
226 if (faceID == BC_ID_PERIODIC_2_DONOR)
227 (*cell2nodePbiSerial)(iCell, ic2n).setP2True();
228 if (faceID == BC_ID_PERIODIC_3_DONOR)
229 (*cell2nodePbiSerial)(iCell, ic2n).setP3True();
230 }
231 }
232 }
233 }
234 }
235 bnd2nodePbiSerial = make_ssp<decltype(bnd2nodePbiSerial)::element_type>(ObjName{"Deduplicate1to1Periodic::bnd2nodePbiSerial"}, mesh->getMPI());
236 bnd2nodePbiSerial->Resize(bnd2nodeSerial->Size());
237 for (index iBnd = 0; iBnd < bnd2nodeSerial->Size(); iBnd++)
238 bnd2nodePbiSerial->ResizeRow(iBnd, bnd2nodeSerial->RowSize(iBnd));
239 // scan the periodic bnds
240 for (index iBnd = 0; iBnd < bnd2nodeSerial->Size(); iBnd++)
241 {
242 auto faceID = bndElemInfoSerial->operator()(iBnd, 0).zone;
243 if (FaceIDIsPeriodicDonor(faceID))
244 {
245 for (auto iNodeFace : (*bnd2nodeSerial)[iBnd])
246 {
247 DNDS_assert(donorNode2Bnd.count(iNodeFace));
248 for (auto iBndAdj : donorNode2Bnd[iNodeFace])
249 for (auto iNode : (*bnd2nodeSerial)[iBnd])
250 for (rowsize ic2n = 0; ic2n < (*bnd2nodeSerial).RowSize(iBndAdj); ic2n++)
251 {
252 if ((*bnd2nodeSerial)(iBndAdj, ic2n) == iNode)
253 {
254 if (faceID == BC_ID_PERIODIC_1_DONOR)
255 (*bnd2nodePbiSerial)(iBndAdj, ic2n).setP1True();
256 if (faceID == BC_ID_PERIODIC_2_DONOR)
257 (*bnd2nodePbiSerial)(iBndAdj, ic2n).setP2True();
258 if (faceID == BC_ID_PERIODIC_3_DONOR)
259 (*bnd2nodePbiSerial)(iBndAdj, ic2n).setP3True();
260 }
261 }
262 }
263 }
264 }
265 /******************************************************************/
266
267 std::vector<index> iNodeOld2New(coordSerial->Size(), 0);
268 index nNodeNew;
269
270 auto getNewNodeMap = [&](const decltype(iNodeDonorToMain1) &iNodeDonorToMainI) -> void
271 {
272 nNodeNew = 0;
273 for (index iNode = 0; iNode < coordSerial->Size(); iNode++)
274 {
275 if (iNodeOld2New[iNode] < 0)
276 continue;
277 bool isDonor = iNodeDonorToMainI.count(iNode);
278 if (!isDonor)
279 iNodeOld2New[iNode] = nNodeNew++;
280 }
281 for (const auto &p : iNodeDonorToMainI)
282 {
283 index iMain = p.second;
284 iNodeOld2New.at(p.first) = -1 - iMain; // minus means duplication here; leaves a self-pointing pointer in the minus region
285 }
286 };
287 getNewNodeMap(iNodeDonorToMain1);
288 getNewNodeMap(iNodeDonorToMain2);
289 getNewNodeMap(iNodeDonorToMain3);
290
291 auto track = [&](index iNode) -> index
292 {
293 while (iNodeOld2New.at(iNode) < 0)
294 iNode = -1 - iNodeOld2New[iNode];
295 return iNodeOld2New.at(iNode); // recursive tracker
296 };
297
298 for (index i = 0; i < cell2nodeSerial->Size(); i++)
299 for (auto &ii : (*cell2nodeSerial)[i])
300 ii = track(ii);
301 for (index i = 0; i < bnd2nodeSerial->Size(); i++)
302 for (auto &ii : (*bnd2nodeSerial)[i])
303 ii = track(ii);
304
305 // some assertions:65
306 // for (index iBnd = 0; iBnd < bnd2nodeSerial->Size(); iBnd++)
307 // {
308 // auto faceID = bndElemInfoSerial->operator()(iBnd, 0).zone;
309 // if (FaceIDIsPeriodicDonor(faceID))
310 // {
311
312 // }
313 // } //TODO
314
315 cell2nodePbiSerial->Compress();
316 bnd2nodePbiSerial->Compress();
317
318 decltype(coordSerial) coordSerialOld = coordSerial;
319 coordSerial = make_ssp<decltype(coordSerial)::element_type>(ObjName{"Deduplicate1to1Periodic::coordSerial"}, mesh->getMPI());
320 coordSerial->Resize(nNodeNew);
321 for (index i = 0; i < coordSerialOld->Size(); i++)
322 if (iNodeOld2New[i] >= 0)
323 (*coordSerial)[iNodeOld2New[i]] = (*coordSerialOld)[i];
324 coordSerialOld.reset();
325 }
326
327 void UnstructuredMeshSerialRW::
328 BuildCell2Cell() // currently does not handle parallel input mode
329 {
330 DNDS_assert(this->dataIsSerialIn);
331 real t0 = MPI_Wtime();
332 if (mesh->getMPI().rank == mRank)
333 DNDS::log() << "UnstructuredMeshSerialRW === Doing BuildCell2Cell "
334#ifdef DNDS_USE_OMP
335 << fmt::format("Using OMP [{}]", omp_get_max_threads())
336#endif
337 << std::endl;
338 cell2cellSerial = make_ssp<decltype(cell2cellSerial)::element_type>(ObjName{"BuildCell2Cell::cell2cellSerial"}, mesh->getMPI());
339 // if (mRank != mesh->getMPI().rank)
340 // return;
341 /// TODO: abstract these: invert cone (like node 2 cell -> cell 2 node) (also support operating on pair)
342 /**********************************************************************************************************************/
343 // getting node2cell // only primary vertices
344 std::vector<std::vector<DNDS::index>> node2cell(coordSerial->Size());
345 // std::vector<DNDS::rowsize> node2cellSz(coordSerial->Size(), 0);
346 // for (DNDS::index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
347 // for (DNDS::rowsize iN = 0; iN < Elem::Element{(*cellElemInfoSerial)(iCell, 0).getElemType()}.GetNumVertices(); iN++)
348 // node2cellSz[(*cell2nodeSerial)(iCell, iN)]++;
349 // for (DNDS::index iNode = 0; iNode < node2cell.size(); iNode++)
350 // node2cell[iNode].reserve(node2cellSz[iNode]);
351 for (DNDS::index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
352 for (DNDS::rowsize iN = 0; iN < Elem::Element{(*cellElemInfoSerial)(iCell, 0).getElemType()}.GetNumVertices(); iN++)
353 node2cell[(*cell2nodeSerial)(iCell, iN)].push_back(iCell);
354 // if (mesh->getMPI().rank == mRank)
355 // for (auto &row : node2cell)
356 // if (row.size() != 3)
357 // std::cout << "fucked not 4" << std::endl;
358 // node2cellSz.clear();
359 /**********************************************************************************************************************/
360 if (mesh->getMPI().rank == mRank)
361 DNDS::log() << "UnstructuredMeshSerialRW === Doing BuildCell2Cell Part 1" << std::endl;
362 cell2cellSerial->Resize(cell2nodeSerial->Size(), 32); // 32 is a guess for the number of nodal neighbors
363
364 index nCells = cell2nodeSerial->Size();
365 index nCellsDone = 0;
366#ifdef DNDS_USE_OMP
367# pragma omp parallel for
368#endif
369 for (DNDS::index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
370 {
371 auto CellElem = Elem::Element{(*cellElemInfoSerial)[iCell]->getElemType()};
372 std::vector<DNDS::index> cell2nodeRow{
373 (*cell2nodeSerial)[iCell].begin(),
374 (*cell2nodeSerial)[iCell].begin() + CellElem.GetNumVertices()};
375 std::sort(cell2nodeRow.begin(), cell2nodeRow.end());
376 // only primary vertices
377 std::unordered_set<DNDS::index> c_neighbors_hash;
378 c_neighbors_hash.reserve(64);
379 std::vector<DNDS::index>
380 c_neighbors;
381 c_neighbors.reserve(64);
382 // std::set<DNDS::index> c_neighbors; // could optimize?
383 // // std::vector<DNDS::index> c_neighbors;
384 // // c_neighbors.reserve(30);
385 // /****/
386#ifdef DNDS_USE_OMP
387 // #pragma omp single
388# pragma omp critical
389#endif
390 {
391 if (nCellsDone % (nCells / 1000 + 1) == 0)
392 {
393 auto fmt = DNDS::log().flags();
394 // DNDS::log() << "\r\033[K" << std::setw(5) << double(nCellsDone) / nCells * 100 << "%";
395 // DNDS::log().flush();
396 print_progress(log(), double(nCellsDone) / nCells);
397 DNDS::log().setf(fmt);
398 }
399
400 if (nCellsDone == nCells - 1)
401 DNDS::log()
402 // << "\r\033[K"
403 << std::endl;
404 }
405 /****/
406
407 for (auto iNode : cell2nodeRow)
408 for (auto iCellOther : node2cell[iNode])
409 {
410 if (iCellOther == iCell || c_neighbors_hash.count(iCellOther) == 1)
411 continue;
412 //! ** override: point-complete cell2cell info!
413 c_neighbors.push_back(iCellOther);
414 c_neighbors_hash.insert(iCellOther);
415 continue;
416 //! ** override: point-complete cell2cell info!
417
418 auto CellElemO = Elem::Element{(*cellElemInfoSerial)[iCellOther]->getElemType()};
419 std::vector<DNDS::index> cell2nodeRowO{
420 (*cell2nodeSerial)[iCellOther].begin(),
421 (*cell2nodeSerial)[iCellOther].begin() + CellElemO.GetNumVertices()};
422 std::sort(cell2nodeRowO.begin(), cell2nodeRowO.end());
423 std::vector<DNDS::index> interSect;
424 interSect.reserve(32);
425 std::set_intersection(cell2nodeRowO.begin(), cell2nodeRowO.end(), cell2nodeRow.begin(), cell2nodeRow.end(),
426 std::back_inserter(interSect));
427 if ((iCellOther != iCell) && interSect.size() > (mesh->dim - 1)) //! mesh->dim - 1 is a simple method
428
429 // (false ||
430 // Elem::cellsAreFaceConnected(
431 // cell2nodeSerial->operator[](iCell),
432 // cell2nodeSerial->operator[](iCellOther),
433 // CellElem,
434 // Elem::Element{(*cellElemInfoSerial)[iCellOther]->getElemType()}))
435 {
436 c_neighbors_hash.insert(iCellOther);
437 c_neighbors.push_back(iCellOther);
438 }
439 }
440 // if (!c_neighbors.erase(iCell))
441 // DNDS_assert_info(false, "neighbors should include myself now");
442 // std::sort(c_neighbors.begin(), c_neighbors.end());
443 // auto last = std::unique(c_neighbors.begin(), c_neighbors.end());
444 // c_neighbors.erase(last, c_neighbors.end());
445
446 std::sort(c_neighbors.begin(), c_neighbors.end());
447 cell2cellSerial->ResizeRow(iCell, c_neighbors.size());
448 DNDS::rowsize ic2c = 0;
449 for (auto iCellOther : c_neighbors)
450 (*cell2cellSerial)(iCell, ic2c++) = iCellOther;
451#ifdef DNDS_USE_OMP
452# pragma omp atomic
453#endif
454 nCellsDone++;
455 }
456 (*cell2cellSerial).Compress();
457 if (mesh->getMPI().rank == mRank)
458 log() << std::endl;
459
460 /*************************************************************************************************/
461 cell2cellSerialFacial = make_ssp<decltype(cell2cellSerialFacial)::element_type>(ObjName{"BuildCell2Cell::cell2cellSerialFacial"}, mesh->getMPI());
462 if (mesh->getMPI().rank == mRank)
463 DNDS::log() << "UnstructuredMeshSerialRW === Doing BuildCell2Cell Part 2" << std::endl;
464 cell2cellSerialFacial->Resize(cell2cellSerial->Size(), 6);
465 nCellsDone = 0;
466#ifdef DNDS_USE_OMP
467# pragma omp parallel for
468#endif
469 for (index iCell = 0; iCell < cell2cellSerial->Size(); iCell++)
470 {
471#ifdef DNDS_USE_OMP
472# pragma omp critical
473#endif
474 {
475 if (nCellsDone % (nCells / 1000 + 1) == 0)
476 {
477 auto fmt = DNDS::log().flags();
478 // DNDS::log() << "\r\033[K" << std::setw(5) << double(nCellsDone) / nCells * 100 << "%";
479 // DNDS::log().flush();
480 print_progress(log(), double(nCellsDone) / nCells);
481 DNDS::log().setf(fmt);
482 }
483
484 if (nCellsDone == nCells - 1)
485 DNDS::log()
486 // << "\r\033[K"
487 << std::endl;
488 }
489 std::vector<index> facialNeighbors;
490 facialNeighbors.reserve(6);
491 std::vector<index> c2ni((*cell2nodeSerial)[iCell].begin(),
492 (*cell2nodeSerial)[iCell].begin() +
493 Elem::Element{(*cellElemInfoSerial)[iCell]->getElemType()}.GetNumVertices());
494 std::sort(c2ni.begin(), c2ni.end());
495 for (auto iCellOther : (*cell2cellSerial)[iCell])
496 {
497 std::vector<index> c2nj((*cell2nodeSerial)[iCellOther].begin(),
498 (*cell2nodeSerial)[iCellOther].begin() +
499 Elem::Element{(*cellElemInfoSerial)[iCellOther]->getElemType()}.GetNumVertices());
500 std::sort(c2nj.begin(), c2nj.end());
501 std::vector<index> intersect;
502 intersect.reserve(9);
503 std::set_intersection(c2ni.begin(), c2ni.end(), c2nj.begin(), c2nj.end(), std::back_inserter(intersect));
504 if (intersect.size() >= mesh->dim) // for 2d, exactly 2; for 3d, 3 or 4
505 facialNeighbors.push_back(iCellOther);
506 }
507 (*cell2cellSerialFacial).ResizeRow(iCell, facialNeighbors.size());
508 (*cell2cellSerialFacial)[iCell] = facialNeighbors;
509#ifdef DNDS_USE_OMP
510# pragma omp atomic
511#endif
512 nCellsDone++;
513 }
514 if (mesh->getMPI().rank == mRank)
515 log() << std::endl;
516 /*************************************************************************************************/
517
518 if (mesh->getMPI().rank == mRank)
519 DNDS::log() << "UnstructuredMeshSerialRW === Done BuildCell2Cell"
520 << fmt::format(" time [{:.4g}]", MPI_Wtime() - t0) << std::endl;
521 }
522
523 void UnstructuredMesh::RecreatePeriodicNodes()
524 {
525 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
526 DNDS_assert(cell2node.trans.pLGhostMapping);
527 DNDS_assert(coords.trans.pLGhostMapping && coords.trans.pLGlobalMapping);
528 if (!isPeriodic)
529 {
530 return;
531 }
532 tAdj1Pair nodeNeedCreate;
533 nodeNeedCreate.InitPair("RecreatePeriodicNodes::nodeNeedCreate", mpi);
534 nodeNeedCreate.TransAttach();
535 nodeNeedCreate.father->Resize(coords.father->Size());
536 nodeNeedCreate.trans.BorrowGGIndexing(coords.trans);
537 nodeNeedCreate.trans.createMPITypes();
538 // std::cout << "here X0.5" << std::endl;
539 for (index i = 0; i < nodeNeedCreate.Size(); i++)
540 nodeNeedCreate(i, 0) = 0x01U;
541 for (index iC = 0; iC < cell2node.father->Size(); iC++)
542 for (rowsize ic2n = 0; ic2n < cell2node.RowSize(iC); ic2n++)
543 nodeNeedCreate(cell2node(iC, ic2n), 0) |= (0x01U << uint8_t(cell2nodePbi(iC, ic2n)));
544 DNDS::ArrayTransformerType<tAdj1::element_type>::Type nodeNeedCreatePastTrans;
545 tAdj1 nodeNeedCreatePast;
546 nodeNeedCreatePast = make_ssp<decltype(nodeNeedCreatePast)::element_type>(ObjName{"RecreatePeriodicNodes::nodeNeedCreatePast"}, mpi);
547 nodeNeedCreatePastTrans.setFatherSon(nodeNeedCreate.son, nodeNeedCreatePast);
548 nodeNeedCreatePastTrans.createFatherGlobalMapping();
549 std::vector<index> pushSonSeries(nodeNeedCreate.son->Size());
550 for (index i = 0; i < pushSonSeries.size(); i++)
551 pushSonSeries[i] = i;
552 nodeNeedCreatePastTrans.createGhostMapping(pushSonSeries, nodeNeedCreate.trans.pLGhostMapping->ghostStart);
553 nodeNeedCreatePastTrans.createMPITypes();
554 nodeNeedCreatePastTrans.pullOnce();
555 for (index i = 0; i < nodeNeedCreatePast->Size(); i++)
556 {
557 index iNodeG = nodeNeedCreate.trans.pLGhostMapping->pushingIndexGlobal[i]; //?should be right
558 auto [ret, rank, iNode] = nodeNeedCreate.trans.pLGhostMapping->search_indexAppend(iNodeG);
559 DNDS_assert(rank == -1); // must be in main
560 nodeNeedCreate.father->operator()(iNode, 0) |= nodeNeedCreatePast->operator()(i, 0);
561 }
562 // std::cout << "here X1" << std::endl;
563
564 index nCreatedNodes{0};
565 tAdj8Pair node2recreatedNodes;
566 node2recreatedNodes.InitPair("RecreatePeriodicNodes::node2recreatedNodes", mpi);
567 node2recreatedNodes.TransAttach();
568 node2recreatedNodes.father->Resize(coords.father->Size());
569 node2recreatedNodes.trans.BorrowGGIndexing(coords.trans);
570 node2recreatedNodes.trans.createMPITypes();
571 for (index i = 0; i < coords.father->Size(); i++)
572 {
573 DNDS_assert(nodeNeedCreate(i, 0) & 0x01U);
574 node2recreatedNodes(i, 0) = i;
575 for (uint8_t bitCur = 1; bitCur < 8; bitCur++)
576 {
577 if (nodeNeedCreate(i, 0) & (0x01U << bitCur))
578 node2recreatedNodes(i, rowsize(bitCur)) = coords.father->Size() + nCreatedNodes++;
579 else
580 node2recreatedNodes(i, rowsize(bitCur)) = UnInitIndex;
581 }
582 }
583 // node2recreatedNodes now points to local new coords
584
585 coordsPeriodicRecreated.InitPair("RecreatePeriodicNodes::coordsPeriodicRecreated", mpi);
586 coordsPeriodicRecreated.TransAttach();
587 coordsPeriodicRecreated.father->Resize(nCreatedNodes + coords.father->Size());
588 coordsPeriodicRecreated.trans.createFatherGlobalMapping();
589 nodeRecreated2nodeLocal.resize(coordsPeriodicRecreated.father->Size());
590 for (index iN = 0; iN < coords.father->Size(); iN++) // fill the coords
591 {
592 coordsPeriodicRecreated[iN] = coords[iN];
593 nodeRecreated2nodeLocal[iN] = iN;
594 for (uint8_t bitCur = 1; bitCur < 8; bitCur++)
595 {
596 if (node2recreatedNodes(iN, rowsize(bitCur)) != UnInitIndex)
597 {
598 coordsPeriodicRecreated[node2recreatedNodes(iN, rowsize(bitCur))] =
599 periodicInfo.GetCoordByBits(coords[iN], NodePeriodicBits{bitCur});
600 nodeRecreated2nodeLocal[node2recreatedNodes(iN, rowsize(bitCur))] = iN;
601 }
602 }
603 }
604 for (index iN = 0; iN < coords.father->Size(); iN++) // convert node2recreatedNodes into global
605 for (int i = 0; i < node2recreatedNodes.RowSize(iN); i++)
606 if (node2recreatedNodes(iN, i) != UnInitIndex)
607 node2recreatedNodes(iN, i) = coordsPeriodicRecreated.trans.pLGlobalMapping->operator()(
608 mpi.rank, node2recreatedNodes(iN, i));
609 node2recreatedNodes.trans.pullOnce(); // for cell2node query
610 // std::cout << "here X2" << std::endl;
611 cell2nodePeriodicRecreated.InitPair("RecreatePeriodicNodes::cell2nodePeriodicRecreated", mpi);
612 cell2nodePeriodicRecreated.TransAttach();
613 cell2nodePeriodicRecreated.father->Resize(cell2node.father->Size());
614 for (index iC = 0; iC < cell2node.father->Size(); iC++)
615 {
616 cell2nodePeriodicRecreated.ResizeRow(iC, cell2node.RowSize(iC));
617 for (rowsize ic2n = 0; ic2n < cell2node.RowSize(iC); ic2n++)
618 {
619 index iNode = cell2node(iC, ic2n);
620 index iNodeNewGlobal = node2recreatedNodes(iNode, 0);
621 if (cell2nodePbi(iC, ic2n))
622 {
623 iNodeNewGlobal = node2recreatedNodes(iNode, rowsize(uint8_t(cell2nodePbi(iC, ic2n))));
624 DNDS_assert(iNodeNewGlobal != UnInitIndex);
625 }
626 cell2nodePeriodicRecreated(iC, ic2n) = iNodeNewGlobal;
627 }
628 }
629 // std::cout << "here X3" << std::endl;
630 // here, cell2nodePeriodicRecreated point to global coordsPeriodicRecreated
631 // coordsPeriodicRecreated has no son
632
633 index nNodeRecreatedGlobal = coordsPeriodicRecreated.father->globalSize();
634 if (mpi.rank == mRank)
635 log() << fmt::format("=== RecreatePeriodicNodes: new nNodes [{}]", nNodeRecreatedGlobal) << std::endl;
636 }
637
638 void UnstructuredMesh::BuildVTKConnectivity()
639 {
640 DNDS_assert(cell2node.father);
641 DNDS_assert(coords.trans.pLGhostMapping); // for this->NodeIndexLocal2Global()
642 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
643 if (isPeriodic)
644 {
645 DNDS_assert(cell2nodePeriodicRecreated.father);
646 DNDS_assert(coordsPeriodicRecreated.father);
647 }
648 vtkCell2nodeOffsets.resize(cell2node.father->Size() + 1);
649 vtkCell2nodeOffsets[0] = 0;
650 vtkCellType.resize(cell2node.father->Size());
651 for (index iC = 0; iC < cell2node.father->Size(); iC++)
652 {
653 auto eCell = GetCellElement(iC);
654 auto [vtkType, vtkC2n] = Elem::ToVTKVertsAndData(eCell, cell2node[iC]);
655 vtkCellType[iC] = vtkType;
656 vtkCell2nodeOffsets[iC + 1] = vtkCell2nodeOffsets[iC] + vtkC2n.size();
657 }
658 vtkCell2node.reserve(vtkCell2nodeOffsets.back()); // now offsets are local
659 for (index iC = 0; iC < cell2node.father->Size(); iC++)
660 {
661 auto eCell = GetCellElement(iC);
662 if (!isPeriodic)
663 {
664 auto [vtkType, vtkC2n] = Elem::ToVTKVertsAndData(eCell, cell2node[iC]);
665 for (auto iN : vtkC2n)
666 vtkCell2node.push_back(this->NodeIndexLocal2Global(iN));
667 }
668 else
669 { // cell2nodePeriodicRecreated points to global
670 auto [vtkType, vtkC2n] = Elem::ToVTKVertsAndData(eCell, cell2nodePeriodicRecreated[iC]);
671 for (auto iN : vtkC2n)
672 vtkCell2node.push_back(iN);
673 }
674 }
675 index localVtkCellOffset = vtkCell2nodeOffsets.back();
676 index sumPrevVtkCellOffset{0};
677 MPI::Scan(&localVtkCellOffset, &sumPrevVtkCellOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
678 sumPrevVtkCellOffset -= localVtkCellOffset;
679 for (auto &v : vtkCell2nodeOffsets)
680 v += sumPrevVtkCellOffset;
681
682 // get some offsets
683 index curNNode = coords.father->Size();
684 if (isPeriodic)
685 curNNode = coordsPeriodicRecreated.father->Size();
686 MPI::Scan(&curNNode, &vtkNodeOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
687 vtkNodeOffset -= curNNode;
688 index curNCell = cell2node.father->Size();
689 MPI::Scan(&curNCell, &vtkCellOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
690 vtkCellOffset -= curNCell;
691
692 index vtkCell2NodeSize = vtkCell2node.size();
693 MPI::Allreduce(&vtkCell2NodeSize, &vtkCell2NodeGlobalSiz, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
694
695 index nCellsLocal = cell2node.father->Size();
696 MPI_Allreduce(&nCellsLocal, &vtkNCellGlobal, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
697 tCoord coordOut = isPeriodic ? coordsPeriodicRecreated.father : coords.father;
698 index nNodesLocal = coordOut->Size();
699 MPI_Allreduce(&nNodesLocal, &vtkNNodeGlobal, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
700
701 // std::cout << mpi.rank << " " << sumPrevVtkCellOffset << std::endl;
702 }
703}
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::ArrayTransformer
Ghost-communication engine for a father / son ParArray pair.
Definition ArrayTransformer.hpp:428
DNDS::Geom::Elem::ToVTKVertsAndData
std::pair< int, std::vector< index > > ToVTKVertsAndData(Element e, const TIn &vin)
Definition Elements.hpp:566
DNDS::Geom::tAdj1
decltype(tAdj1Pair::father) tAdj1
Definition Mesh_DeviceView.hpp:69
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::Geom::FaceIDIsPeriodicDonor
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicDonor(t_index id)
Definition BoundaryCondition.hpp:130
DNDS::Geom::FaceIDIsPeriodicMain
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicMain(t_index id)
Definition BoundaryCondition.hpp:123
DNDS::Geom::tCoord
decltype(tCoordPair::father) tCoord
Definition Mesh_DeviceView.hpp:79
DNDS::MPI::Scan
MPI_int Scan(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Scan (inclusive prefix reduction).
Definition MPI.cpp:220
DNDS::MPI::Allreduce
MPI_int Allreduce(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Allreduce.
Definition MPI.cpp:203
DNDS::TermColor::Reset
constexpr std::string_view Reset
ANSI escape: reset all attributes.
Definition Defines.hpp:883
DNDS::TermColor::Red
constexpr std::string_view Red
ANSI escape: bright red foreground.
Definition Defines.hpp:869
DNDS::DNDS_MPI_INDEX
const MPI_Datatype DNDS_MPI_INDEX
MPI datatype matching index (= MPI_INT64_T).
Definition MPI.hpp:90
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::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:107
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:517
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::print_progress
void print_progress(std::ostream &os, double progress)
Render a textual progress bar to os for progress in [0, 1].
Definition Defines.cpp:75
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:49
DNDS::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:190
DNDS::ArrayPair
Convenience bundle of a father, son, and attached ArrayTransformer.
Definition ArrayPair.hpp:163
DNDS::ArrayPair::TransAttach
void TransAttach()
Bind the transformer to the current father / son pointers.
Definition ArrayPair.hpp:301
DNDS::ArrayPair::father
ssp< TArray > father
Owned-side array (must be resized before ghost setup).
Definition ArrayPair.hpp:170
DNDS::ArrayPair::InitPair
void InitPair(const std::string &name, Args &&...args)
Allocate both father and son arrays, forwarding all args to TArray constructor.
Definition ArrayPair.hpp:322
DNDS::ArrayPair::RowSize
auto RowSize() const
Uniform row width (delegates to father).
Definition ArrayPair.hpp:258
DNDS::ArrayPair::ResizeRow
void ResizeRow(index i, rowsize rs)
Resize a single row in the combined address space.
Definition ArrayPair.hpp:273
DNDS::ArrayPair::trans
TTrans trans
Ghost-communication engine bound to father and son.
Definition ArrayPair.hpp:175
DNDS::Geom::Elem::Element::GetNumVertices
DNDS_DEVICE_CALLABLE constexpr t_index GetNumVertices() const
Definition Elements.hpp:209
DNDS::Geom::Periodicity::GetCoordByBits
DNDS_DEVICE_CALLABLE tPoint GetCoordByBits(const tPoint &c, const NodePeriodicBits &bits) const
Definition PeriodicInfo.hpp:320
DNDS::Geom::PointCloudKDTree::pts
std::vector< tPoint > pts
Definition PointCloud.hpp:11
DNDS::Geom::UnstructuredMeshSerialRW::mesh
DNDS::ssp< UnstructuredMesh > mesh
Definition Mesh.hpp:856
DNDS::Geom::UnstructuredMeshSerialRW::BuildCell2Cell
void BuildCell2Cell()
build cell2cell topology, with node-neighbors included
Definition Mesh_Serial_BuildCell2Cell.cpp:328
DNDS::Geom::UnstructuredMesh::t_deviceView
UnstructuredMeshDeviceView< B > t_deviceView
Definition Mesh.hpp:777
DNDS::Geom::UnstructuredMesh::coordsPeriodicRecreated
tCoordPair coordsPeriodicRecreated
Definition Mesh.hpp:141
DNDS::Geom::UnstructuredMesh::GetCellElement
Elem::Element GetCellElement(index iC)
Definition Mesh.hpp:480
DNDS::Geom::UnstructuredMesh::cell2nodePbi
tPbiPair cell2nodePbi
periodic only, after reader
Definition Mesh.hpp:65
DNDS::Geom::UnstructuredMesh::vtkCellType
std::vector< uint8_t > vtkCellType
Definition Mesh.hpp:133
DNDS::Geom::UnstructuredMesh::cell2nodePeriodicRecreated
tAdjPair cell2nodePeriodicRecreated
Definition Mesh.hpp:140
DNDS::Geom::UnstructuredMesh::vtkCell2nodeOffsets
std::vector< index > vtkCell2nodeOffsets
for parallel out
Definition Mesh.hpp:132
DNDS::Geom::UnstructuredMesh::coords
tCoordPair coords
reader
Definition Mesh.hpp:54
DNDS::Geom::UnstructuredMesh::adjPrimaryState
MeshAdjState adjPrimaryState
Definition Mesh.hpp:41
DNDS::Geom::UnstructuredMesh::vtkCell2node
std::vector< index > vtkCell2node
Definition Mesh.hpp:134
DNDS::Geom::UnstructuredMesh::nodeRecreated2nodeLocal
std::vector< index > nodeRecreated2nodeLocal
Definition Mesh.hpp:142
DNDS::MPIInfo::rank
int rank
This rank's 0-based index within comm (-1 until initialised).
Definition MPI.hpp:219
DNDS::MPIInfo::comm
MPI_Comm comm
The underlying MPI communicator handle.
Definition MPI.hpp:217
DNDS::ObjName
Tag type for naming objects created via make_ssp.
Definition Defines.hpp:249
Eigen::MatrixFMTSafe
Eigen::Matrix wrapper that hides begin/end from fmt.
Definition EigenUtil.hpp:62
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468