DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Mesh_Elevation_SmoothV2.cpp
Go to the documentation of this file.
2#include "Geom/Quadrature.hpp"
3
4#include <fmt/core.h>
5#include <Solver/Linear.hpp>
6
7#include "DNDS/EigenPCH.hpp"
8
9namespace DNDS::Geom
10{
11 // =================================================================
12 // Block sparse matrix element for the elasticity system.
13 // =================================================================
14 namespace
15 {
16 struct MatElem
17 {
19 // NOLINTNEXTLINE(readability-redundant-member-init): Eigen Matrix default-init is uninitialized
20 Eigen::Matrix<real, 3, 3> m{};
21 };
22
23 // Build node-to-node adjacency from cell connectivity.
24 std::vector<std::unordered_set<index>> BuildNode2NodeFromCells(
25 const tAdjPair &cell2node, index nLocalNodes)
26 {
27 std::vector<std::unordered_set<index>> node2nodeV(nLocalNodes);
28 for (index iCell = 0; iCell < cell2node.Size(); iCell++)
29 for (auto iN : cell2node[iCell])
30 for (auto iNOther : cell2node[iCell])
31 {
32 if (iN == iNOther)
33 continue;
34 if (iN < nLocalNodes)
35 node2nodeV[iN].insert(iNOther);
36 }
37 return node2nodeV;
38 }
39
40 // Allocate block sparse matrix from node-to-node adjacency.
41 std::vector<std::vector<MatElem>> AllocateBlockSparseMatrix(
42 const std::vector<std::unordered_set<index>> &node2nodeV,
43 index nLocalNodes)
44 {
45 std::vector<std::vector<MatElem>> A(nLocalNodes);
46 for (index iN = 0; iN < nLocalNodes; iN++)
47 {
48 A[iN].resize(node2nodeV[iN].size() + 1);
49 A[iN][0].j = iN;
50 int curRowEnd = 1;
51 for (auto iNOther : node2nodeV[iN])
52 A[iN][curRowEnd++].j = iNOther;
53 for (auto &ME : A[iN])
54 ME.m.setConstant(0);
55 }
56 return A;
57 }
58
59 // Block-sparse matrix-vector product with MPI ghost pull.
60 void BlockMatVec(
61 const std::vector<std::vector<MatElem>> &A,
63 {
64 for (index iN = 0; iN < x.father->Size(); iN++)
65 {
66 Ax[iN].setZero();
67 for (const auto &ME : A[iN])
68 Ax[iN] += ME.m * x[ME.j];
69 }
70 Ax.trans.startPersistentPull();
71 Ax.trans.waitPersistentPull();
72 }
73
74 // Block-Jacobi preconditioner with MPI ghost pull.
75 void BlockJacobiPrecondition(
76 const std::vector<std::vector<MatElem>> &A,
77 tCoordPair &x, tCoordPair &ADIx)
78 {
79 for (index iN = 0; iN < x.father->Size(); iN++)
80 ADIx[iN] = A[iN][0].m.fullPivLu().solve(x[iN]);
81 ADIx.trans.startPersistentPull();
82 ADIx.trans.waitPersistentPull();
83 }
84 } // anonymous namespace
85
86 // =================================================================
87 // ElevatedNodesSolveInternalSmoothV2
88 // =================================================================
90 {
92 if (!setupOpt)
93 return;
94 auto &setup = *setupOpt;
95 auto &nodesBoundInterpolated = setup.nodesBoundInterpolated;
96 auto &boundInterpCoo = setup.boundInterpCoo;
97
98 // ----- Build node-to-node adjacency -----
100
101 if (mpi.rank == mRank)
102 log() << "RBF set: " << boundInterpCoo.son->Size() << std::endl;
103
104 // ----- KD-tree for nearest boundary point lookup -----
105 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
106 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTreeCoordPair>,
108 auto coordsI = PointCloudKDTreeCoordPair(boundInterpCoo.son);
110
111 DNDS_MPI_InsertCheck(mpi, "Got node2node");
112
113 // ----- Initialize DOF vectors -----
115 auto initDOF = [&](CoordPairDOF &v)
116 {
117 v.InitPair("V2::dof", mpi);
118 v.father->Resize(coords.father->Size());
119 v.TransAttach();
120 v.trans.BorrowGGIndexing(coords.trans);
121 v.trans.createMPITypes();
122 v.trans.initPersistentPull();
123 };
128 initDOF(bO2);
130
131 // ----- Allocate block sparse matrix -----
133
134 // Initialize DOF values from boundary displacement data.
135 for (index iN = 0; iN < coords.father->Size(); iN++)
136 {
137 bO2[iN].setZero();
138 dispO2[iN].setZero();
139 if (coordsElevDisp[iN](0) != largeReal)
141 if (coordsElevDisp[iN](2) == 2 * largeReal ||
142 coordsElevDisp[iN](2) == 3 * largeReal)
143 dispO2[iN].setZero();
144 }
145 dispO2.trans.startPersistentPull();
146 dispO2.trans.waitPersistentPull();
147
148 // =============================================================
149 // Assemble elasticity stiffness matrix and RHS.
150 // Uses bisect-FEM: each O2 cell is bisected into O1 sub-elements,
151 // linear elasticity stiffness is assembled on each sub-element.
152 // =============================================================
153 auto AssembleRHSMat = [&](tCoordPair &uCur)
154 {
155 index nFix{0}, nFixB{0};
156
157 for (index iCell = 0; iCell < cell2node.Size(); iCell++)
158 {
159 auto c2n = cell2node[iCell];
160 rowsize nnLoc = c2n.size();
164 tPoint cellCent = coordsC.rowwise().mean();
165
166 // Find distance to nearest boundary point.
167 real wDist{0};
168 {
169 std::vector<real> sqrDists(1);
170 std::vector<index> idxFound(1);
171 auto nFound = bndInterpTree.knnSearch(
172 cellCent.data(), 1, idxFound.data(), sqrDists.data());
173 DNDS_assert(nFound == 1);
174 wDist = std::sqrt(sqrDists[0]);
175 }
176
178 ALoc.setZero(static_cast<index>(nnLoc) * 3, static_cast<index>(nnLoc) * 3 + 1);
179 auto localNodeIdx = Eigen::ArrayXi::LinSpaced(c2n.size(), 0, c2n.size() - 1);
180
181 // Bisect the O2 element into O1 sub-elements.
182 auto eCell = GetCellElement(iCell);
183 int nBi = eCell.GetO2NumBisect();
184 int iBiVariant = GetO2ElemBisectVariant(eCell, coordsC);
185 for (int iBi = 0; iBi < nBi; iBi++)
186 {
187 auto eCellSub = eCell.ObtainO2BisectElem(iBi);
188 Eigen::ArrayXi c2nSubLocal;
189 c2nSubLocal.resize(eCellSub.GetNumNodes());
190 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, localNodeIdx, c2nSubLocal);
192 coordsCSub.resize(3, eCellSub.GetNumNodes());
193 coordsCuSub.resize(3, eCellSub.GetNumNodes());
194 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsC, coordsCSub);
195 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsCu, coordsCuSub);
197 auto nnLocSub = rowsize(c2nSubLocal.size());
199 ALocSub.setZero(static_cast<index>(nnLocSub) * 3, static_cast<index>(nnLocSub) * 3 + 1);
200 mLoc.resize(6, static_cast<index>(nnLocSub) * 3);
201 Eigen::ArrayXi c2nSubLocal3;
202 c2nSubLocal3.resize(c2nSubLocal.size() * 3);
203 c2nSubLocal3(Eigen::seq(c2nSubLocal.size() * 0, c2nSubLocal.size() * 1 - 1)) = c2nSubLocal + nnLoc * 0;
204 c2nSubLocal3(Eigen::seq(c2nSubLocal.size() * 1, c2nSubLocal.size() * 2 - 1)) = c2nSubLocal + nnLoc * 1;
205 c2nSubLocal3(Eigen::seq(c2nSubLocal.size() * 2, c2nSubLocal.size() * 3 - 1)) = c2nSubLocal + nnLoc * 2;
206
207 // Integrate elasticity stiffness on the sub-element.
208 qCellSub.Integration(
209 ALocSub,
210 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
211 {
212 auto getJ = [&](auto coos)
213 {
215 tJacobi JInv = tJacobi::Identity();
216 real JDet = UnInitReal;
217 if (dim == 2)
218 JDet = J(EigenAll, 0).cross(J(EigenAll, 1)).stableNorm();
219 else
220 JDet = J.determinant();
221 if (dim == 2)
222 JInv({0, 1}, {0, 1}) = J({0, 1}, {0, 1}).inverse().eval();
223 else
224 JInv = J.inverse().eval();
225 return std::make_tuple(J, JInv, JDet);
226 };
227 auto [J, JInv, JDet] = getJ(coordsCSub);
228 auto [Jn, JInvn, JDetn] = getJ(coordsCSub + coordsCuSub);
229 DNDS_assert(JDet > 0);
230
231 // Strain-displacement matrix (Voigt notation).
232 tSmallCoords m3 = JInv.transpose() * DiNj({1, 2, 3}, EigenAll);
233 mLoc.setZero();
234 mLoc(0, Eigen::seq(nnLocSub * 0, nnLocSub * 0 + nnLocSub - 1)) = m3(0, EigenAll);
235 mLoc(1, Eigen::seq(nnLocSub * 1, nnLocSub * 1 + nnLocSub - 1)) = m3(1, EigenAll);
236 mLoc(2, Eigen::seq(nnLocSub * 2, nnLocSub * 2 + nnLocSub - 1)) = m3(2, EigenAll);
237 mLoc(3, Eigen::seq(nnLocSub * 1, nnLocSub * 1 + nnLocSub - 1)) = 0.5 * m3(2, EigenAll);
238 mLoc(3, Eigen::seq(nnLocSub * 2, nnLocSub * 2 + nnLocSub - 1)) = 0.5 * m3(1, EigenAll);
239 mLoc(4, Eigen::seq(nnLocSub * 2, nnLocSub * 2 + nnLocSub - 1)) = 0.5 * m3(0, EigenAll);
240 mLoc(4, Eigen::seq(nnLocSub * 0, nnLocSub * 0 + nnLocSub - 1)) = 0.5 * m3(2, EigenAll);
241 mLoc(5, Eigen::seq(nnLocSub * 0, nnLocSub * 0 + nnLocSub - 1)) = 0.5 * m3(1, EigenAll);
242 mLoc(5, Eigen::seq(nnLocSub * 1, nnLocSub * 1 + nnLocSub - 1)) = 0.5 * m3(0, EigenAll);
243
244 // Constitutive matrix (isotropic linear elasticity).
245 Eigen::Matrix<real, 6, 6> DStruct;
246 DStruct.setIdentity();
247 real lam = 1;
248 real muu = 100;
249 DStruct *= muu;
250 DStruct({0, 1, 2}, {0, 1, 2}) += muu * tJacobi::Identity();
251 DStruct({0, 1, 2}, {0, 1, 2}).array() += lam;
252 DStruct *= std::pow(JDetn / JDet, -0.0) + 1;
253
254 vInc.resizeLike(ALocSub);
255 vInc(EigenAll, Eigen::seq(0, 3 * nnLocSub - 1)) = mLoc.transpose() * DStruct * mLoc;
256 vInc(EigenAll, EigenLast).setZero();
257 vInc *= JDet;
258 });
259 ALoc(c2nSubLocal3, c2nSubLocal3) += ALocSub(EigenAll, Eigen::seq(0, EigenLast - 1));
260 ALoc(c2nSubLocal3, EigenLast) += ALocSub(EigenAll, EigenLast);
261 }
262
263 // Zero out z-component rows/cols for 2D problems.
264 if (dim != 3)
265 {
266 ALoc(Eigen::seq(2 * nnLoc, 3 * nnLoc - 1), Eigen::seq(2 * nnLoc, 3 * nnLoc - 1)).setZero();
267 ALoc(EigenAll, Eigen::seq(2 * nnLoc, 3 * nnLoc - 1)).setZero();
268 ALoc(Eigen::seq(2 * nnLoc, 3 * nnLoc - 1), EigenAll).setZero();
269 }
270
271 // Scatter local matrix into global block-sparse matrix.
272 for (rowsize ic2n = 0; ic2n < nnLoc; ic2n++)
273 {
274 index iN = c2n[ic2n];
275 if (iN >= coords.father->Size())
276 continue;
277 for (rowsize jc2n = 0; jc2n < nnLoc; jc2n++)
278 {
279 index jN = c2n[jc2n];
280 int nMatrixFound = 0;
281 for (auto &ME : A[iN])
282 if (ME.j == jN)
283 {
284 nMatrixFound++;
285 ME.m(0, 0) += ALoc(0 * nnLoc + ic2n, 0 * nnLoc + jc2n);
286 ME.m(0, 1) += ALoc(0 * nnLoc + ic2n, 1 * nnLoc + jc2n);
287 ME.m(0, 2) += ALoc(0 * nnLoc + ic2n, 2 * nnLoc + jc2n);
288 ME.m(1, 0) += ALoc(1 * nnLoc + ic2n, 0 * nnLoc + jc2n);
289 ME.m(1, 1) += ALoc(1 * nnLoc + ic2n, 1 * nnLoc + jc2n);
290 ME.m(1, 2) += ALoc(1 * nnLoc + ic2n, 2 * nnLoc + jc2n);
291 ME.m(2, 0) += ALoc(2 * nnLoc + ic2n, 0 * nnLoc + jc2n);
292 ME.m(2, 1) += ALoc(2 * nnLoc + ic2n, 1 * nnLoc + jc2n);
293 ME.m(2, 2) += ALoc(2 * nnLoc + ic2n, 2 * nnLoc + jc2n);
294 if (isPeriodic)
295 {
296 auto ipbi = cell2nodePbi(iCell, ic2n);
297 auto jpbi = cell2nodePbi(iCell, jc2n);
298 tJacobi iTrans = periodicInfo.GetVectorByBits<3, 3>(tJacobi::Identity(), ipbi);
299 tJacobi jTrans = periodicInfo.GetVectorByBits<3, 3>(tJacobi::Identity(), jpbi);
300 ME.m = jTrans.transpose() * ME.m * iTrans;
301 }
302 }
304 }
305
306 // RHS contribution.
307 bO2[iN](0) += ALoc(0 * nnLoc + ic2n, EigenLast);
308 bO2[iN](1) += ALoc(1 * nnLoc + ic2n, EigenLast);
309 bO2[iN](2) += ALoc(2 * nnLoc + ic2n, EigenLast);
310
311 // Apply Dirichlet BC for boundary-displaced nodes.
312 if (coordsElevDisp[iN](0) != largeReal)
313 {
314 nFix++;
315 real nDiag = 0;
316 for (auto &ME : A[iN])
317 if (ME.j == iN)
318 nDiag = ME.m.array().abs().maxCoeff();
320 for (auto &ME : A[iN])
321 if (ME.j == iN)
322 ME.m = tJacobi::Identity() * nDiag;
323 else
324 ME.m.setZero();
325 }
326 // Apply Dirichlet BC for O1 nodes and other boundary markers.
327 if (coordsElevDisp[iN](2) == 2 * largeReal ||
328 coordsElevDisp[iN](2) == 3 * largeReal)
329 {
330 nFixB++;
331 real nDiag = 0;
332 for (auto &ME : A[iN])
333 if (ME.j == iN)
334 nDiag = ME.m.array().abs().maxCoeff();
335 bO2[iN].setZero();
336 for (auto &ME : A[iN])
337 if (ME.j == iN)
338 ME.m = tJacobi::Identity() * nDiag;
339 else
340 ME.m.setZero();
341 }
342 // Constrain z-direction in 2D.
343 if (dim != 3)
344 {
345 real nDiag = 0;
346 for (auto &ME : A[iN])
347 if (ME.j == iN)
348 nDiag = ME.m.array().abs().maxCoeff();
349 for (auto &ME : A[iN])
350 if (ME.j == iN)
351 ME.m(2, 2) += nDiag;
352 }
353 }
354 }
355
356 bO2.trans.startPersistentPull();
357 bO2.trans.waitPersistentPull();
360 if (mpi.rank == mRank)
361 {
362 log() << fmt::format(
363 "UnstructuredMesh === ElevatedNodesSolveInternalSmooth(): "
364 "Matrix Assembled, nFix {}, nFixB {} ",
365 nFix, nFixB)
366 << std::endl;
367 }
368 };
369
370 // ----- Assemble the linear system -----
371 dispO2.setConstant(0.0);
373
374 // ----- Solve via preconditioned GMRES -----
376 gmres.solve(
378 { BlockMatVec(A, x, Ax); },
381 [&](CoordPairDOF &a, CoordPairDOF &b) -> real
382 { return a.dot(b); },
383 bO2,
384 dispO2,
385 elevationInfo.nIter,
386 [&](int iRestart, real res, real resB)
387 {
388 if (mpi.rank == mRank)
389 log() << fmt::format("iRestart [{}] res [{:3e}] -> [{:3e}]",
390 iRestart, resB, res)
391 << std::endl;
392 return res < resB * 1e-10;
393 });
394
395 // ----- Apply solved displacement to coordinates -----
396 for (index iN = 0; iN < coords.father->Size(); iN++)
397 {
398 if (dim == 2)
399 dispO2[iN](2) = 0;
400 coords[iN] += dispO2[iN];
401 }
402 coords.trans.pullOnce();
403 }
404
405} // namespace DNDS::Geom
Pre-compiled-header style shim that includes the heavy Eigen headers under DNDSR's warning suppressio...
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:112
#define DNDS_MPI_InsertCheck(mpi, info)
Debug helper: barrier + print "CHECK <info> RANK r @ fn @ file:line".
Definition MPI.hpp:343
Eigen::Matrix< real, 3, 3 > m
bool solve(TFA &&FA, TFML &&FML, TFDot &&fDot, TDATA &b, TDATA &x, uint32_t nRestart, TFstop &&FStop)
Definition Linear.hpp:47
tJacobi ShapeJacobianCoordD01Nj(const tCoordsIn &cs, Eigen::Ref< const tD01Nj > DiNj)
Definition Elements.hpp:537
Eigen::Matrix< t_real, 4, Eigen::Dynamic > tD01Nj
Definition Elements.hpp:203
Eigen::Matrix3d tJacobi
Definition Geometric.hpp:10
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::ArrayPair< DNDS::ArrayEigenVector< 3 > > tCoordPair
Eigen::Matrix< real, 3, Eigen::Dynamic > tSmallCoords
Definition Geometric.hpp:50
std::optional< SmoothSolverSetup > PrepareSmoothSolverSetup(UnstructuredMesh &mesh)
Common preamble for all smooth solver variants.
DNDS::ArrayAdjacencyPair< DNDS::NonUniformSize > tAdjPair
void AllreduceOneIndex(index &v, MPI_Op op, const MPIInfo &mpi)
Single-scalar Allreduce helper for indices (in-place, count = 1).
Definition MPI.hpp:727
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:181
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:114
DNDS_CONSTANT const real UnInitReal
Sentinel "not initialised" real value (NaN). Cheap to detect with std::isnan or IsUnInitReal; survive...
Definition Defines.hpp:179
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:210
DNDS_CONSTANT const real largeReal
Loose upper bound (e.g., for non-dimensional limits).
Definition Defines.hpp:197
int64_t index
Global row / DOF index type (signed 64-bit; handles multi-billion-cell meshes).
Definition Defines.hpp:112
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:110
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:50
ssp< TArray > father
Owned-side array (must be resized before ghost setup).
index Size() const
Combined row count (father->Size() + son->Size()).
TTrans trans
Ghost-communication engine bound to father and son.
DNDS_DEVICE_CALLABLE auto GetVectorByBits(const Eigen::Matrix< real, dim, nVec > &v, const NodePeriodicBits &bits)
UnstructuredMeshDeviceView< B > t_deviceView
Definition Mesh.hpp:1036
Elem::Element GetCellElement(index iC)
Definition Mesh.hpp:689
tPbiPair cell2nodePbi
periodic only, after reader
Definition Mesh.hpp:68
tCoordPair coordsElevDisp
only elevation
Definition Mesh.hpp:156
void GetCoordsOnCell(index iCell, tSmallCoords &cs)
Definition Mesh.hpp:787
tCoordPair coords
reader
Definition Mesh.hpp:57
struct DNDS::Geom::UnstructuredMesh::ElevationInfo elevationInfo
AdjPairTracked< tAdjPair > cell2node
Definition Mesh.hpp:58
int rank
This rank's 0-based index within comm (-1 until initialised).
Definition MPI.hpp:235
dof1 setConstant(0.0)
Eigen::Matrix< real, 5, 1 > v
tVec Ax(NCells)
tVec x(NCells)
tVec b(NCells)
if(g_mpi.rank==0) std auto[rhsDensity, rhsEnergy, incL2]
for(int i=0;i< 10;i++) theta(i