DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Mesh_Elevation.cpp
Go to the documentation of this file.
1#include "DNDS/Defines.hpp"
2#include "Mesh.hpp"
3#include "Geom/Quadrature.hpp"
4#include "DNDS/ArrayDerived/ArrayEigenUniMatrixBatch.hpp"
5#include "Geom/RadialBasisFunction.hpp"
6
7#include <fmt/core.h>
8#include <functional>
9#include <unordered_set>
10#include <Solver/Linear.hpp>
11
12#include "Geom/PointCloud.hpp"
13#include <nanoflann.hpp>
14
15#include "DNDS/EigenPCH.hpp"
16#ifdef DNDS_USE_SUPERLU
17# include <superlu_ddefs.h>
18#endif
19
20namespace DNDS::Geom
21{
22 void UnstructuredMesh::BuildO2FromO1Elevation(UnstructuredMesh &meshO1)
23 {
24 real epsSqrDist = 1e-20;
25
26 bool O1MeshIsO1 = meshO1.IsO1();
27 DNDS_assert(O1MeshIsO1);
28 // std::cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXX Phase0" << std::endl;
29 DNDS_assert(meshO1.adjPrimaryState == Adj_PointToLocal);
30 DNDS_assert(meshO1.cell2node.isLocal() && meshO1.bnd2node.isLocal());
31 mRank = meshO1.mRank;
32 mpi = meshO1.mpi;
33 dim = meshO1.dim;
34 isPeriodic = meshO1.isPeriodic;
35 periodicInfo = meshO1.periodicInfo;
36 nNodeO1 = meshO1.coords.father->Size();
37 elevState = Elevation_O1O2;
38
39 tAdjPair cellNewNodes; // records iNode - iNodeO1
40 index localNewNodeNum{0};
41 cellNewNodes.InitPair("BuildO2FromO1Elevation::cellNewNodes", mpi);
42 cellNewNodes.father->Resize(meshO1.cell2node.father->Size());
43 std::vector<Eigen::Vector<real, 3>> newCoords;
44
45 // std::cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXX Phase1" << std::endl;
46
47 /**********************************/ // TODO
48 // add new nodes to cellNewNodes, edge and face O2 nodes belong to smallest neighboring cell, cellNewNodes record local node indices
49 for (index iCell = 0; iCell < meshO1.cell2node.father->Size(); iCell++)
50 {
51 Elem::Element eCell = meshO1.GetCellElement(iCell);
52 auto c2n = meshO1.cell2node[iCell];
53 auto c2c = meshO1.cell2cell[iCell];
54 index iCellGlob = meshO1.cell2node.trans.pLGhostMapping->operator()(-1, iCell);
55
56 std::vector<index> cellNewNodeRowV;
57 cellNewNodeRowV.resize(eCell.GetNumElev_O1O2());
58 SmallCoordsAsVector coordsC;
59 meshO1.GetCoordsOnCell(iCell, coordsC);
60
61 for (index iNNode = 0; iNNode < eCell.GetNumElev_O1O2(); iNNode++)
62 {
63 auto eSpan = eCell.ObtainElevNodeSpan(iNNode);
64 std::array<index, Elem::CellNumNodeMax> spanNodes{UnInitIndex};
65 eCell.ExtractElevNodeSpanNodes(iNNode, c2n, spanNodes);
66 auto spanNodesSrt = spanNodes;
67 std::sort(spanNodesSrt.begin(), spanNodesSrt.begin() + eSpan.GetNumNodes());
68 tPoint newNodeCoord;
69
70 { // the new node method, should cope with shape func
71 SmallCoordsAsVector coordsSpan;
72 coordsSpan.resize(3, eSpan.GetNumNodes());
73 eCell.ExtractElevNodeSpanNodes(iNNode, coordsC, coordsSpan);
74 newNodeCoord = coordsSpan.rowwise().mean();
75
76 if (isPeriodic)
77 {
78 std::vector<NodePeriodicBits> pbi;
79 pbi.resize(eSpan.GetNumNodes());
80 eCell.ExtractElevNodeSpanNodes(iNNode, meshO1.cell2nodePbi[iCell], pbi);
81 NodePeriodicBits pbiR;
82 pbiR.setP1True(), pbiR.setP2True(), pbiR.setP3True();
83 for (auto pbie : pbi)
84 pbiR = pbiR & pbie;
85 newNodeCoord = periodicInfo.GetCoordBackByBits(newNodeCoord, pbiR);
86 }
87 }
88
89 index curMinGlobiCell = iCellGlob;
90 index curMinic2c = -1;
91 for (int ic2c = 0; ic2c < c2c.size(); ic2c++)
92 {
93 index iCellOther = c2c[ic2c];
94 if (iCellOther == iCell)
95 continue;
96 std::vector<index> c2nOther = meshO1.cell2node[iCellOther];
97 index iCellGlobOther = meshO1.cell2node.trans.pLGhostMapping->operator()(-1, iCellOther);
98 std::sort(c2nOther.begin(), c2nOther.end());
99 if (iCellGlobOther < curMinGlobiCell &&
100 std::includes(c2nOther.begin(), c2nOther.end(), spanNodesSrt.begin(), spanNodesSrt.begin() + eSpan.GetNumNodes()))
101 {
102 curMinGlobiCell = iCellGlobOther;
103 curMinic2c = ic2c;
104 }
105 }
106 if (curMinGlobiCell == iCellGlob)
107 {
108 newCoords.push_back(newNodeCoord);
109 cellNewNodeRowV.at(iNNode) = (localNewNodeNum++); //* recorded as local idx in newCoords
110 }
111 else
112 {
113 cellNewNodeRowV.at(iNNode) = (-1 - curMinic2c); //* recorded as -1 - ic2c that points to the owner
114 }
115 }
116 cellNewNodes.ResizeRow(iCell, cellNewNodeRowV.size());
117 cellNewNodes[iCell] = cellNewNodeRowV;
118 // std::cout << "iCell " << iCell << ": ";
119 // for (auto i : cellNewNodeRowV)
120 // std::cout << i << ", ";
121 // std::cout << std::endl;
122 }
123 // std::cout << fmt::format("Num NewNode {} ", localNewNodeNum) << std::endl;
124 // for (auto v : newCoords)
125 // std::cout << v.transpose() << std::endl;
126 index numNewNode = UnInitIndex;
127 MPI::Allreduce(&localNewNodeNum, &numNewNode, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
128 if (mpi.rank == mRank)
129 log() << fmt::format("=== Mesh Elevation: Num NewNode {} ", numNewNode) << std::endl;
130 index localNewNodeNumOffset{0};
131 MPI::Scan(&localNewNodeNum, &localNewNodeNumOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
132 if (mpi.rank == mpi.size - 1)
133 DNDS_assert(localNewNodeNumOffset == numNewNode);
134 localNewNodeNumOffset -= localNewNodeNum;
135 index numNodeO1Global = meshO1.NumNodeGlobal();
136
137 /**********************************/
138 // build each proc coordO2 from cellNewNodes
139 coords.InitPair("BuildO2FromO1Elevation::coords", mpi);
140 coords.father->Resize(meshO1.coords.father->Size() + localNewNodeNum);
141 for (index iC = 0; iC < meshO1.coords.father->Size(); iC++)
142 coords[iC] = meshO1.coords[iC];
143 for (index iC = meshO1.coords.father->Size(); iC < meshO1.coords.father->Size() + localNewNodeNum; iC++)
144 coords[iC] = newCoords.at(iC - meshO1.coords.father->Size());
145 coords.father->createGlobalMapping();
146
147 node2nodeOrig.InitPair("BuildO2FromO1Elevation::node2nodeOrig", mpi);
148 node2nodeOrig.father->Resize(meshO1.coords.father->Size() + localNewNodeNum);
149 for (index iC = 0; iC < meshO1.coords.father->Size(); iC++)
150 node2nodeOrig[iC] = meshO1.node2nodeOrig[iC];
151 for (index iC = meshO1.coords.father->Size(); iC < meshO1.coords.father->Size() + localNewNodeNum; iC++)
152 node2nodeOrig[iC][0] = // new node does not correspond to any "original", so assigned sequential new space
153 numNodeO1Global +
154 (iC - meshO1.coords.father->Size() + localNewNodeNumOffset);
155
156 // tAdj1Pair nodeLocalIdxOld;
157 // DNDS_MAKE_SSP(nodeLocalIdxOld.father, mpi);
158 // DNDS_MAKE_SSP(nodeLocalIdxOld.son, mpi);
159 // nodeLocalIdxOld.father->Resize(meshO1.coords.father->Size());
160 // for (index iC = 0; iC < meshO1.coords.father->Size(); iC++)
161 // nodeLocalIdxOld[iC][0] = iC;
162 // nodeLocalIdxOld.TransAttach();
163 // nodeLocalIdxOld.trans.BorrowGGIndexing(mesh->cell2node.trans);
164 // nodeLocalIdxOld.trans.createMPITypes();
165 // nodeLocalIdxOld.trans.pullOnce();
166
167 /**********************************/
168 // from coordO2 get cellNewNodesGlobal, and get cellNewNodesGlobalPair
169 for (index iCell = 0; iCell < meshO1.cell2node.father->Size(); iCell++)
170 {
171 auto cellNewNodesRow = cellNewNodes[iCell];
172 for (auto &iNodeNew : cellNewNodesRow)
173 {
174 // nodes here must be at local proc; now use global idx
175 iNodeNew =
176 iNodeNew < 0
177 ? iNodeNew // negative field was filled with indicating which cell to seek the node
178 : coords.father->pLGlobalMapping->operator()(mpi.rank, iNodeNew + meshO1.coords.father->Size());
179 }
180 }
181 cellNewNodes.TransAttach();
182 cellNewNodes.trans.BorrowGGIndexing(meshO1.cell2node.trans);
183 cellNewNodes.trans.createMPITypes();
184 cellNewNodes.trans.pullOnce();
185
186 coords.TransAttach();
187 std::vector<DNDS::index> ghostNodesTmp;
188 for (index iCell = 0; iCell < meshO1.cell2node.Size(); iCell++)
189 {
190 auto cellNewNodesRow = cellNewNodes[iCell];
191 for (auto iNodeNew : cellNewNodesRow)
192 {
193 if (iNodeNew < 0)
194 continue;
195 MPI_int rank = UnInitMPIInt;
196 index val = UnInitIndex;
197 if (!coords.trans.pLGlobalMapping->search(iNodeNew, rank, val))
198 DNDS_assert_info(false, "search failed");
199 if (rank != mpi.rank)
200 ghostNodesTmp.push_back(iNodeNew);
201 }
202 }
203 coords.trans.createGhostMapping(ghostNodesTmp);
204 coords.trans.createMPITypes();
205 coords.trans.pullOnce();
206
207 node2nodeOrig.TransAttach();
208 node2nodeOrig.trans.BorrowGGIndexing(coords.trans);
209 node2nodeOrig.trans.createMPITypes();
210 node2nodeOrig.trans.pullOnce();
211
212 /**********************************/
213 // each cell obtain new global cell2node global state with cellNewNodesGlobalPair
214 cell2node.InitPair("BuildO2FromO1Elevation::cell2node", mpi);
215 cellElemInfo.InitPair("BuildO2FromO1Elevation::cellElemInfo", ElemInfo::CommType(), ElemInfo::CommMult(), mpi);
216 if (isPeriodic)
217 {
218 cell2nodePbi.InitPair("BuildO2FromO1Elevation::cell2nodePbi", NodePeriodicBits::CommType(), NodePeriodicBits::CommMult(), getMPI());
219 }
220 cell2cellOrig.InitPair("BuildO2FromO1Elevation::cell2cellOrig", mpi);
221
222 cell2node.father->Resize(meshO1.cell2node.father->Size());
223 cellElemInfo.father->Resize(meshO1.cell2node.father->Size());
224 if (isPeriodic)
225 cell2nodePbi.father->Resize(meshO1.cell2node.father->Size());
226 cell2cellOrig.father->Resize(meshO1.cell2cellOrig.father->Size());
227 for (index iCell = 0; iCell < meshO1.cell2node.father->Size(); iCell++)
228 {
229 Elem::Element eCell = meshO1.GetCellElement(iCell);
230 Elem::Element eCellO2 = eCell.ObtainElevatedElem();
231 cell2node.father->ResizeRow(iCell, eCellO2.GetNumNodes());
232 if (isPeriodic)
233 cell2nodePbi.father->ResizeRow(iCell, eCellO2.GetNumNodes());
234 auto c2n = meshO1.cell2node[iCell];
235 for (int ic2n = 0; ic2n < c2n.size(); ic2n++)
236 {
237 // fill in the O1 nodes, note that global indices for O1 nodes have changed
238 index iNodeOldGlobal = meshO1.coords.trans.pLGhostMapping->operator()(-1, c2n[ic2n]);
239 index nodeOldOrigLocalIdx{-1};
240 int nodeOldOrigRank{-1};
241 if (!meshO1.coords.trans.pLGlobalMapping->search(iNodeOldGlobal, nodeOldOrigRank, nodeOldOrigLocalIdx))
242 DNDS_assert_info(false, "search failed");
243 // nodeOldOrigRank and nodeOldOrigLocalIdx is same in new
244 cell2node(iCell, ic2n) =
245 coords.father->pLGlobalMapping->operator()(nodeOldOrigRank, nodeOldOrigLocalIdx); // now point to global
246
247 // fill in node pbi
248 if (isPeriodic)
249 cell2nodePbi(iCell, ic2n) = meshO1.cell2nodePbi(iCell, ic2n);
250 }
251 for (int ic2n = c2n.size(); ic2n < cell2node[iCell].size(); ic2n++)
252 cell2node(iCell, ic2n) = -1;
253
254 cellElemInfo(iCell, 0) = meshO1.cellElemInfo(iCell, 0);
255 cellElemInfo(iCell, 0).setElemType(eCellO2.type); // update cell elem info
256 cell2cellOrig[iCell] = meshO1.cell2cellOrig[iCell];
257 }
258 // std::cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXX" << std::endl;
259 for (index iCell = 0; iCell < meshO1.cell2node.father->Size(); iCell++)
260 {
261 Elem::Element eCell = meshO1.GetCellElement(iCell);
262 auto c2n = meshO1.cell2node[iCell];
263 auto c2c = meshO1.cell2cell[iCell];
264 index iCellGlob = meshO1.cell2node.trans.pLGhostMapping->operator()(-1, iCell);
265
266 // std::vector<index> cellNewNodeRowV;
267 // cellNewNodeRowV.resize(eCell.GetNumElev_O1O2());
268 SmallCoordsAsVector coordsC;
269 meshO1.GetCoordsOnCell(iCell, coordsC);
270
271 for (int iNNode = 0; iNNode < eCell.GetNumElev_O1O2(); iNNode++)
272 {
273 auto eSpan = eCell.ObtainElevNodeSpan(iNNode);
274 std::array<index, Elem::CellNumNodeMax> spanNodes{};
275 eCell.ExtractElevNodeSpanNodes(iNNode, c2n, spanNodes);
276 auto spanNodesSrt = spanNodes;
277 std::sort(spanNodesSrt.begin(), spanNodesSrt.begin() + eSpan.GetNumNodes());
278 tPoint newNodeCoord;
279 { // the new node method, should cope with shape func
280 SmallCoordsAsVector coordsSpan;
281 coordsSpan.resize(3, eSpan.GetNumNodes());
282 eCell.ExtractElevNodeSpanNodes(iNNode, coordsC, coordsSpan);
283 newNodeCoord = coordsSpan.rowwise().mean();
284
285 // std::cout << fmt::format("############\niCell {}, iNNode {}", iCell, iNNode) << std::endl;
286 // std::cout << newNodeCoord.transpose() << std::endl;
287
288 if (isPeriodic)
289 {
290 std::vector<NodePeriodicBits> pbi;
291 pbi.resize(eSpan.GetNumNodes());
292 eCell.ExtractElevNodeSpanNodes(iNNode, meshO1.cell2nodePbi[iCell], pbi);
293 NodePeriodicBits pbiR;
294 pbiR.setP1True(), pbiR.setP2True(), pbiR.setP3True();
295 for (auto pbie : pbi)
296 pbiR = pbiR & pbie;
297 // std::cout << uint(uint8_t(pbiR)) << std::endl;
298 newNodeCoord = periodicInfo.GetCoordBackByBits(newNodeCoord, pbiR);
299 cell2nodePbi[iCell][c2n.size() + iNNode] = pbiR; //* fill In the O2 cell2nodePbi
300 // ** cell2nodePbiO2: edge using edge common, face using face common
301 // ** note that this technique requires at least 2 layer each periodic direction
302 }
303 }
304 if (cellNewNodes[iCell][iNNode] >= 0)
305 {
306 cell2node[iCell][c2n.size() + iNNode] = cellNewNodes[iCell][iNNode];
307 continue;
308 }
309
310 int nFound = 0;
311 real minSqrDist = veryLargeReal;
312 for (auto iNNodeC : cellNewNodes[c2c[-1 - cellNewNodes[iCell][iNNode]]])
313 {
314 if (iNNodeC < 0)
315 continue;
316 // iNNodeC is global iNode
317 MPI_int rank{-1};
318 index val{-1};
319 if (!coords.trans.pLGhostMapping->search_indexAppend(iNNodeC, rank, val))
320 DNDS_assert_info(false, "search failed");
321 // std::cout << fmt::format("val {}, iNNodeC {}", val, iNNodeC) << std::endl;
322 // std::cout << (coords[val] - newNodeCoord).squaredNorm() << std::endl;
323 // std::cout << coords[val].transpose() << ", " << newNodeCoord.transpose() << std::endl;
324 real sqrDist = (coords[val] - newNodeCoord).squaredNorm();
325 minSqrDist = std::min(sqrDist, minSqrDist);
326 if (sqrDist < epsSqrDist)
327 {
328 nFound++;
329 cell2node[iCell][c2n.size() + iNNode] = iNNodeC;
330 }
331 }
332 DNDS_assert_info(nFound == 1, fmt::format("geometric search for elevated point failed, nFound = {}, minSqrDist = {}, pos ({},{},{})", nFound, minSqrDist,
333 newNodeCoord(0), newNodeCoord(1), newNodeCoord(2)));
334 //* comment: way of ridding of the geometric search:
335 //* use interpolated topo: face/edge
336 //* or record the full vertex string for each O2 nodes' span (O2 nodes do not have siblings in the same span)
337 }
338 }
339 /**********************************/
340 //! now cell2cell, bnd2cell are lost
341 {
342 //! jumped because cell2cell, bnd2cell is to be recovered mandatorily
343 }
344 bnd2node.InitPair("BuildO2FromO1Elevation::bnd2node", mpi);
345 bndElemInfo.InitPair("BuildO2FromO1Elevation::bndElemInfo", mpi);
346 if (isPeriodic)
347 {
348 bnd2nodePbi.InitPair("BuildO2FromO1Elevation::bnd2nodePbi", getMPI());
349 }
350 bnd2bndOrig.InitPair("BuildO2FromO1Elevation::bnd2bndOrig", mpi);
351 bnd2node.father->Resize(meshO1.bnd2node.father->Size());
352 bndElemInfo.father->Resize(meshO1.bnd2node.father->Size());
353 if (isPeriodic)
354 bnd2nodePbi.father->Resize(meshO1.bnd2node.father->Size());
355 bnd2bndOrig.father->Resize(meshO1.bnd2node.father->Size());
356 for (index iBnd = 0; iBnd < meshO1.bnd2node.father->Size(); iBnd++)
357 {
358 index iCell = meshO1.bnd2cell(iBnd, 0); // my own bnd2cell is to global!
359 auto eCellO2 = this->GetCellElement(iCell);
360 auto b2n = meshO1.bnd2node[iBnd];
361 auto c2nO2 = this->cell2node[iCell];
362 std::vector<NodeIndexPBI> b2nPbi = meshO1.getBnd2NodeIndexPbiRow(iBnd);
363 std::vector<NodeIndexPBI> c2nPbiO2 = this->getCell2NodeIndexPbiRow(iCell);
364
365 std::vector<index> b2nv = b2n;
366 std::vector<NodeIndexPBI> b2nPbiv = b2nPbi;
367 for (int ib2n = 0; ib2n < b2nv.size(); ib2n++)
368 {
369 index iN = b2nv[ib2n];
370 //* note that b2nv holds O1 nodes' old and local indices
371 index iNodeOldGlobal = meshO1.coords.trans.pLGhostMapping->operator()(-1, iN);
372 index nodeOldOrigLocalIdx{-1};
373 int nodeOldOrigRank{-1};
374 if (!meshO1.coords.trans.pLGlobalMapping->search(iNodeOldGlobal, nodeOldOrigRank, nodeOldOrigLocalIdx))
375 DNDS_assert_info(false, "search failed");
376 // nodeOldOrigRank and nodeOldOrigLocalIdx is same in new
377 iN = coords.father->pLGlobalMapping->operator()(nodeOldOrigRank, nodeOldOrigLocalIdx); // now point to global
378 b2nv[ib2n] = iN;
379 b2nPbiv[ib2n].i = iN; // also update this
380 }
381 std::vector<index> b2nvSorted = b2nv;
382 std::vector<NodeIndexPBI> b2nPbivSorted = b2nPbiv;
383 std::sort(b2nvSorted.begin(), b2nvSorted.end());
384 std::sort(b2nPbivSorted.begin(), b2nPbivSorted.end());
385
386 int nFound{0};
387 int c2fFound{-1};
388 std::vector<index> b2nO2Found;
389 std::vector<NodeIndexPBI> b2nPbiO2Found;
390 Elem::Element eBndFound;
391 for (int ic2f = 0; ic2f < eCellO2.GetNumFaces(); ic2f++)
392 {
393 auto eFaceO2 = eCellO2.ObtainFace(ic2f);
394 std::vector<index> f2nO2, f2nO2Sorted;
395 std::vector<NodeIndexPBI> f2nPbiO2, f2nPbiO2Sorted;
396 f2nO2.resize(eFaceO2.GetNumNodes());
397 eCellO2.ExtractFaceNodes(ic2f, c2nO2, f2nO2);
398 if (isPeriodic)
399 {
400 f2nPbiO2.resize(eFaceO2.GetNumNodes());
401 eCellO2.ExtractFaceNodes(ic2f, c2nPbiO2, f2nPbiO2);
402 }
403 f2nO2Sorted = f2nO2;
404 std::sort(f2nO2Sorted.begin(), f2nO2Sorted.end()); //! cannot use sorted
405 f2nPbiO2Sorted = f2nPbiO2;
406 std::sort(f2nPbiO2Sorted.begin(), f2nPbiO2Sorted.end()); //! cannot use sorted
407 if (std::includes(f2nO2Sorted.begin(), f2nO2Sorted.end(), b2nvSorted.begin(), b2nvSorted.end()))
408 {
409 if (isPeriodic) // need to doublecheck if pbis match
410 {
411 if (std::includes(f2nPbiO2Sorted.begin(), f2nPbiO2Sorted.end(), b2nPbivSorted.begin(), b2nPbivSorted.end()))
412 ;
413 else
414 continue;
415 }
416 nFound++;
417 c2fFound = ic2f;
418 b2nO2Found = f2nO2;
419 b2nPbiO2Found = f2nPbiO2;
420 eBndFound = eFaceO2;
421 }
422 }
423 DNDS_assert(nFound == 1);
424 bnd2node.father->ResizeRow(iBnd, b2nO2Found.size());
425 bnd2node[iBnd] = b2nO2Found;
426 if (isPeriodic)
427 {
428 bnd2nodePbi.father->ResizeRow(iBnd, b2nO2Found.size());
429 for (int ib2n = 0; ib2n < b2nPbiO2Found.size(); ib2n++)
430 bnd2nodePbi[iBnd][ib2n] = b2nPbiO2Found[ib2n].pbi;
431 }
432 bndElemInfo(iBnd, 0) = meshO1.bndElemInfo(iBnd, 0);
433 bndElemInfo(iBnd, 0).setElemType(eBndFound.type);
434 bnd2bndOrig[iBnd] = meshO1.bnd2bndOrig[iBnd];
435 }
436 adjPrimaryState = Adj_PointToGlobal;
437 cell2node.idx.markGlobal();
438 bnd2node.idx.markGlobal();
439
440 coords.son = make_ssp<decltype(coords.son)::element_type>(ObjName{"BuildO2FromO1Elevation::coords.son"}, mpi); // delete because reconstructed later
441
442 // this->BuildGhostPrimary();
443
444 // this->AdjGlobal2LocalPrimary();
445 // if (mpi.rank == 0)
446 // {
447 // for (index iCell = 0; iCell < meshO1.cell2node.Size(); iCell++)
448 // {
449 // std::vector<index> v1 = meshO1.cell2node[iCell];
450 // std::vector<index> v2 = cell2node[iCell];
451 // std::cout << "iCell " << iCell << std::endl;
452 // for (auto i : v1)
453 // std::cout << i << "(" << meshO1.coords[i].transpose() << ")"
454 // << ", ";
455 // std::cout << std::endl;
456 // for (auto i : v2)
457 // std::cout << i << "(" << coords[i].transpose() << ")"
458 // << ", ";
459 // std::cout << std::endl;
460
461 // std::vector<index> v3 = meshO1.cell2cell[iCell];
462 // std::vector<index> v4 = cell2cell[iCell];
463 // for (auto i : v3)
464 // std::cout << i
465 // << ", ";
466 // std::cout << std::endl;
467 // for (auto i : v4)
468 // std::cout << i
469 // << ", ";
470 // std::cout << std::endl;
471 // }
472 // }
473 // // this->AdjLocal2GlobalPrimary();
474 }
475
476 void UnstructuredMesh::BuildBisectO1FormO2(UnstructuredMesh &meshO2)
477 {
478 bool O2MeshIsO2 = meshO2.IsO2();
479 DNDS_assert(O2MeshIsO2);
480 DNDS_assert(meshO2.adjPrimaryState == Adj_PointToGlobal);
481 DNDS_assert(meshO2.cell2node.isGlobal() && meshO2.bnd2node.isGlobal());
482 adjPrimaryState = meshO2.adjPrimaryState;
483 cell2node.idx.markGlobal();
484 bnd2node.idx.markGlobal();
485 mRank = meshO2.mRank;
486 mpi = meshO2.mpi;
487 dim = meshO2.dim;
488 isPeriodic = meshO2.isPeriodic;
489 periodicInfo = meshO2.periodicInfo;
490 nNodeO1 = meshO2.nNodeO1;
491 elevState = MeshElevationState::Elevation_Untouched;
492
493 coords.son = make_ssp<decltype(coords.son)::element_type>(ObjName{"BuildBisectO1FormO2::coords.son"}, mpi);
494 coords.father = meshO2.coords.father; // coords are the same, taken without change
495 node2nodeOrig.son = make_ssp<decltype(node2nodeOrig.son)::element_type>(ObjName{"BuildBisectO1FormO2::node2nodeOrig.son"}, mpi);
496 node2nodeOrig.father = meshO2.node2nodeOrig.father; // node2nodeOrig corresponds to coords, and is taken without change
497
498 /**********************************/
499 //* cell
500 std::vector<std::vector<index>> cell2nodeV;
501 std::vector<std::vector<NodePeriodicBits>> cell2nodePbiV;
502 std::vector<ElemInfo> cellElemInfoV;
503 // std::vector<index> cell2cellOrigV;
504
505 index newNumCell = 0;
506 for (index iCell = 0; iCell < meshO2.cell2node.father->Size(); iCell++)
507 {
508 auto eCell = meshO2.GetCellElement(iCell);
509 newNumCell += eCell.GetO2NumBisect();
510 }
511 cell2nodeV.reserve(newNumCell);
512 cell2nodePbiV.reserve(newNumCell);
513 cellElemInfoV.reserve(newNumCell);
514 // cell2cellOrigV.reserve(newNumCell);
515
516 // index newNumCellOffset = 0;
517 // MPI::Scan(&newNumCell, &newNumCellOffset, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
518 // newNumCellOffset -= newNumCell;
519
520 for (index iCell = 0; iCell < meshO2.cell2node.father->Size(); iCell++)
521 {
522 SmallCoordsAsVector coordsC;
523 auto eCell = meshO2.GetCellElement(iCell);
524 meshO2.GetCoordsOnCell(iCell, coordsC);
525 int nBi = eCell.GetO2NumBisect();
526 int iBiVariant = GetO2ElemBisectVariant(eCell, coordsC);
527 for (int iBi = 0; iBi < nBi; iBi++)
528 {
529 auto eCellSub = eCell.ObtainO2BisectElem(iBi);
530 std::vector<index> c2nSub;
531 c2nSub.resize(eCellSub.GetNumNodes());
532 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, meshO2.cell2node[iCell], c2nSub);
533 cell2nodeV.push_back(c2nSub);
534 ElemInfo eInfo = meshO2.cellElemInfo(iCell, 0);
535 eInfo.setElemType(eCellSub.type);
536 cellElemInfoV.push_back(eInfo);
537 if (isPeriodic)
538 {
539 std::vector<NodePeriodicBits> c2nPbiSub;
540 c2nPbiSub.resize(eCellSub.GetNumNodes());
541 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, meshO2.cell2nodePbi[iCell], c2nPbiSub);
542 cell2nodePbiV.push_back(c2nPbiSub);
543 }
544 //! cell2cellOrig uses new global size later, cell2cellOrigV discarded
545 // cell2cellOrigV.push_back(newNumCellOffset + cell2cellOrigV.size());
546 // cell2cellOrigV.push_back(UnInitIndex);
547 }
548 }
549 // std::cout << "here1" << std::endl;
550 cell2node.InitPair("BuildBisectO1FormO2::cell2node", mpi);
551 cellElemInfo.InitPair("BuildBisectO1FormO2::cellElemInfo", mpi);
552 if (isPeriodic)
553 { // we preserve a sample of detailed array constructor
554 cell2nodePbi.InitPair("BuildBisectO1FormO2::cell2nodePbi", NodePeriodicBits::CommType(), NodePeriodicBits::CommMult(), mpi);
555 }
556 cell2cellOrig.InitPair("BuildBisectO1FormO2::cell2cellOrig", mpi);
557 cell2node.father->Resize(cell2nodeV.size());
558 cellElemInfo.father->Resize(cellElemInfoV.size());
559 DNDS_assert(cellElemInfoV.size() == cell2nodeV.size());
560 if (isPeriodic)
561 cell2nodePbi.father->Resize(cell2nodePbiV.size()),
562 DNDS_assert(cell2nodePbiV.size() == cell2nodeV.size());
563 cell2cellOrig.father->Resize(cell2nodeV.size());
564 cell2node.father->createGlobalMapping();
565 for (index i = 0; i < cell2nodeV.size(); i++)
566 {
567 cell2node.father->ResizeRow(i, cell2nodeV[i].size());
568 for (rowsize ic2n = 0; ic2n < cell2nodeV[i].size(); ic2n++)
569 cell2node.father->operator()(i, ic2n) = cell2nodeV[i][ic2n];
570 cellElemInfo(i, 0) = cellElemInfoV[i];
571 if (isPeriodic)
572 {
573 DNDS_assert(cell2nodePbiV[i].size() == cell2nodeV[i].size());
574 cell2nodePbi.father->ResizeRow(i, cell2nodePbiV[i].size());
575 for (rowsize ic2n = 0; ic2n < cell2nodePbiV[i].size(); ic2n++)
576 cell2nodePbi.father->operator()(i, ic2n) = cell2nodePbiV[i][ic2n];
577 }
578 //! cell2cellOrig uses new global size now
579 cell2cellOrig(i, 0) = cell2node.father->pLGlobalMapping->operator()(mpi.rank, i);
580 }
581 // std::cout << "here2" << std::endl;
582
583 /**********************************/
584 //* bnd
585
586 std::vector<std::vector<index>> bnd2nodeV;
587 std::vector<std::vector<NodePeriodicBits>> bnd2nodePbiV;
588 std::vector<ElemInfo> bndElemInfoV;
589 for (index iBnd = 0; iBnd < meshO2.bnd2node.father->Size(); iBnd++)
590 {
591 auto eBnd = meshO2.GetBndElement(iBnd);
592
593 index iCellG = meshO2.bnd2cell(iBnd, 0);
594 index iCell = meshO2.CellIndexGlobal2Local_NoSon(iCellG);
595 // if (mpi.rank == 0)
596 // std::cout << "iCell " << iCell << "; " << meshO2.cell2node.father->Size() << std::endl;
597 DNDS_assert_info(iCell >= 0, "bnd's main cell is not in current process main, need reordering of bnd");
598 SmallCoordsAsVector coordsC;
599 meshO2.GetCoordsOnCell(iCell, coordsC);
600 auto c2n = meshO2.cell2node[iCell];
601 SmallCoordsAsVector coordsB;
602 coordsB.resize(Eigen::NoChange, meshO2.bnd2node.RowSize(iBnd));
603 for (rowsize ib2n = 0; ib2n < meshO2.bnd2node[iBnd].size(); ib2n++)
604 {
605 bool found{false};
606 for (rowsize ic2n = 0; ic2n < coordsC.cols(); ic2n++)
607 if (c2n[ic2n] == meshO2.bnd2node[iBnd][ib2n])
608 coordsB[ib2n] = coordsC[ic2n], found = true;
609 DNDS_assert(found);
610 }
611
612 int nBi = eBnd.GetO2NumBisect();
613 int iBiVariant = GetO2ElemBisectVariant(eBnd, coordsB);
614 for (int iBi = 0; iBi < nBi; iBi++)
615 {
616 auto eBndSub = eBnd.ObtainO2BisectElem(iBi);
617 std::vector<index> b2nSub;
618 b2nSub.resize(eBndSub.GetNumNodes());
619 std::vector<NodePeriodicBits> b2nPbiSub;
620 b2nPbiSub.resize(eBndSub.GetNumNodes());
621 eBnd.ExtractO2BisectElemNodes(iBi, iBiVariant, meshO2.bnd2node[iBnd], b2nSub);
622 bnd2nodeV.push_back(b2nSub);
623 if (isPeriodic)
624 {
625 eBnd.ExtractO2BisectElemNodes(iBi, iBiVariant, meshO2.bnd2nodePbi[iBnd], b2nPbiSub);
626 bnd2nodePbiV.push_back(b2nPbiSub);
627 }
628 ElemInfo eInfo = meshO2.bndElemInfo(iBnd, 0);
629 eInfo.setElemType(eBndSub.type);
630 bndElemInfoV.push_back(eInfo);
631 }
632 }
633 bnd2node.InitPair("BuildBisectO1FormO2::bnd2node", mpi);
634 bndElemInfo.InitPair("BuildBisectO1FormO2::bndElemInfo", ElemInfo::CommType(), ElemInfo::CommMult(), mpi);
635 if (isPeriodic)
636 {
637 bnd2nodePbi.InitPair("BuildBisectO1FormO2::bnd2nodePbi", mpi);
638 }
639 bnd2bndOrig.InitPair("BuildBisectO1FormO2::bnd2bndOrig", mpi);
640 bnd2node.father->Resize(bnd2nodeV.size());
641 bndElemInfo.father->Resize(bndElemInfoV.size());
642 DNDS_assert(bndElemInfoV.size() == bnd2nodeV.size());
643 if (isPeriodic)
644 bnd2nodePbi.father->Resize(bnd2nodePbiV.size()), DNDS_assert(bnd2nodePbiV.size() == bnd2nodeV.size());
645 bnd2bndOrig.father->Resize(bnd2nodeV.size());
646 bnd2node.father->createGlobalMapping();
647
648 for (index i = 0; i < bnd2nodeV.size(); i++)
649 {
650 bnd2node.father->ResizeRow(i, bnd2nodeV[i].size());
651 for (rowsize ic2n = 0; ic2n < bnd2nodeV[i].size(); ic2n++)
652 bnd2node.father->operator()(i, ic2n) = bnd2nodeV[i][ic2n];
653 if (isPeriodic)
654 {
655 bnd2nodePbi.father->ResizeRow(i, bnd2nodePbiV[i].size());
656 DNDS_assert(bnd2nodePbiV[i].size() == bnd2nodeV[i].size());
657 for (rowsize ic2n = 0; ic2n < bnd2nodePbiV[i].size(); ic2n++)
658 bnd2nodePbi.father->operator()(i, ic2n) = bnd2nodePbiV[i][ic2n];
659 }
660 bndElemInfo(i, 0) = bndElemInfoV[i];
661 //! bnd2bndOrig uses new global size now
662 bnd2bndOrig(i, 0) = bnd2node.father->pLGlobalMapping->operator()(mpi.rank, i);
663 }
664
665 //* unhandled info
666 bnd2cell.InitPair("BuildBisectO1FormO2::bnd2cell", mpi);
667 cell2cell.InitPair("BuildBisectO1FormO2::cell2cell", mpi);
668 }
669
670 static tPoint HermiteInterpolateMidPointOnLine2WithNorm(tPoint c0, tPoint c1, tPoint n0, tPoint n1)
671 {
672 tPoint c01 = c1 - c0;
673 tPoint c01U = c01.stableNormalized();
674 real c01Len = c01.stableNorm();
675 tPoint t0 = -c01U.cross(n0).cross(n0);
676 tPoint t1 = -c01U.cross(n1).cross(n1);
677 t0 = t0.stableNormalized() * c01Len;
678 t1 = t1.stableNormalized() * c01Len;
679 if (n0.norm() == 0 && n1.norm() != 0)
680 return 0.25 * c0 + 0.75 * c1 - 0.25 * t1;
681 if (n0.norm() != 0 && n1.norm() == 0)
682 return 0.75 * c0 + 0.25 * c1 + 0.25 * t0;
683 if (n0.norm() == 0 && n1.norm() == 0)
684 return 0.5 * (c0 + c1);
685 return 0.5 * (c0 + c1) + 0.125 * (t0 - t1);
686 }
687
688 static tPoint HermiteInterpolateMidPointOnQuad4WithNorm(
689 tPoint c0, tPoint c1, tPoint c2, tPoint c3,
690 tPoint n0, tPoint n1, tPoint n2, tPoint n3)
691 {
692 tPoint c01 = HermiteInterpolateMidPointOnLine2WithNorm(c0, c1, n0, n1);
693 tPoint c12 = HermiteInterpolateMidPointOnLine2WithNorm(c1, c2, n1, n2);
694 tPoint c23 = HermiteInterpolateMidPointOnLine2WithNorm(c2, c3, n2, n3);
695 tPoint c30 = HermiteInterpolateMidPointOnLine2WithNorm(c3, c0, n3, n0);
696 return 0.5 * (c01 + c12 + c23 + c30) - 0.25 * (c0 + c1 + c2 + c3);
697 }
698
699 void UnstructuredMesh::ElevatedNodesGetBoundarySmooth(const std::function<bool(t_index)> &FiFBndIdNeedSmooth)
700 {
701 // if (mpi.rank == 1)
702 // for (index iCell = 0; iCell < cell2face.Size(); iCell++)
703 // {
704 // std::cout << iCell << ": ";
705 // for (auto i : cell2face[iCell])
706 // std::cout << i << ", ";
707 // std::cout << std::endl;
708 // }
709 DNDS_assert(elevState == Elevation_O1O2);
710 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
711 DNDS_assert(cell2node.isLocal() && bnd2node.isLocal());
712 DNDS_assert(adjFacialState == Adj_PointToLocal);
713 DNDS_assert(face2cell.isLocal() && face2node.isLocal());
714 DNDS_assert(adjC2FState == Adj_PointToLocal);
715 DNDS_assert(cell2face.isLocal() && bnd2face.isLocal());
716 DNDS_assert(face2node.father);
717
718 coordsElevDisp.InitPair("ElevatedNodesGetBoundarySmooth::coordsElevDisp", mpi);
719 coordsElevDisp.father->Resize(coords.father->Size());
720 coordsElevDisp.TransAttach();
721 coordsElevDisp.trans.BorrowGGIndexing(coords.trans);
722 coordsElevDisp.trans.createMPITypes();
723
724 for (index iN = 0; iN < coords.Size(); iN++)
725 coordsElevDisp[iN].setConstant(largeReal);
726
727 // tAdj1Pair nCoordBndNum;
728 // DNDS_MAKE_SSP(nCoordBndNum.father, mpi);
729 // DNDS_MAKE_SSP(nCoordBndNum.son, mpi);
730 // nCoordBndNum.father->Resize(coords.father->Size());
731 // nCoordBndNum.TransAttach();
732 // nCoordBndNum.trans.BorrowGGIndexing(coords.trans);
733 // nCoordBndNum.trans.createMPITypes();
734
735 // for (index iN = 0; iN < coords.size(); iN++)
736 // nCoordBndNum[iN][0] = 0;
737
738 /***********************************/
739 // build faceExteded info
740 tAdjPair face2nodeExtended;
741 tElemInfoArrayPair faceElemInfoExtended;
742 tPbiPair face2nodePbiExtended;
743 face2nodeExtended.father = face2node.father;
744 faceElemInfoExtended.father = faceElemInfo.father;
745 if (isPeriodic)
746 face2nodePbiExtended.father = face2nodePbi.father;
747 face2nodeExtended.son = make_ssp<decltype(face2nodeExtended.son)::element_type>(ObjName{"ElevatedNodesGetBoundarySmooth::face2nodeExtended.son"}, mpi);
748 faceElemInfoExtended.son = make_ssp<decltype(faceElemInfoExtended.son)::element_type>(ObjName{"ElevatedNodesGetBoundarySmooth::faceElemInfoExtended.son"}, ElemInfo::CommType(), ElemInfo::CommMult(), mpi);
749 if (isPeriodic)
750 face2nodePbiExtended.son = make_ssp<decltype(face2nodePbiExtended.son)::element_type>(ObjName{"ElevatedNodesGetBoundarySmooth::face2nodePbiExtended.son"}, NodePeriodicBits::CommType(), NodePeriodicBits::CommMult(), mpi);
751
752 this->AdjLocal2GlobalFacial();
753 this->AdjLocal2GlobalC2F();
754
755 std::vector<index> faceGhostExt;
756 for (index iCell = 0; iCell < cell2face.Size(); iCell++)
757 for (auto iFace : cell2face[iCell])
758 {
759 DNDS_assert_info(iFace >= 0 || iFace == UnInitIndex, fmt::format("iFace {}", iFace));
760 if (iFace == UnInitIndex) // old cell2face could contain void pointing
761 continue;
762 DNDS::MPI_int rank = UnInitMPIInt;
763 DNDS::index val = UnInitIndex;
764 if (!face2node.trans.pLGlobalMapping->search(iFace, rank, val))
765 DNDS_assert_info(false, "search failed");
766 if (rank != mpi.rank)
767 faceGhostExt.push_back(iFace);
768 // if (mpi.rank == 1)
769 // std::cout << "added Face " << iFace << std::endl;
770 }
771 face2nodeExtended.TransAttach();
772 faceElemInfoExtended.TransAttach();
773 if (isPeriodic)
774 face2nodePbiExtended.TransAttach();
775 face2nodeExtended.trans.createGhostMapping(faceGhostExt);
776 faceElemInfoExtended.trans.BorrowGGIndexing(face2nodeExtended.trans);
777 if (isPeriodic)
778 face2nodePbiExtended.trans.BorrowGGIndexing(face2nodeExtended.trans);
779 face2nodeExtended.trans.createMPITypes();
780 faceElemInfoExtended.trans.createMPITypes();
781 if (isPeriodic)
782 face2nodePbiExtended.trans.createMPITypes();
783 face2nodeExtended.trans.pullOnce();
784 faceElemInfoExtended.trans.pullOnce();
785 if (isPeriodic)
786 face2nodePbiExtended.trans.pullOnce();
787
788 // std::cout << fmt::format("rank {} faceElemInfo {} {}, face2node {} {}",
789 // mpi.rank,
790 // faceElemInfo.father->Size(), faceElemInfo.son->Size(),
791 // face2node.father->Size(), face2node.son->Size())
792 // << std::endl;
793 // std::cout << fmt::format("rank {} faceElemInfoExt {} {}, face2nodeExt {} {}",
794 // mpi.rank,
795 // faceElemInfoExtended.father->Size(), faceElemInfoExtended.son->Size(),
796 // face2nodeExtended.father->Size(), face2nodeExtended.son->Size())
797 // << std::endl;
798
799 this->AdjGlobal2LocalFacial();
800 this->AdjGlobal2LocalC2F();
801 //** a direct copy from AdjGlobal2LocalPrimary()
802 auto NodeIndexGlobal2Local = [&](DNDS::index &iNodeOther)
803 {
804 if (iNodeOther == UnInitIndex)
805 return;
806 DNDS::MPI_int rank = UnInitMPIInt;
807 DNDS::index val = UnInitIndex;
808 // if (!cell2cell.trans.pLGlobalMapping->search(iCellOther, rank, val))
809 // DNDS_assert_info(false, "search failed");
810 // if (rank != mpi.rank)
811 // iCellOther = -1 - iCellOther;
812 auto result = coords.trans.pLGhostMapping->search_indexAppend(iNodeOther, rank, val);
813 if (result)
814 iNodeOther = val;
815 else
816 iNodeOther = -1 - iNodeOther; // mapping to un-found in father-son
817 };
818
819 // std::cout << fmt::format("rank {} Sizes ext : extFaceProc: {}->{}", mpi.rank, faceElemInfo.Size(), faceElemInfoExtended.Size()) << std::endl;
820 // if (mpi.rank == 1)
821 // for (index iFace = 0; iFace < face2nodeExtended.Size(); iFace++)
822 // {
823 // std::cout << fmt::format("rank {}, face {}", mpi.rank, face2nodeExtended.trans.pLGhostMapping->operator()(-1, iFace)) << std::endl;
824 // }
825
826 // this->AdjGlobal2LocalPrimary();
827
828 for (index iFace = face2nodeExtended.father->Size(); iFace < face2nodeExtended.Size(); iFace++)
829 {
830 for (auto &iN : face2nodeExtended[iFace])
831 {
832 DNDS_assert_info(iN >= 0, fmt::format("rank {}, iN {}", mpi.rank, iN)); // can't be unfound node
833 NodeIndexGlobal2Local(iN);
834 DNDS_assert_info(iN >= 0, fmt::format("rank {}, iN {}", mpi.rank, iN)); // can't be unfound node
835 }
836 }
837 // std::cout << mpi.rank << " rank Here XXXX 1" << std::endl;
838 /***********************************/
839 // build nodeNormClusters
840 using t3VecsPair = ArrayPair<ArrayEigenUniMatrixBatch<3, 1>>;
841 t3VecsPair nodeNormClusters;
842 nodeNormClusters.InitPair("ElevatedNodesGetBoundarySmooth::nodeNormClusters", mpi);
843 nodeNormClusters.father->Resize(coords.father->Size(), 3, 1);
844
845 std::vector<int> nodeBndNum(coords.father->Size(), 0); //? need row-appending methods for NonUniform arrays?
846 for (index iFace = 0; iFace < face2nodeExtended.Size(); iFace++)
847 {
848 auto faceBndID = faceElemInfoExtended(iFace, 0).zone;
849 if (!FiFBndIdNeedSmooth(faceBndID))
850 continue;
851 for (auto iNode : face2nodeExtended[iFace])
852 if (iNode < coords.father->Size() && iNode < nNodeO1) //* being local node and O1 node
853 nodeBndNum.at(iNode)++;
854 }
855 // if (mpi.rank == 0)
856 // {
857 // // std::cout << "faceBndID: \n";
858 // // for (index iFace = 0; iFace < face2nodeExtended.Size(); iFace++)
859 // // {
860 // // auto faceBndID = faceElemInfoExtended(iFace, 0).zone;
861 // // std::cout << faceBndID << ", " << FiFBndIdNeedSmooth(faceBndID) << "; "
862 // // << coords[face2node(iFace, 0)].transpose()
863 // // << " ||| " << coords[face2node(iFace, 1)].transpose() << std::endl;
864 // // }
865 // // std::cout << "XXXXXXXXXXXX" << std::endl;
866 // std::cout << "nodeBndNum: ";
867 // for (auto i : nodeBndNum)
868 // std::cout
869 // << i << ";";
870 // std::cout << "XXXXXXXXXXXX" << std::endl;
871 // }
872
873 // std::cout << mpi.rank << " rank Here XXXX 2" << std::endl;
874 for (index iNode = 0; iNode < coords.father->Size(); iNode++)
875 nodeNormClusters.father->ResizeBatch(iNode, std::max(nodeBndNum.at(iNode), 0));
876 for (auto &v : nodeBndNum)
877 v = 0; // set to zero
878 // std::cout << mpi.rank << " rank Here XXXX 3" << std::endl;
879 auto GetCoordsOnFaceExtended = [&](index iFace, tSmallCoords &cs)
880 {
881 if (!isPeriodic)
882 _detail_GetCoords(face2nodeExtended[iFace], cs);
883 else
884 _detail_GetCoordsOnElem(face2nodeExtended[iFace], face2nodePbiExtended[iFace], cs);
885 };
886
887 for (index iFace = 0; iFace < face2nodeExtended.Size(); iFace++)
888 {
889 auto faceBndID = faceElemInfoExtended(iFace, 0).zone;
890 if (!FiFBndIdNeedSmooth(faceBndID))
891 continue;
892 auto eFace = Elem::Element{faceElemInfoExtended(iFace, 0).getElemType()}; // O1 faces could use only one norm (not strict for Quad4)
893 auto qFace = Elem::Quadrature{eFace, 1};
894 real faceArea{0};
895 tPoint uNorm;
896 tSmallCoords coordsF;
897 GetCoordsOnFaceExtended(iFace, coordsF);
898 qFace.Integration(
899 faceArea,
900 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
901 {
902 tJacobi J = Elem::ShapeJacobianCoordD01Nj(coordsF, DiNj);
903 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coordsF, DiNj);
904 tPoint np;
905 if (dim == 2)
906 np = FacialJacobianToNormVec<2>(J);
907 else
908 np = FacialJacobianToNormVec<3>(J);
909 // np.stableNormalize(); // do not normalize to preserve face area
910 uNorm = np;
911 vInc = np.norm();
912 });
913
914 for (int if2n = 0; if2n < face2nodeExtended[iFace].size(); if2n++)
915 {
916 index iNode = face2nodeExtended[iFace][if2n];
917 if (iNode < coords.father->Size() && iNode < nNodeO1) //* being local node and O1 node
918 {
919 tPoint uNormC = uNorm;
920 if (isPeriodic)
921 uNormC = periodicInfo.GetVectorBackByBits<3, 1>(uNorm, face2nodePbiExtended[iFace][if2n]);
922 nodeNormClusters(iNode, nodeBndNum.at(iNode)++) = uNormC;
923 }
924 }
925 }
926 nodeNormClusters.TransAttach();
927 nodeNormClusters.trans.BorrowGGIndexing(coords.trans);
928 nodeNormClusters.trans.createMPITypes();
929 MPI::Barrier(mpi.comm);
930 nodeNormClusters.trans.pullOnce();
931
932 /***********************************/
933 // do interpolation
934 tPoint bboxMove;
935 bboxMove.setZero();
936 index nMoved{0};
937 std::unordered_set<index> moveds;
938 for (index iFace = 0; iFace < face2nodeExtended.Size(); iFace++)
939 {
940 auto faceBndID = faceElemInfoExtended(iFace, 0).zone;
941 // TODO: add marking for non-moving bounaries
942 if (!FiFBndIdNeedSmooth(faceBndID))
943 {
944 if (!FaceIDIsInternal(faceBndID)) // * note don't operate on internal faces
945 for (auto iN : face2nodeExtended[iFace])
946 if (coordsElevDisp[iN](0) == largeReal && coordsElevDisp[iN](2) != 2 * largeReal)
947 coordsElevDisp[iN](2) = 3 * largeReal; // marking other boundary nodes
948 continue;
949 }
950 // for (auto iN : face2nodeExtended[iFace])
951 // coordsElevDisp[iN](2) = 2 * largeReal; // marking O1 nodes
952 auto eFace = Elem::Element{faceElemInfoExtended(iFace, 0).getElemType()}; // O1 faces could use only one norm (not strict for Quad4)
953 auto qFace = Elem::Quadrature{eFace, 1};
954 Elem::SummationNoOp noOp;
955 tPoint uNorm;
956 SmallCoordsAsVector coordsF;
957 GetCoordsOnFaceExtended(iFace, coordsF);
958 qFace.Integration(
959 noOp,
960 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
961 {
962 tJacobi J = Elem::ShapeJacobianCoordD01Nj(coordsF, DiNj);
963 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coordsF, DiNj);
964 tPoint np;
965 if (dim == 2)
966 np = FacialJacobianToNormVec<2>(J);
967 else
968 np = FacialJacobianToNormVec<3>(J);
969 np.stableNormalize();
970 uNorm = np;
971 });
972
973 auto eFaceO1 = eFace.ObtainO1Elem();
974
975 auto f2n = face2nodeExtended[iFace];
976
977 std::vector<int> f2nVecSeq(f2n.size());
978 for (int i = 0; i < f2n.size(); i++)
979 f2nVecSeq[i] = i;
980 for (int if2n = 0; if2n < eFaceO1.GetNumNodes(); if2n++)
981 {
982 index iNode = f2n[if2n];
983 coordsElevDisp[iNode](2) = 2 * largeReal;
984 }
985
986 for (int if2n = eFaceO1.GetNumNodes(); if2n < f2n.size(); if2n++)
987 {
988 index iNode = f2n[if2n];
989 if (iNode >= coords.father->Size())
990 continue;
991 DNDS_assert(iNode >= nNodeO1);
992
993 std::vector<index> spanO2Node;
994 std::vector<int> spanO2if2n;
995 SmallCoordsAsVector coordsSpan;
996 int iElev = if2n - eFaceO1.GetNumNodes();
997 int spanNnode = eFaceO1.ObtainElevNodeSpan(iElev).GetNumNodes();
998 spanO2Node.resize(spanNnode);
999 coordsSpan.resize(3, spanNnode);
1000 spanO2if2n.resize(spanNnode);
1001
1002 eFaceO1.ExtractElevNodeSpanNodes(iElev, f2n, spanO2Node); // should use O1f2n, but O2f2n is equivalent
1003 eFaceO1.ExtractElevNodeSpanNodes(iElev, coordsF, coordsSpan);
1004 eFaceO1.ExtractElevNodeSpanNodes(iElev, f2nVecSeq, spanO2if2n);
1005 tPoint cooUnsmooth = coordsF(EigenAll, if2n);
1006
1007 std::vector<tPoint> edges;
1008 if (spanO2Node.size() == 2)
1009 {
1010 edges.push_back((coordsSpan[1] - coordsSpan[0]).stableNormalized());
1011 }
1012 else if (spanO2Node.size() == 4)
1013 {
1014 edges.push_back((coordsSpan[1] - coordsSpan[0]).stableNormalized());
1015 edges.push_back((coordsSpan[2] - coordsSpan[1]).stableNormalized());
1016 edges.push_back((coordsSpan[3] - coordsSpan[2]).stableNormalized());
1017 edges.push_back((coordsSpan[0] - coordsSpan[3]).stableNormalized());
1018 }
1019 else
1020 DNDS_assert(false);
1021
1022 Eigen::Matrix<real, 3, 4> norms;
1023 Eigen::Vector<real, 4> nAdd;
1024 nAdd.setZero();
1025 norms.setZero();
1026 // if (mpi.rank == 1)
1027 // std::cout << nodeNormClusters.RowSize(spanO2Node[0]) << " nodeNormClustersRowSize" << std::endl;
1028
1029 for (int iN = 0; iN < spanO2Node.size(); iN++)
1030 {
1031 real sinValFFMax{0};
1032 for (int iNorm = 0; iNorm < nodeNormClusters.RowSize(spanO2Node[iN]); iNorm++)
1033 for (int jNorm = 0; jNorm < nodeNormClusters.RowSize(spanO2Node[iN]); jNorm++)
1034 sinValFFMax = std::max(
1035 sinValFFMax,
1036 std::sqrt(1 - sqr(std::abs(nodeNormClusters(spanO2Node[iN], iNorm)
1037 .stableNormalized()
1038 .dot(nodeNormClusters(spanO2Node[iN], jNorm).stableNormalized())))));
1039 for (int iNorm = 0; iNorm < nodeNormClusters.RowSize(spanO2Node[iN]); iNorm++)
1040 {
1041 tPoint normN = nodeNormClusters(spanO2Node[iN], iNorm);
1042 if (isPeriodic)
1043 normN = periodicInfo.GetVectorByBits<3, 1>(normN, face2nodePbiExtended[iFace][if2n]);
1044 real sinValMax = 0;
1045 for (auto &e : edges)
1046 sinValMax = std::max(sinValMax, std::abs(e.dot(normN.stableNormalized())));
1047 //! angle is here
1048 //* now using FFMaxAngle also
1049 if (std::max(sinValMax, sinValFFMax * 0) > std::sin(pi / 180. * elevationInfo.MaxIncludedAngle))
1050 continue;
1051 nAdd[iN] += 1.;
1052 // norms(EigenAll, iN) += normN * 1.;
1053 norms(EigenAll, iN) += normN.stableNormalized();
1054
1055 // std::cout << fmt::format("iFace {}, if2n {}, iN {}, iNorm{}, sinValMax {}; ====",
1056 // iFace, if2n, iN, iNorm, sinValMax) << normN.transpose() << std::endl;
1057 }
1058 }
1059 // norms.array().rowwise() /= nAdd.transpose().array();
1060 norms.colwise().normalize();
1061 for (int iN = 0; iN < spanO2Node.size(); iN++)
1062 {
1063 if (nAdd[iN] < 0.1) // no found, then no mov
1064 norms(EigenAll, iN).setZero();
1065 // DNDS_assert(nAdd[iN] > 0.1);
1066 }
1067
1068 tPoint cooSmooth;
1069 tPoint cooInc;
1070 if (spanO2Node.size() == 2)
1071 {
1072 cooSmooth = HermiteInterpolateMidPointOnLine2WithNorm(
1073 coordsSpan[0], coordsSpan[1], norms(EigenAll, 0), norms(EigenAll, 1));
1074 cooInc = cooSmooth - cooUnsmooth;
1075 }
1076 else // definitely is spanO2Node.size() == 4
1077 {
1078 f2n[if2n - 1];
1079 DNDS_assert(f2n.size() == 9); // has to be a Quad9
1080 cooInc = 0.5 * (coordsElevDisp[f2n[if2n - 4]] +
1081 coordsElevDisp[f2n[if2n - 3]] +
1082 coordsElevDisp[f2n[if2n - 2]] +
1083 coordsElevDisp[f2n[if2n - 1]]);
1084 if (isPeriodic)
1085 {
1086 cooInc =
1087 0.5 *
1088 (periodicInfo.GetVectorByBits<3, 1>(coordsElevDisp[f2n[if2n - 4]], face2nodePbiExtended[iFace][if2n - 4]) +
1089 periodicInfo.GetVectorByBits<3, 1>(coordsElevDisp[f2n[if2n - 3]], face2nodePbiExtended[iFace][if2n - 3]) +
1090 periodicInfo.GetVectorByBits<3, 1>(coordsElevDisp[f2n[if2n - 2]], face2nodePbiExtended[iFace][if2n - 2]) +
1091 periodicInfo.GetVectorByBits<3, 1>(coordsElevDisp[f2n[if2n - 1]], face2nodePbiExtended[iFace][if2n - 1]));
1092 }
1093
1094 cooSmooth = HermiteInterpolateMidPointOnQuad4WithNorm(
1095 coordsSpan[0], coordsSpan[1], coordsSpan[2], coordsSpan[3],
1096 norms(EigenAll, 0), norms(EigenAll, 1), norms(EigenAll, 2), norms(EigenAll, 3));
1097
1098 cooInc = cooSmooth - cooUnsmooth;
1099 }
1100
1101 if (isPeriodic)
1102 cooInc = periodicInfo.GetVectorBackByBits<3, 1>(cooInc, face2nodePbiExtended[iFace][if2n]);
1103 // if (cooInc.stableNorm() > 0) //! could use a threshold
1104 {
1105
1106 coordsElevDisp[iNode] = cooInc * 1; // maybe output the increment directly?
1107 nMoved++;
1108 moveds.insert(iNode);
1109 bboxMove = bboxMove.array().max(cooInc.array().abs());
1110 }
1111
1112 // std::cout << mpi.rank << " rank; iNode " << iNode << " alterwith " << coordsElevDisp[iNode].transpose() << std::endl;
1113 }
1114 }
1115 tPoint bboxMoveT;
1116 MPI::Allreduce(bboxMove.data(), bboxMoveT.data(), 3, DNDS_MPI_REAL, MPI_MAX, mpi.comm);
1117
1118 coordsElevDisp.trans.pullOnce();
1119 nMoved = moveds.size();
1120 MPI::Allreduce(&nMoved, &nTotalMoved, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
1121 if (mpi.rank == mRank)
1122 log() << fmt::format(
1123 "UnstructuredMesh === ElevatedNodesGetBoundarySmooth: Smoothing Complete, total Moved [{}], moving Vec Bnd {:.2g},{:.2g},{:.2g}",
1124 nTotalMoved, bboxMoveT(0), bboxMoveT(1), bboxMoveT(2))
1125 << std::endl;
1126
1127 // std::cout << mpi.rank << " rank Here XXXX -1" << std::endl;
1128 }
1129}
ArrayEigenUniMatrixBatch.hpp
Batch of uniform-sized Eigen matrices per row, with variable batch count.
Defines.hpp
Core type aliases, constants, and metaprogramming utilities for the DNDS framework.
EigenPCH.hpp
Pre-compiled-header style shim that includes the heavy Eigen headers under DNDSR's warning suppressio...
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
DNDS::Geom::Elem::ShapeJacobianCoordD01Nj
tJacobi ShapeJacobianCoordD01Nj(const tCoordsIn &cs, Eigen::Ref< const tD01Nj > DiNj)
Definition Elements.hpp:537
DNDS::Geom::Elem::tD01Nj
Eigen::Matrix< t_real, 4, Eigen::Dynamic > tD01Nj
Definition Elements.hpp:203
DNDS::Geom::Elem::PPhysicsCoordD01Nj
tPoint PPhysicsCoordD01Nj(const tCoordsIn &cs, Eigen::Ref< const tD01Nj > DiNj)
Definition Elements.hpp:543
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::FaceIDIsInternal
DNDS_DEVICE_CALLABLE bool FaceIDIsInternal(t_index id)
Definition BoundaryCondition.hpp:85
DNDS::Geom::tJacobi
Eigen::Matrix3d tJacobi
Definition Geometric.hpp:10
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::Geom::tSmallCoords
Eigen::Matrix< real, 3, Eigen::Dynamic > tSmallCoords
Definition Geometric.hpp:50
DNDS::MPI::Scan
MPI_int Scan(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Scan (inclusive prefix reduction).
Definition MPI.cpp:219
DNDS::MPI::Allreduce
MPI_int Allreduce(const void *sendbuf, void *recvbuf, MPI_int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
Wrapper over MPI_Allreduce.
Definition MPI.cpp:202
DNDS::MPI::Barrier
MPI_int Barrier(MPI_Comm comm)
Wrapper over MPI_Barrier.
Definition MPI.cpp:247
DNDS::DNDS_MPI_INDEX
const MPI_Datatype DNDS_MPI_INDEX
MPI datatype matching index (= MPI_INT64_T).
Definition MPI.hpp:106
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::pi
DNDS_CONSTANT const real pi
π in double precision (matches DNDS_E_PI macro).
Definition Defines.hpp:204
DNDS::UnInitMPIInt
constexpr MPI_int UnInitMPIInt
Sentinel "not initialised" MPI_int value (= -1, invalid rank).
Definition MPI.hpp:66
DNDS::largeReal
DNDS_CONSTANT const real largeReal
Loose upper bound (e.g., for non-dimensional limits).
Definition Defines.hpp:197
DNDS::index
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:112
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:522
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:110
DNDS::make_ssp
ssp< T > make_ssp(Args &&...args)
Type-safe replacement for DNDS_MAKE_SSP. Creates ssp<T> with forwarded args.
Definition Defines.hpp:260
DNDS::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:50
DNDS::veryLargeReal
DNDS_CONSTANT const real veryLargeReal
Catch-all upper bound ("practically infinity") for physical scalars.
Definition Defines.hpp:195
DNDS::MPI_int
int MPI_int
MPI counterpart type for MPI_int (= C int). Used for counts and ranks in MPI calls.
Definition MPI.hpp:54
DNDS::DNDS_MPI_REAL
const MPI_Datatype DNDS_MPI_REAL
MPI datatype matching real (= MPI_REAL8).
Definition MPI.hpp:108
DNDS::ArrayPairDeviceView_Base::RowSize
DNDS_DEVICE_CALLABLE auto RowSize() const
Uniform row width (delegates to father; father/son share it).
Definition ArrayPair.hpp:40
DNDS::ArrayPairDeviceView_Base::Size
DNDS_DEVICE_CALLABLE index Size() const
Combined father + son row count.
Definition ArrayPair.hpp:33
DNDS::ArrayPairDeviceView::father
t_arrayDeviceView father
Definition ArrayPair.hpp:108
DNDS::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:305
DNDS::ArrayPair::father
ssp< TArray > father
Owned-side array (must be resized before ghost setup).
Definition ArrayPair.hpp:170
DNDS::ArrayPair::Size
index Size() const
Combined row count (father->Size() + son->Size()).
Definition ArrayPair.hpp:296
DNDS::ArrayPair::InitPair
void InitPair(const std::string &name, Args &&...args)
Allocate both father and son arrays, forwarding all args to TArray constructor.
Definition ArrayPair.hpp:326
DNDS::ArrayPair::son
ssp< TArray > son
Ghost-side array (sized automatically by createMPITypes / BorrowAndPull).
Definition ArrayPair.hpp:172
DNDS::ArrayPair::trans
TTrans trans
Ghost-communication engine bound to father and son.
Definition ArrayPair.hpp:175
DNDS::Geom::ElemInfo::CommType
static MPI_Datatype CommType()
Definition Mesh_DeviceView.hpp:43
DNDS::Geom::NodePeriodicBits::CommType
static MPI_Datatype CommType()
Definition PeriodicInfo.hpp:41
DNDS::Geom::Periodicity::GetCoordBackByBits
DNDS_DEVICE_CALLABLE tPoint GetCoordBackByBits(const tPoint &c, const NodePeriodicBits &bits) const
Definition PeriodicInfo.hpp:360
DNDS::Geom::Periodicity::GetVectorByBits
DNDS_DEVICE_CALLABLE auto GetVectorByBits(const Eigen::Matrix< real, dim, nVec > &v, const NodePeriodicBits &bits)
Definition PeriodicInfo.hpp:346
DNDS::Geom::Periodicity::GetVectorBackByBits
DNDS_DEVICE_CALLABLE auto GetVectorBackByBits(const Eigen::Matrix< real, dim, nVec > &v, const NodePeriodicBits &bits)
Definition PeriodicInfo.hpp:375
DNDS::Geom::UnstructuredMeshDeviceView::bnd2node
AdjPairTrackedDeviceView< B, tAdjPair::t_arr > bnd2node
Definition Mesh_DeviceView.hpp:185
DNDS::Geom::UnstructuredMeshDeviceView::coords
tCoordPair::t_deviceView< B > coords
reader
Definition Mesh_DeviceView.hpp:183
DNDS::Geom::UnstructuredMeshDeviceView::GetBndElement
DNDS_DEVICE_CALLABLE Elem::Element GetBndElement(index iB)
Definition Mesh_DeviceView.hpp:329
DNDS::Geom::UnstructuredMeshDeviceView::cell2cell
AdjPairTrackedDeviceView< B, tAdjPair::t_arr > cell2cell
Definition Mesh_DeviceView.hpp:187
DNDS::Geom::UnstructuredMeshDeviceView::cellElemInfo
tElemInfoArrayPair::t_deviceView< B > cellElemInfo
Definition Mesh_DeviceView.hpp:188
DNDS::Geom::UnstructuredMeshDeviceView::cell2nodePbi
tPbiPair::t_deviceView< B > cell2nodePbi
periodic only, after reader
Definition Mesh_DeviceView.hpp:191
DNDS::Geom::UnstructuredMeshDeviceView::bnd2cell
AdjPairTrackedDeviceView< B, tAdj2Pair::t_arr > bnd2cell
Definition Mesh_DeviceView.hpp:186
DNDS::Geom::UnstructuredMeshDeviceView::GetCoordsOnCell
DNDS_DEVICE_CALLABLE void GetCoordsOnCell(index iCell, tSmallCoords &cs)
Definition Mesh_DeviceView.hpp:425
DNDS::Geom::UnstructuredMeshDeviceView::bnd2nodePbi
tPbiPair::t_deviceView< B > bnd2nodePbi
Definition Mesh_DeviceView.hpp:192
DNDS::Geom::UnstructuredMeshDeviceView::GetCellElement
DNDS_DEVICE_CALLABLE Elem::Element GetCellElement(index iC)
Definition Mesh_DeviceView.hpp:327
DNDS::Geom::UnstructuredMeshDeviceView::bndElemInfo
tElemInfoArrayPair::t_deviceView< B > bndElemInfo
Definition Mesh_DeviceView.hpp:189
DNDS::Geom::UnstructuredMeshDeviceView::cell2node
AdjPairTrackedDeviceView< B, tAdjPair::t_arr > cell2node
Definition Mesh_DeviceView.hpp:184
DNDS::Geom::UnstructuredMesh::adjC2FState
MeshAdjState adjC2FState
Definition Mesh.hpp:48
DNDS::Geom::UnstructuredMesh::t_deviceView
UnstructuredMeshDeviceView< B > t_deviceView
Definition Mesh.hpp:1036
DNDS::Geom::UnstructuredMesh::GetCellElement
Elem::Element GetCellElement(index iC)
Definition Mesh.hpp:689
DNDS::Geom::UnstructuredMesh::cell2nodePbi
tPbiPair cell2nodePbi
periodic only, after reader
Definition Mesh.hpp:68
DNDS::Geom::UnstructuredMesh::coordsElevDisp
tCoordPair coordsElevDisp
only elevation
Definition Mesh.hpp:156
DNDS::Geom::UnstructuredMesh::face2nodePbi
tPbiPair face2nodePbi
periodic only, after interpolated
Definition Mesh.hpp:107
DNDS::Geom::UnstructuredMesh::face2cell
AdjPairTracked< tAdj2Pair > face2cell
Definition Mesh.hpp:99
DNDS::Geom::UnstructuredMesh::NodeIndexGlobal2Local
index NodeIndexGlobal2Local(DNDS::index i)
Definition Mesh.hpp:320
DNDS::Geom::UnstructuredMesh::_detail_GetCoordsOnElem
void _detail_GetCoordsOnElem(const tC2n &c2n, const tC2nPbi &c2nPbi, tSmallCoords &cs)
specially for periodicity
Definition Mesh.hpp:753
DNDS::Geom::UnstructuredMesh::ElevatedNodesGetBoundarySmooth
void ElevatedNodesGetBoundarySmooth(const std::function< bool(t_index)> &FiFBndIdNeedSmooth)
Definition Mesh_Elevation.cpp:699
DNDS::Geom::UnstructuredMesh::coords
tCoordPair coords
reader
Definition Mesh.hpp:57
DNDS::Geom::UnstructuredMesh::cell2face
AdjPairTracked< tAdjPair > cell2face
interpolated
Definition Mesh.hpp:97
DNDS::Geom::UnstructuredMesh::BuildBisectO1FormO2
void BuildBisectO1FormO2(UnstructuredMesh &meshO2)
Definition Mesh_Elevation.cpp:476
DNDS::Geom::UnstructuredMesh::getCell2NodeIndexPbiRow
auto getCell2NodeIndexPbiRow(index iCell)
Definition Mesh.hpp:184
DNDS::Geom::UnstructuredMesh::cell2cell
AdjPairTracked< tAdjPair > cell2cell
Definition Mesh.hpp:61
DNDS::Geom::UnstructuredMesh::_detail_GetCoords
void _detail_GetCoords(const tC2n &c2n, tSmallCoords &cs)
directly load coords; gets faulty if isPeriodic!
Definition Mesh.hpp:715
DNDS::Geom::UnstructuredMesh::face2node
AdjPairTracked< tAdjPair > face2node
Definition Mesh.hpp:98
DNDS::Geom::UnstructuredMesh::adjPrimaryState
MeshAdjState adjPrimaryState
Definition Mesh.hpp:44
DNDS::Geom::UnstructuredMesh::bnd2face
AdjPairTracked< tAdj1Pair > bnd2face
Definition Mesh.hpp:101
DNDS::Geom::UnstructuredMesh::faceElemInfo
tElemInfoArrayPair faceElemInfo
Definition Mesh.hpp:100
DNDS::Geom::UnstructuredMesh::elevationInfo
struct DNDS::Geom::UnstructuredMesh::ElevationInfo elevationInfo
DNDS::Geom::UnstructuredMesh::adjFacialState
MeshAdjState adjFacialState
Definition Mesh.hpp:46
DNDS::Geom::UnstructuredMesh::BuildO2FromO1Elevation
void BuildO2FromO1Elevation(UnstructuredMesh &meshO1)
Definition Mesh_Elevation.cpp:22
DNDS::Geom::UnstructuredMesh::cell2node
AdjPairTracked< tAdjPair > cell2node
Definition Mesh.hpp:58
DNDS::Geom::UnstructuredMesh::cellElemInfo
tElemInfoArrayPair cellElemInfo
Definition Mesh.hpp:62
DNDS::Geom::UnstructuredMesh::bnd2node
AdjPairTracked< tAdjPair > bnd2node
Definition Mesh.hpp:59
DNDS::Geom::UnstructuredMesh::bndElemInfo
tElemInfoArrayPair bndElemInfo
Definition Mesh.hpp:63
DNDS::Geom::UnstructuredMesh::elevState
MeshElevationState elevState
Definition Mesh.hpp:55
DNDS::Geom::UnstructuredMesh::bnd2cell
AdjPairTracked< tAdj2Pair > bnd2cell
Definition Mesh.hpp:60
DNDS::MPIInfo::size
int size
Number of ranks in comm (-1 until initialised).
Definition MPI.hpp:237
DNDS::MPIInfo::rank
int rank
This rank's 0-based index within comm (-1 until initialised).
Definition MPI.hpp:235
DNDS::MPIInfo::comm
MPI_Comm comm
The underlying MPI communicator handle.
Definition MPI.hpp:233
DNDS::ObjName
Tag type for naming objects created via make_ssp.
Definition Defines.hpp:254
setConstant
dof1 setConstant(0.0)
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
result
auto result
Definition test_Limiters.cpp:110
tPoint
Eigen::Vector3d tPoint
Definition test_MeshConnectivity_Interpolate.cpp:1771
np
DNDS::index np
Definition test_MeshConnectivity_Interpolate.cpp:828
i
res i
Definition test_MeshConnectivity_Interpolate.cpp:205