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