DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Mesh_Elevation_SmoothSolver.cpp
Go to the documentation of this file.
1#include "Mesh.hpp"
2#include "Quadrature.hpp"
3#include "DNDS/ArrayDerived/ArrayEigenUniMatrixBatch.hpp"
4#include "RadialBasisFunction.hpp"
5
6#include <omp.h>
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 struct CoordPairDOF : public tCoordPair
23 {
24 real dot(CoordPairDOF &R)
25 {
26 real ret = 0;
27 for (index i = 0; i < this->father->Size(); i++)
28 ret += (*this)[i].dot(R[i]);
29 real retSum;
30 MPI::Allreduce(&ret, &retSum, 1, DNDS_MPI_REAL, MPI_SUM, this->father->getMPI().comm);
31 return retSum;
32 }
33
34 real norm2()
35 {
36 return std::sqrt(this->dot(*this));
37 }
38
39 void addTo(CoordPairDOF &R, real alpha)
40 {
41 for (index i = 0; i < this->Size(); i++)
42 (*this)[i] += R[i] * alpha;
43 }
44
45 void setConstant(real v)
46 {
47 for (index i = 0; i < this->Size(); i++)
48 (*this)[i].setConstant(v);
49 }
50
51 void operator=(CoordPairDOF &R)
52 {
53 for (index i = 0; i < this->Size(); i++)
54 (*this)[i] = R[i];
55 }
56
57 void operator*=(real r)
58 {
59 for (index i = 0; i < this->Size(); i++)
60 (*this)[i] *= r;
61 }
62 };
63
64 struct PointCloudKDTreeCoordPair
65 {
66 tCoord ref;
67 using coord_t = real; //!< The type of each coordinate
72
73 // Must return the number of data points
74 [[nodiscard]] inline size_t
75 kdtree_get_point_count() const
76 {
77 DNDS_assert(ref);
78 return ref->Size();
79 }
80
81 [[nodiscard]] inline real kdtree_get_pt(const size_t idx, const size_t dim) const
82 {
83 DNDS_assert(ref);
84 return ref->operator[](idx)(dim);
85 }
86
87 template <class BBOX>
88 bool kdtree_get_bbox(BBOX & /* bb */) const
89 {
90 return false;
91 }
92 };
93
94 void UnstructuredMesh::ElevatedNodesSolveInternalSmooth()
95 {
96 DNDS_assert(elevState == Elevation_O1O2);
97 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
98 DNDS_assert(adjFacialState == Adj_PointToLocal);
99 DNDS_assert(adjC2FState == Adj_PointToLocal);
100 DNDS_assert(face2node.father);
101 DNDS_assert(nTotalMoved >= 0);
102 if (!nTotalMoved)
103 {
104 if (mpi.rank == mRank)
105 log() << "UnstructuredMesh === ElevatedNodesSolveInternalSmooth() early exit for no nodes were moved";
106 return;
107 }
108
109 DNDS_MPI_InsertCheck(mpi, "Got node2node");
110
111 coordsElevDisp.trans.initPersistentPull(); // only holds local nodes
112
113 std::unordered_set<index> nodesBoundInterpolated;
114 for (index iN = 0; iN < coords.father->Size(); iN++)
115 {
116 if (coordsElevDisp[iN](0) != largeReal || coordsElevDisp[iN](2) == 2 * largeReal)
117 {
118 nodesBoundInterpolated.insert(iN);
119 }
120 }
121
122 tCoordPair boundInterpCoo;
123 tCoordPair boundInterpVal;
124 boundInterpCoo.InitPair("ElevatedNodesSolveInternalSmooth::boundInterpCoo", mpi);
125 boundInterpVal.InitPair("ElevatedNodesSolveInternalSmooth::boundInterpVal", mpi);
126
127 boundInterpCoo.father->Resize(nodesBoundInterpolated.size());
128 boundInterpVal.father->Resize(nodesBoundInterpolated.size());
129
130 index top{0};
131 for (auto iN : nodesBoundInterpolated)
132 {
133 boundInterpCoo[top] = coords[iN];
134 boundInterpVal[top] = (coordsElevDisp[iN](0) != largeReal) ? tPoint{coordsElevDisp[iN]} : tPoint::Zero();
135 top++;
136 }
137 std::vector<index> boundInterpPullIdx;
138
139 index boundInterpGlobSize = boundInterpCoo.father->globalSize();
140 boundInterpPullIdx.resize(boundInterpGlobSize);
141 for (index i = 0; i < boundInterpPullIdx.size(); i++)
142 boundInterpPullIdx[i] = i;
143 boundInterpCoo.father->createGlobalMapping();
144 boundInterpCoo.TransAttach();
145 boundInterpCoo.trans.createGhostMapping(boundInterpPullIdx);
146 boundInterpCoo.trans.createMPITypes();
147 boundInterpCoo.trans.pullOnce();
148
149 boundInterpVal.TransAttach();
150 boundInterpVal.trans.BorrowGGIndexing(boundInterpCoo.trans);
151 boundInterpVal.trans.createMPITypes();
152 boundInterpVal.trans.pullOnce();
153
154 if (mpi.rank == mRank)
155 log() << "RBF set: " << boundInterpCoo.son->Size() << std::endl;
156
157 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
158 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTreeCoordPair>,
159 PointCloudKDTreeCoordPair,
160 3,
161 index>;
162 auto coordsI = PointCloudKDTreeCoordPair(boundInterpCoo.son);
163 // for (index iF = 0; iF < boundInterpCoo.son->Size(); iF++)
164 // {
165 // std::cout << boundInterpVal.sfaon->operator[](iF).transpose() << std::endl;
166 // }
167 kdtree_t bndInterpTree(3, coordsI);
168
169 for (index iN = 0; iN < coords.father->Size(); iN++)
170 {
171 if (nodesBoundInterpolated.count(iN))
172 continue;
173
174 index nFind = elevationInfo.nSearch;
175 tPoint cooC = coords[iN];
176
177 std::vector<index> idxFound;
178 Eigen::Vector<real, Eigen::Dynamic> outDistancesSqr;
179 idxFound.resize(nFind);
180 outDistancesSqr.resize(nFind);
181
182 index nFound = bndInterpTree.knnSearch(cooC.data(), nFind, idxFound.data(), outDistancesSqr.data());
183 DNDS_assert(nFound >= 1);
184 idxFound.resize(nFound);
185 outDistancesSqr.resize(nFound);
186 tSmallCoords coordBnd;
187 tSmallCoords dispBnd;
188 coordBnd.resize(3, nFound);
189 dispBnd.resize(3, nFound);
190 for (index iF = 0; iF < nFound; iF++)
191 {
192 coordBnd(EigenAll, iF) = (*boundInterpCoo.son)[idxFound[iF]];
193 dispBnd(EigenAll, iF) = (*boundInterpVal.son)[idxFound[iF]];
194 }
195 real RMin = std::sqrt(outDistancesSqr.minCoeff());
196 tPoint dispC;
197 dispC.setZero();
198
199 if (RMin < sqr(elevationInfo.RBFRadius * (KernelIsCompact(elevationInfo.kernel) ? 1 : 5)))
200 {
201 tPoint coordBndC = coordBnd.rowwise().mean();
202 tSmallCoords coordBndRel = coordBnd.colwise() - coordBndC;
203 MatrixXR
204 coefs = RBF::RBFInterpolateSolveCoefsNoPoly(coordBndRel, dispBnd.transpose(), elevationInfo.RBFRadius, elevationInfo.kernel);
205 tSmallCoords qs;
206 qs.resize(3, 1);
207 qs = cooC - coordBndC;
208 tPoint dCur =
209 (RBF::RBFCPC2(coordBndRel, qs, elevationInfo.RBFRadius, elevationInfo.kernel).transpose() * coefs.topRows(coordBndRel.cols()) +
210 coefs(EigenLast - 3, EigenAll) +
211 qs.transpose() * coefs.bottomRows(3))
212 .transpose();
213 dispC = dCur;
214 }
215
216 coordsElevDisp[iN] = dispC;
217 }
218
219 coordsElevDisp.trans.startPersistentPull();
220 coordsElevDisp.trans.waitPersistentPull();
221
222 for (index iN = 0; iN < coords.father->Size(); iN++)
223 {
224 // if(dim == 2)
225 // dispO2[iN](2) = 0;
226 // if(dispO2[iN].norm() != 0)
227 // std::cout << dispO2[iN].transpose() << std::endl;
228 if (coordsElevDisp[iN](0) != largeReal)
229 coords[iN] += coordsElevDisp[iN];
230 }
231 coords.trans.pullOnce();
232 }
233
234 void UnstructuredMesh::ElevatedNodesSolveInternalSmoothV1Old()
235 {
236 DNDS_assert(elevState == Elevation_O1O2);
237 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
238 DNDS_assert(adjFacialState == Adj_PointToLocal);
239 DNDS_assert(adjC2FState == Adj_PointToLocal);
240 DNDS_assert(face2node.father);
241 DNDS_assert(nTotalMoved >= 0);
242 if (!nTotalMoved)
243 {
244 if (mpi.rank == mRank)
245 log() << "UnstructuredMesh === ElevatedNodesSolveInternalSmooth() early exit for no nodes were moved";
246 return;
247 }
248
249 DNDS_MPI_InsertCheck(mpi, "Got node2node");
250
251 coordsElevDisp.trans.initPersistentPull(); // only holds local nodes
252
253 std::unordered_set<index> nodesBoundInterpolated;
254 for (index iN = 0; iN < coords.father->Size(); iN++)
255 {
256 if (coordsElevDisp[iN](0) != largeReal || coordsElevDisp[iN](2) == 2 * largeReal)
257 {
258 nodesBoundInterpolated.insert(iN);
259 }
260 }
261
262 tCoordPair boundInterpCoo;
263 tCoordPair boundInterpVal;
264 boundInterpCoo.InitPair("ElevatedNodesSolveInternalSmoothV1Old::boundInterpCoo", mpi);
265 boundInterpVal.InitPair("ElevatedNodesSolveInternalSmoothV1Old::boundInterpVal", mpi);
266
267 boundInterpCoo.father->Resize(nodesBoundInterpolated.size());
268 boundInterpVal.father->Resize(nodesBoundInterpolated.size());
269
270 index top{0};
271 for (auto iN : nodesBoundInterpolated)
272 {
273 boundInterpCoo[top] = coords[iN];
274 boundInterpVal[top] = (coordsElevDisp[iN](0) != largeReal) ? tPoint{coordsElevDisp[iN]} : tPoint::Zero();
275 top++;
276 }
277 std::vector<index> boundInterpPullIdx;
278
279 index boundInterpGlobSize = boundInterpCoo.father->globalSize();
280 boundInterpPullIdx.resize(boundInterpGlobSize);
281 for (index i = 0; i < boundInterpPullIdx.size(); i++)
282 boundInterpPullIdx[i] = i;
283 boundInterpCoo.father->createGlobalMapping();
284 boundInterpCoo.TransAttach();
285 boundInterpCoo.trans.createGhostMapping(boundInterpPullIdx);
286 boundInterpCoo.trans.createMPITypes();
287 boundInterpCoo.trans.pullOnce();
288
289 boundInterpVal.TransAttach();
290 boundInterpVal.trans.BorrowGGIndexing(boundInterpCoo.trans);
291 boundInterpVal.trans.createMPITypes();
292 boundInterpVal.trans.pullOnce();
293
294 if (mpi.rank == mRank)
295 log() << "RBF set: " << boundInterpCoo.son->Size() << std::endl;
296
297 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
298 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTreeCoordPair>,
299 PointCloudKDTreeCoordPair,
300 3,
301 index>;
302 using kdtree_tcoo = nanoflann::KDTreeSingleIndexAdaptor<
303 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTree>,
304 PointCloudKDTree,
305 3,
306 index>;
307 auto coordsI = PointCloudKDTreeCoordPair(boundInterpCoo.son);
308 kdtree_t bndInterpTree(3, coordsI);
309
310 tCoordPair boundInterpCoef;
311 using tScalarPair = DNDS::ArrayPair<DNDS::ParArray<real, 1>>;
312 tScalarPair boundInterpR;
313
314 boundInterpCoef.InitPair("ElevatedNodesSolveInternalSmoothV1Old::boundInterpCoef", mpi);
315 boundInterpR.InitPair("ElevatedNodesSolveInternalSmoothV1Old::boundInterpR", mpi);
316 boundInterpCoef.father->Resize(mpi.rank == mRank ? boundInterpGlobSize : 0);
317 boundInterpR.father->Resize(mpi.rank == mRank ? boundInterpGlobSize : 0);
318
319 if (mpi.rank == mRank) // that dirty work
320 {
321
322 Eigen::SparseMatrix<real> M(boundInterpGlobSize, boundInterpGlobSize);
323 M.reserve(Eigen::Vector<index, Eigen::Dynamic>::Constant(boundInterpGlobSize, elevationInfo.nSearch));
324 MatrixXR f;
325 f.resize(boundInterpGlobSize, 3);
326 for (index iN = 0; iN < boundInterpGlobSize; iN++)
327 {
328 std::cout << iN << std::endl;
329 tPoint cooC = (*boundInterpCoo.son)[iN];
330
331 index nFind = elevationInfo.nSearch;
332 std::vector<index> idxFound;
333 Eigen::Vector<real, Eigen::Dynamic> outDistancesSqr;
334 idxFound.resize(nFind);
335 outDistancesSqr.resize(nFind);
336 index nFound = bndInterpTree.knnSearch(cooC.data(), nFind, idxFound.data(), outDistancesSqr.data());
337 DNDS_assert(nFound >= 1);
338 idxFound.resize(nFound);
339 outDistancesSqr.resize(nFound);
340
341 real RMin = outDistancesSqr.minCoeff();
342 real RMax = outDistancesSqr.maxCoeff();
343 real RRBF = elevationInfo.RBFRadius * std::sqrt(RMax);
344 boundInterpR(iN, 0) = RRBF;
345 DNDS_assert(RRBF > 0);
346 Eigen::Vector<real, Eigen::Dynamic> outDists = outDistancesSqr.array().sqrt() * (1. / RRBF);
347 auto fBasis = RBF::FRBFBasis(outDists, elevationInfo.kernel);
348 for (index in2n = 0; in2n < nFound; in2n++)
349 M.insert(iN, idxFound[in2n]) = fBasis(in2n);
350
351 f(iN, EigenAll) = (*boundInterpVal.son)[iN].transpose();
352 }
353
354 log() << "RBF assembled: " << boundInterpCoo.son->Size() << std::endl;
355 MatrixXR coefs;
356 M.makeCompressed();
357 Eigen::SparseLU<Eigen::SparseMatrix<real>, Eigen::COLAMDOrdering<int>> LUSolver;
358 LUSolver.analyzePattern(M);
359 LUSolver.factorize(M);
360 coefs = LUSolver.solve(f);
361 for (index iN = 0; iN < boundInterpGlobSize; iN++)
362 boundInterpCoef[iN] = coefs(iN, EigenAll).transpose();
363 }
364
365 boundInterpCoef.father->createGlobalMapping();
366 boundInterpCoef.TransAttach();
367 boundInterpCoef.trans.createGhostMapping(boundInterpPullIdx);
368 boundInterpCoef.trans.createMPITypes();
369 boundInterpCoef.trans.pullOnce();
370
371 boundInterpR.TransAttach();
372 boundInterpR.trans.BorrowGGIndexing(boundInterpCoef.trans);
373 boundInterpR.trans.createMPITypes();
374 boundInterpR.trans.pullOnce();
375
376 if (mpi.rank == mRank)
377 {
378 for (index iN = 0; iN < boundInterpGlobSize; iN++)
379 {
380 std::cout << "coef: " << iN << " " << boundInterpCoef.son->operator[](iN).transpose() << std::endl;
381 }
382 }
383
384 PointCloudKDTree insidePts;
385 insidePts.pts.reserve(coords.father->Size() - nodesBoundInterpolated.size());
386 std::vector<index> insideNodes;
387 insideNodes.reserve(insidePts.pts.size());
388 for (index iN = 0; iN < coords.father->Size(); iN++)
389 if (!nodesBoundInterpolated.count(iN))
390 insidePts.pts.push_back(coords[iN]), insideNodes.push_back(iN), coordsElevDisp[iN].setZero();
391 kdtree_tcoo nodesDstTree(3, insidePts);
392 for (index iN = 0; iN < boundInterpGlobSize; iN++)
393 {
394
395 tPoint cooC = (*boundInterpCoo.son)[iN];
396#if NANOFLANN_VERSION < 0x150
397 std::vector<std::pair<DNDS::index, DNDS::real>> IndicesDists;
398#else
399 std::vector<nanoflann::ResultItem<DNDS::index, DNDS::real>> IndicesDists;
400#endif
401 IndicesDists.reserve(elevationInfo.nSearch * 5);
402 real RRBF = boundInterpR.son->operator()(iN, 0);
403#if NANOFLANN_VERSION < 0x150
404 nanoflann::SearchParams params{}; // default params
405#else
406 nanoflann::SearchParameters params{};
407#endif
408 index nFound = nodesDstTree.radiusSearch(cooC.data(), RRBF, IndicesDists, params);
409 Eigen::Vector<real, Eigen::Dynamic> outDists;
410 outDists.resize(IndicesDists.size());
411 DNDS_assert(RRBF > 0);
412 for (index i = 0; i < IndicesDists.size(); i++)
413 outDists[i] = std::sqrt(IndicesDists[i].second) / RRBF;
414
415 auto fBasis = RBF::FRBFBasis(outDists, elevationInfo.kernel);
416 for (index i = 0; i < IndicesDists.size(); i++)
417 {
418 coordsElevDisp[insideNodes[IndicesDists[i].first]] +=
419 fBasis(i, 0) *
420 boundInterpCoef.son->operator[](iN);
421 }
422 }
423
424 // for (index iN = 0; iN < coords.father->Size(); iN++)
425 // {
426 // if (nodesBoundInterpolated.count(iN))
427 // continue;
428
429 // index nFind = elevationInfo.nSearch * 5; // for safety
430 // tPoint cooC = coords[iN];
431
432 // std::vector<index> idxFound;
433 // Eigen::Vector<real, Eigen::Dynamic> outDistancesSqr;
434 // idxFound.resize(nFind);
435 // outDistancesSqr.resize(nFind);
436
437 // index nFound = bndInterpTree.knnSearch(cooC.data(), nFind, idxFound.data(), outDistancesSqr.data());
438 // DNDS_assert(nFound >= 1);
439 // idxFound.resize(nFound);
440 // outDistancesSqr.resize(nFound);
441
442 // Eigen::Vector<real, Eigen::Dynamic> foundRRBFs;
443 // foundRRBFs.resize(nFound);
444 // tSmallCoords coefsC;
445 // coefsC.resize(3, nFound);
446 // for (index iF = 0; iF < nFound; iF++)
447 // foundRRBFs[iF] = (*boundInterpR.son)(idxFound[iF], 0),
448 // coefsC(EigenAll, iF) = (*boundInterpCoef.son)[idxFound[iF]].transpose();
449 // Eigen::Vector<real, Eigen::Dynamic> outDists = outDistancesSqr.array().sqrt() / foundRRBFs.array();
450 // auto fBasis = RBF::FRBFBasis(outDists, elevationInfo.kernel);
451 // tPoint dispC = coefsC * fBasis;
452 // coordsElevDisp[iN] = dispC;
453 // }
454
455 coordsElevDisp.trans.startPersistentPull();
456 coordsElevDisp.trans.waitPersistentPull();
457
458 for (index iN = nNodeO1; iN < coords.father->Size(); iN++)
459 {
460 // if(dim == 2)
461 // dispO2[iN](2) = 0;
462 // if(dispO2[iN].norm() != 0)
463 // std::cout << dispO2[iN].transpose() << std::endl;
464 if (coordsElevDisp[iN](0) != largeReal)
465 coords[iN] += coordsElevDisp[iN];
466 }
467 coords.trans.pullOnce();
468 }
469
470 void UnstructuredMesh::ElevatedNodesSolveInternalSmoothV1()
471 {
472 DNDS_assert(elevState == Elevation_O1O2);
473 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
474 DNDS_assert(adjFacialState == Adj_PointToLocal);
475 DNDS_assert(adjC2FState == Adj_PointToLocal);
476 DNDS_assert(face2node.father);
477 DNDS_assert(nTotalMoved >= 0);
478 if (!nTotalMoved)
479 {
480 if (mpi.rank == mRank)
481 log() << "UnstructuredMesh === ElevatedNodesSolveInternalSmooth() early exit for no nodes were moved";
482 return;
483 }
484
485 DNDS_MPI_InsertCheck(mpi, "Got node2node");
486
487 coordsElevDisp.trans.initPersistentPull(); // only holds local nodes
488
489 std::unordered_set<index> nodesBoundInterpolated;
490 for (index iN = 0; iN < coords.father->Size(); iN++)
491 {
492 if (coordsElevDisp[iN](0) != largeReal || coordsElevDisp[iN](2) == 2 * largeReal)
493 {
494 nodesBoundInterpolated.insert(iN);
495 }
496 }
497
498 tCoordPair boundInterpCoo;
499 tCoordPair boundInterpVal;
500 boundInterpCoo.InitPair("ElevatedNodesSolveInternalSmoothV1::boundInterpCoo", mpi);
501 boundInterpVal.InitPair("ElevatedNodesSolveInternalSmoothV1::boundInterpVal", mpi);
502
503 boundInterpCoo.father->Resize(nodesBoundInterpolated.size());
504 boundInterpVal.father->Resize(nodesBoundInterpolated.size());
505
506 index top{0};
507 for (auto iN : nodesBoundInterpolated)
508 {
509 boundInterpCoo[top] = coords[iN];
510 boundInterpVal[top] = (coordsElevDisp[iN](0) != largeReal) ? tPoint{coordsElevDisp[iN]} : tPoint::Zero();
511 top++;
512 }
513 std::vector<index> boundInterpPullIdx;
514
515 index boundInterpGlobSize = boundInterpCoo.father->globalSize();
516 index boundInterpLocSize = nodesBoundInterpolated.size();
517 boundInterpPullIdx.resize(boundInterpGlobSize);
518 for (index i = 0; i < boundInterpPullIdx.size(); i++)
519 boundInterpPullIdx[i] = i;
520 boundInterpCoo.father->createGlobalMapping();
521 boundInterpCoo.TransAttach();
522 boundInterpCoo.trans.createGhostMapping(boundInterpPullIdx);
523 boundInterpCoo.trans.createMPITypes();
524 boundInterpCoo.trans.pullOnce();
525
526 boundInterpVal.TransAttach();
527 boundInterpVal.trans.BorrowGGIndexing(boundInterpCoo.trans);
528 boundInterpVal.trans.createMPITypes();
529 boundInterpVal.trans.pullOnce();
530
531 if (mpi.rank == mRank)
532 log() << "RBF set: " << boundInterpCoo.son->Size() << std::endl;
533
534 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
535 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTreeCoordPair>,
536 PointCloudKDTreeCoordPair,
537 3,
538 index>;
539 using kdtree_tcoo = nanoflann::KDTreeSingleIndexAdaptor<
540 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTree>,
541 PointCloudKDTree,
542 3,
543 index>;
544 auto coordsI = PointCloudKDTreeCoordPair(boundInterpCoo.son);
545 kdtree_t bndInterpTree(3, coordsI);
546
547 CoordPairDOF boundInterpCoef, boundInterpCoefRHS;
548 using tScalarPair = DNDS::ArrayPair<DNDS::ParArray<real, 1>>;
549 tScalarPair boundInterpR;
550 std::cout << "HEre-2" << std::endl;
551 boundInterpCoef.InitPair("ElevatedNodesSolveInternalSmoothV1::boundInterpCoef", mpi);
552 boundInterpCoefRHS.InitPair("ElevatedNodesSolveInternalSmoothV1::boundInterpCoefRHS", mpi);
553 boundInterpR.InitPair("ElevatedNodesSolveInternalSmoothV1::boundInterpR", mpi);
554 boundInterpCoef.father->Resize(boundInterpCoo.father->Size());
555 boundInterpCoefRHS.father->Resize(boundInterpCoo.father->Size());
556 boundInterpR.father->Resize(boundInterpCoo.father->Size());
557 std::vector<std::vector<std::pair<index, real>>> MatC;
558 std::vector<index> boundInterpPullingIndexSolving;
559 MatC.resize(boundInterpLocSize);
560 std::cout << "HEre-1" << std::endl;
561 for (index iN = 0; iN < boundInterpLocSize; iN++)
562 {
563 tPoint cooC = (*boundInterpCoo.father)[iN];
564
565 index nFind = elevationInfo.nSearch;
566 std::vector<index> idxFound;
567 Eigen::Vector<real, Eigen::Dynamic> outDistancesSqr;
568
569 /**********************************************/
570 idxFound.resize(nFind);
571 outDistancesSqr.resize(nFind);
572 index nFound = bndInterpTree.knnSearch(cooC.data(), nFind, idxFound.data(), outDistancesSqr.data());
573 DNDS_assert(nFound >= 1);
574 idxFound.resize(nFound);
575 outDistancesSqr.resize(nFound);
576
577 real RMin = std::sqrt(outDistancesSqr.minCoeff());
578 real RMax = std::sqrt(outDistancesSqr.maxCoeff());
579 real RRBF = elevationInfo.RBFRadius * RMax;
580 /***********************************************/
581 // std::vector<std::pair<DNDS::index, DNDS::real>> indicesDists;
582 // indicesDists.reserve(nFind);
583
584 // nanoflann::SearchParams searchParams;
585 // index nFound = bndInterpTree.radiusSearch(cooC.data(), elevationInfo.RBFRadius, indicesDists, searchParams);
586 // idxFound.reserve(nFound);
587 // outDistancesSqr.resize(nFound);
588 // for (auto v : indicesDists)
589 // outDistancesSqr(idxFound.size()) = v.second, idxFound.push_back(v.first);
590 // real RRBF = elevationInfo.RBFRadius;
591 /***********************************************/
592
593 boundInterpR(iN, 0) = RRBF;
594 DNDS_assert(RRBF > 0);
595 Eigen::Vector<real, Eigen::Dynamic> outDists = outDistancesSqr.array().pow(0.5 / elevationInfo.RBFPower) * (1. / std::pow(RRBF, 1. / elevationInfo.RBFPower));
596 auto fBasis = RBF::FRBFBasis(outDists, elevationInfo.kernel);
597 MatC[iN].resize(nFound);
598 for (index in2n = 0; in2n < nFound; in2n++)
599 {
600 MatC[iN][in2n] = std::make_pair(index(idxFound[in2n]), real(fBasis(in2n)));
601 // idxFound[in2n] is a global indexing!
602 auto [search_good, rank, val] = boundInterpCoo.trans.pLGlobalMapping->search(idxFound[in2n]);
603 DNDS_assert(search_good);
604 if (rank != mpi.rank)
605 boundInterpPullingIndexSolving.push_back(idxFound[in2n]);
606 }
607
608 boundInterpCoefRHS[iN] = (*boundInterpVal.father)[iN];
609 boundInterpCoef[iN].setZero(); // init val
610 // if(mpi.rank == 0)
611 // {
612 // std::cout << "mat row at " << cooC.transpose() << "R = " << RRBF << "\n";
613 // for (auto idx : idxFound)
614 // std::cout << idx << " ";
615 // std::cout << "\n"
616 // << fBasis.transpose() << "\n"
617 // << outDists.transpose() << "\n"
618 // << boundInterpCoefRHS[iN].transpose() << std::endl;
619 // }
620 }
621
622 if (false)
623 { // use superlu_dist to solve
624
625 // gridinfo_t grid;
626 // superlu_gridinit(mpi.comm, 1, mpi.size, &grid);
627 // DNDS_assert(grid.iam < mpi.size && grid.iam >= 0);
628
629 // std::vector<double> nzval, b;
630 // std::vector<int_t> colind;
631 // std::vector<int_t> rowptr;
632 // rowptr.resize(boundInterpLocSize + 1);
633 // rowptr[0] = 0;
634 // for (index i = 0; i < boundInterpLocSize; i++)
635 // rowptr[i + 1] = rowptr[i] + MatC[i].size();
636 // colind.resize(rowptr.back());
637 // nzval.resize(rowptr.back());
638 // for (index i = 0; i < boundInterpLocSize; i++)
639 // for (index i2j = 0; i2j < MatC[i].size(); i2j++)
640 // {
641 // colind[rowptr[i] + i2j] = MatC[i][i2j].first;
642 // nzval[rowptr[i] + i2j] = MatC[i][i2j].second;
643 // }
644 // b.resize(boundInterpLocSize * 3);
645 // for (index i = 0; i < boundInterpLocSize; i++)
646 // {
647 // b[boundInterpLocSize * 0 + i] = boundInterpCoefRHS[i](0);
648 // b[boundInterpLocSize * 1 + i] = boundInterpCoefRHS[i](1);
649 // b[boundInterpLocSize * 2 + i] = boundInterpCoefRHS[i](2);
650 // }
651
652 // SuperMatrix Aloc;
653 // dCreate_CompRowLoc_Matrix_dist(
654 // &Aloc, boundInterpGlobSize, boundInterpGlobSize, nzval.size(), boundInterpLocSize,
655 // boundInterpCoo.trans.pLGlobalMapping->operator()(mpi.rank, 0), nzval.data(), colind.data(), rowptr.data(),
656 // SLU_NR_loc, SLU_D, SLU_GE);
657
658 // superlu_dist_options_t options;
659 // set_default_options_dist(&options);
660 // SuperLUStat_t stat;
661 // PStatInit(&stat);
662 // dScalePermstruct_t ScalePermstruct;
663 // dScalePermstructInit(Aloc.nrow, Aloc.ncol, &ScalePermstruct);
664 // dLUstruct_t LUstruct;
665 // dLUstructInit(Aloc.ncol, &LUstruct);
666 // dSOLVEstruct_t SOLVEstruct;
667
668 // tPoint berr;
669 // berr.setZero();
670 // int info;
671
672 // pdgssvx(&options, &Aloc, &ScalePermstruct,
673 // b.data(), boundInterpLocSize, 3,
674 // &grid, &LUstruct, &SOLVEstruct, berr.data(), &stat, &info);
675 // for (index i = 0; i < boundInterpLocSize; i++)
676 // {
677 // boundInterpCoef[i](0) = b[boundInterpLocSize * 0 + i];
678 // boundInterpCoef[i](1) = b[boundInterpLocSize * 1 + i];
679 // boundInterpCoef[i](2) = b[boundInterpLocSize * 2 + i];
680 // }
681 // PStatPrint(&options, &stat, &grid);
682 // DNDS_assert(info == 0);
683
684 // PStatFree(&stat);
685 // // TODO: sanitize with valgrind!
686 // // Destroy_CompRowLoc_Matrix_dist(&Aloc); // use vector for strorage, this causes double free
687 // SUPERLU_FREE(Aloc.Store);
688 // dScalePermstructFree(&ScalePermstruct);
689 // dDestroy_LU(Aloc.ncol, &grid, &LUstruct);
690 // dLUstructFree(&LUstruct);
691 // if (options.SolveInitialized)
692 // {
693 // dSolveFinalize(&options, &SOLVEstruct);
694 // }
695
696 // superlu_gridexit(&grid);
697 }
698 if (true)
699 { // use GMRES to solve
700 std::cout << "HEre0" << std::endl;
701 boundInterpCoef.father->createGlobalMapping();
702 boundInterpCoef.TransAttach();
703 boundInterpCoef.trans.createGhostMapping(boundInterpPullingIndexSolving);
704 boundInterpCoef.trans.createMPITypes();
705 boundInterpCoef.trans.initPersistentPull();
706
707 boundInterpCoefRHS.TransAttach();
708 boundInterpCoefRHS.trans.BorrowGGIndexing(boundInterpCoef.trans);
709 boundInterpCoefRHS.trans.createMPITypes();
710 boundInterpCoefRHS.trans.initPersistentPull();
711
712 std::cout << "HEre1" << std::endl;
713 for (index iN = 0; iN < boundInterpLocSize; iN++)
714 {
715 for (index in2n = 0; in2n < MatC[iN].size(); in2n++)
716 {
717 // MatC[iN][in2n].first;
718 auto [search_good, rank, val] = boundInterpCoef.trans.pLGhostMapping->search_indexAppend(MatC[iN][in2n].first);
719 DNDS_assert(search_good);
720 DNDS_assert(val < boundInterpCoef.Size());
721 MatC[iN][in2n].first = val; // to local
722 }
723 }
724 Eigen::SparseMatrix<real> M(boundInterpLocSize, boundInterpLocSize);
725 M.reserve(Eigen::Vector<index, Eigen::Dynamic>::Constant(boundInterpLocSize, elevationInfo.nSearch));
726 for (index iN = 0; iN < boundInterpLocSize; iN++)
727 {
728 for (index in2n = 0; in2n < MatC[iN].size(); in2n++)
729 {
730 if (MatC[iN][in2n].first < boundInterpLocSize)
731 M.insert(iN, MatC[iN][in2n].first) = MatC[iN][in2n].second;
732 // ,
733 // std::cout << fmt::format("rank {}, inserting {},{}={}", mpi.rank, iN, MatC[iN][in2n].first, MatC[iN][in2n].second) << std::endl;
734 }
735 }
736 Eigen::SparseLU<Eigen::SparseMatrix<real>, Eigen::COLAMDOrdering<int>> LUSolver;
737 if (boundInterpLocSize)
738 {
739 LUSolver.analyzePattern(M);
740 LUSolver.factorize(M); // full LU perconditioner
741 }
742 std::cout << "HEre2" << std::endl;
743
744 Linear::GMRES_LeftPreconditioned<CoordPairDOF> gmres(
745 5,
746 [&](CoordPairDOF &v)
747 {
748 v.InitPair("ElevatedNodesSolveInternalSmoothV1::v", mpi);
749 v.father->Resize(boundInterpCoef.father->Size());
750 v.TransAttach();
751 v.trans.BorrowGGIndexing(boundInterpCoef.trans);
752 v.trans.createMPITypes();
753 v.trans.initPersistentPull();
754 });
755 std::cout << "HEre3" << std::endl;
756 boundInterpCoef.trans.startPersistentPull();
757 boundInterpCoef.trans.waitPersistentPull();
758 boundInterpCoefRHS.trans.startPersistentPull();
759 boundInterpCoefRHS.trans.waitPersistentPull();
760 gmres.solve(
761 [&](CoordPairDOF &x, CoordPairDOF &Ax)
762 {
763 x.trans.startPersistentPull();
764 x.trans.waitPersistentPull();
765 for (index iN = 0; iN < boundInterpLocSize; iN++)
766 {
767 Ax[iN].setZero();
768 for (index in2n = 0; in2n < MatC[iN].size(); in2n++)
769 {
770 Ax[iN] += x[MatC[iN][in2n].first] * MatC[iN][in2n].second;
771 }
772 }
773 Ax.trans.startPersistentPull();
774 Ax.trans.waitPersistentPull();
775 },
776 [&](CoordPairDOF &x, CoordPairDOF &MLX)
777 {
778 x.trans.startPersistentPull();
779 x.trans.waitPersistentPull();
780 Eigen::Matrix<real, Eigen::Dynamic, 3> xVal;
781 xVal.resize(boundInterpLocSize, 3);
782 for (index iN = 0; iN < boundInterpLocSize; iN++)
783 xVal(iN, EigenAll) = x[iN].transpose();
784 Eigen::Matrix<real, Eigen::Dynamic, 3> xValPrec;
785 if (boundInterpLocSize)
786 xValPrec = LUSolver.solve(xVal);
787 for (index iN = 0; iN < boundInterpLocSize; iN++)
788 MLX[iN] = xValPrec(iN, EigenAll).transpose();
789 MLX = x;
790 MLX.trans.startPersistentPull();
791 MLX.trans.waitPersistentPull();
792 },
793 [&](CoordPairDOF &a, CoordPairDOF &b) -> real
794 {
795 return a.dot(b);
796 },
797 boundInterpCoefRHS,
798 boundInterpCoef,
799 elevationInfo.nIter,
800 [&](int iRestart, real res, real resB)
801 {
802 if (mpi.rank == mRank)
803 log() << fmt::format("iRestart {}, res {}, resB {}", iRestart, res, resB) << std::endl;
804 return res < resB * 1e-6;
805 }
806
807 );
808 boundInterpCoef.trans.startPersistentPull();
809 boundInterpCoef.trans.waitPersistentPull();
810 }
811
812 boundInterpCoef.father->createGlobalMapping();
813 boundInterpCoef.TransAttach();
814 boundInterpCoef.trans.createGhostMapping(boundInterpPullIdx);
815 boundInterpCoef.trans.createMPITypes();
816 boundInterpCoef.trans.pullOnce();
817
818 boundInterpR.TransAttach();
819 boundInterpR.trans.BorrowGGIndexing(boundInterpCoo.trans);
820 boundInterpR.trans.createMPITypes();
821 boundInterpR.trans.pullOnce();
822
823 // if (mpi.rank == mRank)
824 // {
825 // for (index iN = 0; iN < boundInterpGlobSize; iN++)
826 // {
827 // std::cout << "pt - coef: R " << boundInterpR.son->operator()(iN, 0)
828 // << " coo " << boundInterpCoo.son->operator[](iN).transpose()
829 // << " -> " << boundInterpCoef.son->operator[](iN).transpose() << std::endl;
830 // }
831 // }
832
833 PointCloudKDTree insidePts;
834 insidePts.pts.reserve(coords.father->Size() - boundInterpLocSize);
835 std::vector<index> insideNodes;
836 insideNodes.reserve(insidePts.pts.size());
837 for (index iN = 0; iN < coords.father->Size(); iN++)
838 if (!nodesBoundInterpolated.count(iN))
839 insidePts.pts.push_back(coords[iN]), insideNodes.push_back(iN), coordsElevDisp[iN].setZero();
840 kdtree_tcoo nodesDstTree(3, insidePts);
841 for (index iN = 0; iN < boundInterpGlobSize; iN++)
842 {
843
844 tPoint cooC = (*boundInterpCoo.son)[iN];
845#if NANOFLANN_VERSION < 0x150
846 std::vector<std::pair<DNDS::index, DNDS::real>> IndicesDists;
847#else
848 std::vector<nanoflann::ResultItem<DNDS::index, DNDS::real>> IndicesDists;
849#endif
850 IndicesDists.reserve(elevationInfo.nSearch * 5);
851 real RRBF = boundInterpR.son->operator()(iN, 0);
852#if NANOFLANN_VERSION < 0x150
853 nanoflann::SearchParams params{}; // default params
854#else
855 nanoflann::SearchParameters params{};
856#endif
857 index nFound = nodesDstTree.radiusSearch(cooC.data(), RRBF, IndicesDists, params);
858 Eigen::Vector<real, Eigen::Dynamic> outDists;
859 outDists.resize(IndicesDists.size());
860 DNDS_assert(RRBF > 0);
861 for (index i = 0; i < IndicesDists.size(); i++)
862 outDists[i] = std::sqrt(IndicesDists[i].second) / std::pow(RRBF, 1. / elevationInfo.RBFPower);
863
864 auto fBasis = RBF::FRBFBasis(outDists, elevationInfo.kernel);
865 for (index i = 0; i < IndicesDists.size(); i++)
866 {
867 coordsElevDisp[insideNodes[IndicesDists[i].first]] +=
868 fBasis(i, 0) *
869 boundInterpCoef.son->operator[](iN);
870 }
871 }
872
873 // for (index iN = 0; iN < coords.father->Size(); iN++)
874 // {
875 // if (nodesBoundInterpolated.count(iN))
876 // continue;
877
878 // index nFind = elevationInfo.nSearch * 5; // for safety
879 // tPoint cooC = coords[iN];
880
881 // std::vector<index> idxFound;
882 // Eigen::Vector<real, Eigen::Dynamic> outDistancesSqr;
883 // idxFound.resize(nFind);
884 // outDistancesSqr.resize(nFind);
885
886 // index nFound = bndInterpTree.knnSearch(cooC.data(), nFind, idxFound.data(), outDistancesSqr.data());
887 // DNDS_assert(nFound >= 1);
888 // idxFound.resize(nFound);
889 // outDistancesSqr.resize(nFound);
890
891 // Eigen::Vector<real, Eigen::Dynamic> foundRRBFs;
892 // foundRRBFs.resize(nFound);
893 // tSmallCoords coefsC;
894 // coefsC.resize(3, nFound);
895 // for (index iF = 0; iF < nFound; iF++)
896 // foundRRBFs[iF] = (*boundInterpR.son)(idxFound[iF], 0),
897 // coefsC(EigenAll, iF) = (*boundInterpCoef.son)[idxFound[iF]].transpose();
898 // Eigen::Vector<real, Eigen::Dynamic> outDists = outDistancesSqr.array().sqrt() / foundRRBFs.array();
899 // auto fBasis = RBF::FRBFBasis(outDists, elevationInfo.kernel);
900 // tPoint dispC = coefsC * fBasis;
901 // coordsElevDisp[iN] = dispC;
902 // }
903
904 coordsElevDisp.trans.startPersistentPull();
905 coordsElevDisp.trans.waitPersistentPull();
906
907 for (index iN = 0; iN < coords.father->Size(); iN++)
908 {
909 // if(dim == 2)
910 // dispO2[iN](2) = 0;
911 // if(dispO2[iN].norm() != 0)
912 // std::cout << dispO2[iN].transpose() << std::endl;
913 if (coordsElevDisp[iN](0) != largeReal)
914 coords[iN] += coordsElevDisp[iN];
915 }
916 coords.trans.pullOnce();
917 }
918
919 void UnstructuredMesh::ElevatedNodesSolveInternalSmoothV2()
920 {
921 DNDS_assert(elevState == Elevation_O1O2);
922 DNDS_assert(adjPrimaryState == Adj_PointToLocal);
923 DNDS_assert(adjFacialState == Adj_PointToLocal);
924 DNDS_assert(adjC2FState == Adj_PointToLocal);
925 DNDS_assert(face2node.father);
926 DNDS_assert(nTotalMoved >= 0);
927 if (!nTotalMoved)
928 {
929 if (mpi.rank == mRank)
930 log() << "UnstructuredMesh === ElevatedNodesSolveInternalSmooth() early exit for no nodes were moved";
931 return;
932 }
933
934 std::unordered_set<index> nodesBoundInterpolated;
935 for (index iN = 0; iN < coords.father->Size(); iN++)
936 {
937 if (coordsElevDisp[iN](0) != largeReal || coordsElevDisp[iN](2) == 2 * largeReal)
938 {
939 nodesBoundInterpolated.insert(iN);
940 }
941 }
942
943 tCoordPair boundInterpCoo;
944 tCoordPair boundInterpVal;
945 boundInterpCoo.InitPair("ElevatedNodesSolveInternalSmoothV2::boundInterpCoo", mpi);
946 boundInterpVal.InitPair("ElevatedNodesSolveInternalSmoothV2::boundInterpVal", mpi);
947
948 boundInterpCoo.father->Resize(nodesBoundInterpolated.size());
949 boundInterpVal.father->Resize(nodesBoundInterpolated.size());
950
951 index top{0};
952 for (auto iN : nodesBoundInterpolated)
953 {
954 boundInterpCoo[top] = coords[iN];
955 boundInterpVal[top] = (coordsElevDisp[iN](0) != largeReal) ? tPoint{coordsElevDisp[iN]} : tPoint::Zero();
956 top++;
957 }
958 std::vector<index> boundInterpPullIdx;
959
960 index boundInterpGlobSize = boundInterpCoo.father->globalSize();
961 index boundInterpLocSize = nodesBoundInterpolated.size();
962 boundInterpPullIdx.resize(boundInterpGlobSize);
963 for (index i = 0; i < boundInterpPullIdx.size(); i++)
964 boundInterpPullIdx[i] = i;
965 boundInterpCoo.father->createGlobalMapping();
966 boundInterpCoo.TransAttach();
967 boundInterpCoo.trans.createGhostMapping(boundInterpPullIdx);
968 boundInterpCoo.trans.createMPITypes();
969 boundInterpCoo.trans.pullOnce();
970
971 boundInterpVal.TransAttach();
972 boundInterpVal.trans.BorrowGGIndexing(boundInterpCoo.trans);
973 boundInterpVal.trans.createMPITypes();
974 boundInterpVal.trans.pullOnce();
975
976 std::vector<std::unordered_set<index>> node2nodeV;
977 node2nodeV.resize(coords.father->Size());
978 for (index iCell = 0; iCell < cell2node.Size(); iCell++)
979 {
980 for (auto iN : cell2node[iCell])
981 for (auto iNOther : cell2node[iCell])
982 {
983 if (iN == iNOther)
984 continue;
985 if (iN < coords.father->Size())
986 {
987 node2nodeV[iN].insert(iNOther);
988 }
989 }
990 }
991
992 if (mpi.rank == mRank)
993 log() << "RBF set: " << boundInterpCoo.son->Size() << std::endl;
994
995 using kdtree_t = nanoflann::KDTreeSingleIndexAdaptor<
996 nanoflann::L2_Simple_Adaptor<real, PointCloudKDTreeCoordPair>,
997 PointCloudKDTreeCoordPair,
998 3,
999 index>;
1000 auto coordsI = PointCloudKDTreeCoordPair(boundInterpCoo.son);
1001 kdtree_t bndInterpTree(3, coordsI);
1002
1003 DNDS_MPI_InsertCheck(mpi, "Got node2node");
1004 // std::cout << "hereXXXXXXXXXXXXXXX 0" << std::endl;
1005
1006 CoordPairDOF dispO2, bO2, dispO2New, dispO2PrecRHS, dispO2Inc, dispO2IncLimited;
1007 auto initDOF = [&](CoordPairDOF &dispO2)
1008 {
1009 dispO2.InitPair("ElevatedNodesSolveInternalSmoothV2::dispO2", mpi);
1010 dispO2.father->Resize(coords.father->Size()); // O1 part is unused
1011 dispO2.TransAttach();
1012 dispO2.trans.BorrowGGIndexing(coords.trans); // should be subset of coord?
1013 dispO2.trans.createMPITypes();
1014 dispO2.trans.initPersistentPull();
1015 };
1016 initDOF(dispO2);
1017 initDOF(dispO2Inc);
1018 initDOF(dispO2New);
1019 initDOF(dispO2IncLimited);
1020 initDOF(bO2);
1021 initDOF(dispO2PrecRHS);
1022
1023 struct MatElem
1024 {
1025 index j;
1026 Eigen::Matrix<real, 3, 3> m;
1027 };
1028 std::vector<std::vector<MatElem>> A; // stiff matrixTo O2 coords
1029 A.resize(coords.father->Size());
1030 for (index iN = 0; iN < coords.father->Size(); iN++)
1031 {
1032 A[iN].resize(node2nodeV[iN].size() + 1);
1033 A[iN][0].j = iN;
1034 int curRowEnd = 1;
1035 for (auto iNOther : node2nodeV[iN])
1036 {
1037 A[iN][curRowEnd++].j = iNOther;
1038 }
1039 for (auto &ME : A[iN])
1040 ME.m.setConstant(0);
1041 bO2[iN].setZero();
1042 dispO2[iN].setZero();
1043 if (coordsElevDisp[iN](0) != largeReal)
1044 dispO2[iN] = coordsElevDisp[iN];
1045 if (coordsElevDisp[iN](2) == 2 * largeReal || coordsElevDisp[iN](2) == 3 * largeReal)
1046 dispO2[iN].setZero();
1047 }
1048 dispO2.trans.startPersistentPull();
1049 dispO2.trans.waitPersistentPull();
1050
1051 auto AssembleRHSMat = [&](tCoordPair &uCur)
1052 {
1053 index nFix{0}, nFixB{0};
1054 if (false) // fem method
1055 for (index iCell = 0; iCell < cell2node.Size(); iCell++)
1056 {
1057 auto c2n = cell2node[iCell];
1058 rowsize nnLoc = c2n.size();
1059 tSmallCoords coordsC;
1060 GetCoordsOnCell(iCell, coordsC);
1061 tPoint cellCent = coordsC.rowwise().mean();
1062 real wDist{0};
1063 {
1064 std::vector<real> sqrDists(1);
1065 std::vector<index> idxFound(1);
1066 auto nFound = bndInterpTree.knnSearch(cellCent.data(), 1, idxFound.data(), sqrDists.data());
1067 DNDS_assert(nFound == 1);
1068 wDist = std::sqrt(sqrDists[0]);
1069 }
1070
1071 MatrixXR ALoc, mLoc;
1072 ALoc.setZero(3 * nnLoc, 3 * nnLoc + 1);
1073 mLoc.resize(6, 3 * nnLoc);
1074 auto eCell = GetCellElement(iCell);
1075 auto qCell = Elem::Quadrature{eCell, 2}; // for O2 FEM
1076 qCell.Integration(
1077 ALoc,
1078 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
1079 {
1080 // J = dx_i/dxii_j, Jinv = dxii_i/dx_j
1081 tJacobi J = Elem::ShapeJacobianCoordD01Nj(coordsC, DiNj);
1082 tJacobi JInv = tJacobi::Identity();
1083 real JDet;
1084 if (dim == 2)
1085 JDet = J(EigenAll, 0).cross(J(EigenAll, 1)).stableNorm();
1086 else
1087 JDet = J.determinant();
1088 DNDS_assert(JDet > 0);
1089 if (dim == 2)
1090 JInv({0, 1}, {0, 1}) = J({0, 1}, {0, 1}).inverse().eval();
1091 else
1092 JInv = J.inverse().eval();
1093 // JInv.transpose() * DiNj({1,2,3}, EigenAll) * u (size = nnode,3) => du_j/dxi
1094 tSmallCoords m3 = JInv.transpose() * DiNj({1, 2, 3}, EigenAll);
1095 mLoc(0, Eigen::seq(nnLoc * 0, nnLoc * 0 + nnLoc - 1)) = m3(0, EigenAll);
1096 mLoc(1, Eigen::seq(nnLoc * 1, nnLoc * 1 + nnLoc - 1)) = m3(1, EigenAll);
1097 mLoc(2, Eigen::seq(nnLoc * 2, nnLoc * 2 + nnLoc - 1)) = m3(2, EigenAll);
1098
1099 mLoc(3, Eigen::seq(nnLoc * 1, nnLoc * 1 + nnLoc - 1)) = 0.5 * m3(2, EigenAll);
1100 mLoc(3, Eigen::seq(nnLoc * 2, nnLoc * 2 + nnLoc - 1)) = 0.5 * m3(1, EigenAll);
1101 mLoc(4, Eigen::seq(nnLoc * 2, nnLoc * 2 + nnLoc - 1)) = 0.5 * m3(0, EigenAll);
1102 mLoc(4, Eigen::seq(nnLoc * 0, nnLoc * 0 + nnLoc - 1)) = 0.5 * m3(2, EigenAll);
1103 mLoc(5, Eigen::seq(nnLoc * 0, nnLoc * 0 + nnLoc - 1)) = 0.5 * m3(1, EigenAll);
1104 mLoc(5, Eigen::seq(nnLoc * 1, nnLoc * 1 + nnLoc - 1)) = 0.5 * m3(0, EigenAll);
1105
1106 Eigen::Matrix<real, 6, 6> DStruct;
1107 DStruct.setIdentity();
1108 // real lam = std::pow((wDist / 1e-6), -1);
1109 real lam = 1;
1110 real muu = 1;
1111 DStruct *= muu;
1112 DStruct({0, 1, 2}, {0, 1, 2}) += muu * tJacobi::Identity();
1113 DStruct({0, 1, 2}, {0, 1, 2}).array() += lam;
1114 // DStruct *= (1. / JDet);
1115 // DStruct *= std::pow((wDist / 1e-6), -10);
1116
1117 vInc.resizeLike(ALoc);
1118 vInc(EigenAll, Eigen::seq(0, 3 * nnLoc - 1)) = mLoc.transpose() * DStruct * mLoc;
1119 vInc(EigenAll, EigenLast).setZero(); // no node force, thus zero
1120 vInc *= JDet;
1121 });
1122 // for (int ic2n = 0; ic2n < nnLoc; ic2n++)
1123 // for (int jc2n = 0; jc2n < nnLoc; jc2n++)
1124 // {
1125
1126 // if (ic2n == jc2n)
1127 // continue;
1128 // real lBA = (coordsC(EigenAll, jc2n) - coordsC(EigenAll, ic2n)).norm();
1129 // tPoint nBA = (coordsC(EigenAll, jc2n) - coordsC(EigenAll, ic2n)).normalized();
1130 // // real K = std::pow(lBA, 0);
1131 // real K = std::pow((wDist / elevationInfo.refDWall), -1);
1132 // tJacobi mnBA = nBA * nBA.transpose() * K;
1133 // for (int iii = 0; iii < 3; iii++)
1134 // for (int jjj = 0; jjj < 3; jjj++)
1135 // {
1136 // ALoc(iii * nnLoc + ic2n, jjj * nnLoc + ic2n) += mnBA(iii, jjj);
1137 // ALoc(iii * nnLoc + ic2n, jjj * nnLoc + jc2n) -= mnBA(iii, jjj);
1138 // }
1139 //}
1140 if (dim != 3)
1141 {
1142 ALoc(Eigen::seq(2 * nnLoc, 3 * nnLoc - 1), Eigen::seq(2 * nnLoc, 3 * nnLoc - 1)).setZero();
1143 ALoc(EigenAll, Eigen::seq(2 * nnLoc, 3 * nnLoc - 1)).setZero();
1144 ALoc(Eigen::seq(2 * nnLoc, 3 * nnLoc - 1), EigenAll).setZero();
1145 }
1146
1147 for (rowsize ic2n = 0; ic2n < nnLoc; ic2n++)
1148 {
1149 index iN = c2n[ic2n];
1150 if (iN >= coords.father->Size())
1151 continue; // not local row
1152 for (rowsize jc2n = 0; jc2n < nnLoc; jc2n++)
1153 {
1154 index jN = c2n[jc2n];
1155 int nMatrixFound = 0;
1156 for (auto &ME : A[iN])
1157 if (ME.j == jN)
1158 {
1159 nMatrixFound++;
1160 ME.m(0, 0) += ALoc(0 * nnLoc + ic2n, 0 * nnLoc + jc2n);
1161 ME.m(0, 1) += ALoc(0 * nnLoc + ic2n, 1 * nnLoc + jc2n);
1162 ME.m(0, 2) += ALoc(0 * nnLoc + ic2n, 2 * nnLoc + jc2n);
1163 ME.m(1, 0) += ALoc(1 * nnLoc + ic2n, 0 * nnLoc + jc2n);
1164 ME.m(1, 1) += ALoc(1 * nnLoc + ic2n, 1 * nnLoc + jc2n);
1165 ME.m(1, 2) += ALoc(1 * nnLoc + ic2n, 2 * nnLoc + jc2n);
1166 ME.m(2, 0) += ALoc(2 * nnLoc + ic2n, 0 * nnLoc + jc2n);
1167 ME.m(2, 1) += ALoc(2 * nnLoc + ic2n, 1 * nnLoc + jc2n);
1168 ME.m(2, 2) += ALoc(2 * nnLoc + ic2n, 2 * nnLoc + jc2n);
1169 if (isPeriodic)
1170 {
1171 auto ipbi = cell2nodePbi(iCell, ic2n);
1172 auto jpbi = cell2nodePbi(iCell, jc2n);
1173 tJacobi iTrans = periodicInfo.GetVectorByBits<3, 3>(tJacobi::Identity(), ipbi);
1174 tJacobi jTrans = periodicInfo.GetVectorByBits<3, 3>(tJacobi::Identity(), jpbi);
1175 ME.m = jTrans.transpose() * ME.m * iTrans;
1176 }
1177 }
1178 DNDS_assert(nMatrixFound == 1);
1179 }
1180 bO2[iN](0) += ALoc(0 * nnLoc + ic2n, EigenLast);
1181 bO2[iN](1) += ALoc(1 * nnLoc + ic2n, EigenLast);
1182 bO2[iN](2) += ALoc(2 * nnLoc + ic2n, EigenLast); // last col is b
1183 if (coordsElevDisp[iN](0) != largeReal)
1184 {
1185 nFix++;
1186 real nDiag = 0;
1187 for (auto &ME : A[iN])
1188 if (ME.j == iN)
1189 nDiag = ME.m.array().abs().maxCoeff();
1190 bO2[iN] = coordsElevDisp[iN] * nDiag;
1191 // bO2[iN].setZero();
1192 // if (coords[iN].norm() < 2)
1193 // bO2[iN](1) = sqr(coords[iN](0)) * 0.1 * nDiag;
1194 for (auto &ME : A[iN])
1195 if (ME.j == iN)
1196 ME.m = tJacobi::Identity() * nDiag;
1197 else
1198 ME.m.setZero();
1199 }
1200 if (coordsElevDisp[iN](2) == 2 * largeReal || coordsElevDisp[iN](2) == 3 * largeReal)
1201 {
1202 nFixB++;
1203 real nDiag = 0;
1204 for (auto &ME : A[iN])
1205 if (ME.j == iN)
1206 nDiag = ME.m.array().abs().maxCoeff();
1207 bO2[iN].setZero();
1208 // bO2[iN].setZero();
1209 // if (coords[iN].norm() < 2)
1210 // bO2[iN](1) = sqr(coords[iN](0)) * 0.1 * nDiag;
1211 for (auto &ME : A[iN])
1212 if (ME.j == iN)
1213 ME.m = tJacobi::Identity() * nDiag;
1214 else
1215 ME.m.setZero();
1216 }
1217 if (dim != 3)
1218 {
1219 real nDiag = 0;
1220 for (auto &ME : A[iN])
1221 if (ME.j == iN)
1222 nDiag = ME.m.array().abs().maxCoeff();
1223 for (auto &ME : A[iN])
1224 if (ME.j == iN)
1225 ME.m(2, 2) += nDiag; // the z-direction need a constraint
1226 }
1227 }
1228 }
1229
1230 if (true) // bisect fem method
1231 for (index iCell = 0; iCell < cell2node.Size(); iCell++)
1232 {
1233 auto c2n = cell2node[iCell];
1234 rowsize nnLoc = c2n.size();
1235 SmallCoordsAsVector coordsC, coordsCu;
1236 GetCoordsOnCell(iCell, coordsC);
1237 GetCoordsOnCell(iCell, coordsCu, uCur);
1238 tPoint cellCent = coordsC.rowwise().mean();
1239 real wDist{0};
1240 {
1241 std::vector<real> sqrDists(1);
1242 std::vector<index> idxFound(1);
1243 auto nFound = bndInterpTree.knnSearch(cellCent.data(), 1, idxFound.data(), sqrDists.data());
1244 DNDS_assert(nFound == 1);
1245 wDist = std::sqrt(sqrDists[0]);
1246 }
1247
1248 MatrixXR ALoc;
1249 ALoc.setZero(3 * nnLoc, 3 * nnLoc + 1);
1250 auto localNodeIdx = Eigen::ArrayXi::LinSpaced(c2n.size(), 0, c2n.size() - 1);
1251
1252 auto eCell = GetCellElement(iCell);
1253 int nBi = eCell.GetO2NumBisect();
1254 int iBiVariant = GetO2ElemBisectVariant(eCell, coordsC);
1255 for (int iBi = 0; iBi < nBi; iBi++)
1256 {
1257 auto eCellSub = eCell.ObtainO2BisectElem(iBi);
1258 Eigen::ArrayXi c2nSubLocal;
1259 c2nSubLocal.resize(eCellSub.GetNumNodes());
1260 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, localNodeIdx, c2nSubLocal);
1261 SmallCoordsAsVector coordsCSub, coordsCuSub;
1262 coordsCSub.resize(3, eCellSub.GetNumNodes());
1263 coordsCuSub.resize(3, eCellSub.GetNumNodes());
1264 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsC, coordsCSub);
1265 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsCu, coordsCuSub);
1266 auto qCellSub = Elem::Quadrature{eCellSub, 6}; // for O1 FEM
1267 auto nnLocSub = rowsize(c2nSubLocal.size());
1268 MatrixXR ALocSub, mLoc;
1269 ALocSub.setZero(3 * nnLocSub, 3 * nnLocSub + 1);
1270 mLoc.resize(6, 3 * nnLocSub);
1271 Eigen::ArrayXi c2nSubLocal3;
1272 c2nSubLocal3.resize(c2nSubLocal.size() * 3);
1273 c2nSubLocal3(Eigen::seq(c2nSubLocal.size() * 0, c2nSubLocal.size() * 1 - 1)) = c2nSubLocal + nnLoc * 0;
1274 c2nSubLocal3(Eigen::seq(c2nSubLocal.size() * 1, c2nSubLocal.size() * 2 - 1)) = c2nSubLocal + nnLoc * 1;
1275 c2nSubLocal3(Eigen::seq(c2nSubLocal.size() * 2, c2nSubLocal.size() * 3 - 1)) = c2nSubLocal + nnLoc * 2;
1276
1277 qCellSub.Integration(
1278 ALocSub,
1279 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
1280 {
1281 // J = dx_i/dxii_j, Jinv = dxii_i/dx_j
1282 auto getJ = [&](auto coos)
1283 {
1284 tJacobi J = Elem::ShapeJacobianCoordD01Nj(coos, DiNj);
1285 tJacobi JInv = tJacobi::Identity();
1286 real JDet;
1287 if (dim == 2)
1288 JDet = J(EigenAll, 0).cross(J(EigenAll, 1)).stableNorm();
1289 else
1290 JDet = J.determinant();
1291 if (dim == 2)
1292 JInv({0, 1}, {0, 1}) = J({0, 1}, {0, 1}).inverse().eval();
1293 else
1294 JInv = J.inverse().eval();
1295 return std::make_tuple(J, JInv, JDet);
1296 };
1297 auto [J, JInv, JDet] = getJ(coordsCSub);
1298 auto [Jn, JInvn, JDetn] = getJ(coordsCSub + coordsCuSub);
1299 DNDS_assert(JDet > 0);
1300 // JInv.transpose() * DiNj({1,2,3}, EigenAll) * u (size = nnode,3) => du_j/dxi
1301 tSmallCoords m3 = JInv.transpose() * DiNj({1, 2, 3}, EigenAll);
1302 mLoc.setZero();
1303 mLoc(0, Eigen::seq(nnLocSub * 0, nnLocSub * 0 + nnLocSub - 1)) = m3(0, EigenAll);
1304 mLoc(1, Eigen::seq(nnLocSub * 1, nnLocSub * 1 + nnLocSub - 1)) = m3(1, EigenAll);
1305 mLoc(2, Eigen::seq(nnLocSub * 2, nnLocSub * 2 + nnLocSub - 1)) = m3(2, EigenAll);
1306
1307 mLoc(3, Eigen::seq(nnLocSub * 1, nnLocSub * 1 + nnLocSub - 1)) = 0.5 * m3(2, EigenAll);
1308 mLoc(3, Eigen::seq(nnLocSub * 2, nnLocSub * 2 + nnLocSub - 1)) = 0.5 * m3(1, EigenAll);
1309 mLoc(4, Eigen::seq(nnLocSub * 2, nnLocSub * 2 + nnLocSub - 1)) = 0.5 * m3(0, EigenAll);
1310 mLoc(4, Eigen::seq(nnLocSub * 0, nnLocSub * 0 + nnLocSub - 1)) = 0.5 * m3(2, EigenAll);
1311 mLoc(5, Eigen::seq(nnLocSub * 0, nnLocSub * 0 + nnLocSub - 1)) = 0.5 * m3(1, EigenAll);
1312 mLoc(5, Eigen::seq(nnLocSub * 1, nnLocSub * 1 + nnLocSub - 1)) = 0.5 * m3(0, EigenAll);
1313
1314 Eigen::Matrix<real, 6, 6> DStruct;
1315 DStruct.setIdentity();
1316 // real lam = std::pow((wDist / 1e-6), -1);
1317 real lam = 1;
1318 real muu = 100;
1319 // lam *= std::pow(JDetn / JDet, -1) + 1;
1320 DStruct *= muu;
1321 DStruct({0, 1, 2}, {0, 1, 2}) += muu * tJacobi::Identity();
1322 DStruct({0, 1, 2}, {0, 1, 2}).array() += lam;
1323 DStruct *= std::pow(JDetn / JDet, -0.0) + 1;
1324 // DStruct *= (1. / JDet);
1325 // DStruct *= std::pow((wDist / 1e-6), -10);
1326 // DStruct *= std::pow((wDist / elevationInfo.refDWall), -1);
1327
1328 vInc.resizeLike(ALocSub);
1329 vInc(EigenAll, Eigen::seq(0, 3 * nnLocSub - 1)) = mLoc.transpose() * DStruct * mLoc;
1330 vInc(EigenAll, EigenLast).setZero(); // no node force, thus zero
1331 vInc *= JDet;
1332 });
1333 ALoc(c2nSubLocal3, c2nSubLocal3) += ALocSub(EigenAll, Eigen::seq(0, EigenLast - 1));
1334 ALoc(c2nSubLocal3, EigenLast) += ALocSub(EigenAll, EigenLast);
1335 }
1336
1337 if (dim != 3)
1338 {
1339 ALoc(Eigen::seq(2 * nnLoc, 3 * nnLoc - 1), Eigen::seq(2 * nnLoc, 3 * nnLoc - 1)).setZero();
1340 ALoc(EigenAll, Eigen::seq(2 * nnLoc, 3 * nnLoc - 1)).setZero();
1341 ALoc(Eigen::seq(2 * nnLoc, 3 * nnLoc - 1), EigenAll).setZero();
1342 }
1343
1344 for (rowsize ic2n = 0; ic2n < nnLoc; ic2n++)
1345 {
1346 index iN = c2n[ic2n];
1347 if (iN >= coords.father->Size())
1348 continue; // not local row
1349 for (rowsize jc2n = 0; jc2n < nnLoc; jc2n++)
1350 {
1351 index jN = c2n[jc2n];
1352 int nMatrixFound = 0;
1353 for (auto &ME : A[iN])
1354 if (ME.j == jN)
1355 {
1356 nMatrixFound++;
1357 ME.m(0, 0) += ALoc(0 * nnLoc + ic2n, 0 * nnLoc + jc2n);
1358 ME.m(0, 1) += ALoc(0 * nnLoc + ic2n, 1 * nnLoc + jc2n);
1359 ME.m(0, 2) += ALoc(0 * nnLoc + ic2n, 2 * nnLoc + jc2n);
1360 ME.m(1, 0) += ALoc(1 * nnLoc + ic2n, 0 * nnLoc + jc2n);
1361 ME.m(1, 1) += ALoc(1 * nnLoc + ic2n, 1 * nnLoc + jc2n);
1362 ME.m(1, 2) += ALoc(1 * nnLoc + ic2n, 2 * nnLoc + jc2n);
1363 ME.m(2, 0) += ALoc(2 * nnLoc + ic2n, 0 * nnLoc + jc2n);
1364 ME.m(2, 1) += ALoc(2 * nnLoc + ic2n, 1 * nnLoc + jc2n);
1365 ME.m(2, 2) += ALoc(2 * nnLoc + ic2n, 2 * nnLoc + jc2n);
1366 if (isPeriodic)
1367 {
1368 auto ipbi = cell2nodePbi(iCell, ic2n);
1369 auto jpbi = cell2nodePbi(iCell, jc2n);
1370 tJacobi iTrans = periodicInfo.GetVectorByBits<3, 3>(tJacobi::Identity(), ipbi);
1371 tJacobi jTrans = periodicInfo.GetVectorByBits<3, 3>(tJacobi::Identity(), jpbi);
1372 ME.m = jTrans.transpose() * ME.m * iTrans;
1373 }
1374 }
1375 DNDS_assert(nMatrixFound == 1);
1376 }
1377 bO2[iN](0) += ALoc(0 * nnLoc + ic2n, EigenLast);
1378 bO2[iN](1) += ALoc(1 * nnLoc + ic2n, EigenLast);
1379 bO2[iN](2) += ALoc(2 * nnLoc + ic2n, EigenLast); // last col is b
1380 if (coordsElevDisp[iN](0) != largeReal)
1381 {
1382 nFix++;
1383 real nDiag = 0;
1384 for (auto &ME : A[iN])
1385 if (ME.j == iN)
1386 nDiag = ME.m.array().abs().maxCoeff();
1387 bO2[iN] = coordsElevDisp[iN] * nDiag;
1388 // bO2[iN].setZero();
1389 // if (coords[iN].norm() < 2)
1390 // bO2[iN](1) = sqr(coords[iN](0)) * 0.1 * nDiag;
1391 for (auto &ME : A[iN])
1392 if (ME.j == iN)
1393 ME.m = tJacobi::Identity() * nDiag;
1394 else
1395 ME.m.setZero();
1396 }
1397 if (coordsElevDisp[iN](2) == 2 * largeReal || coordsElevDisp[iN](2) == 3 * largeReal || (iN < nNodeO1 && false))
1398 {
1399 nFixB++;
1400 real nDiag = 0;
1401 for (auto &ME : A[iN])
1402 if (ME.j == iN)
1403 nDiag = ME.m.array().abs().maxCoeff();
1404 bO2[iN].setZero();
1405 // bO2[iN].setZero();
1406 // if (coords[iN].norm() < 2)
1407 // bO2[iN](1) = sqr(coords[iN](0)) * 0.1 * nDiag;
1408 for (auto &ME : A[iN])
1409 if (ME.j == iN)
1410 ME.m = tJacobi::Identity() * nDiag;
1411 else
1412 ME.m.setZero();
1413 }
1414 if (dim != 3)
1415 {
1416 real nDiag = 0;
1417 for (auto &ME : A[iN])
1418 if (ME.j == iN)
1419 nDiag = ME.m.array().abs().maxCoeff();
1420 for (auto &ME : A[iN])
1421 if (ME.j == iN)
1422 ME.m(2, 2) += nDiag; // the z-direction need a constraint
1423 }
1424 }
1425 }
1426
1427 if (false)
1428 for (index iN = 0; iN < A.size(); iN++)
1429 {
1430 for (auto jN : node2nodeV[iN])
1431 {
1432 real lBA = (coords[jN] - coords[iN]).norm();
1433 tPoint nBA = (coords[jN] - coords[iN]).normalized();
1434 real K = std::pow(lBA, -0.2);
1435 // real K = std::pow((wDist / elevationInfo.refDWall), -1);
1436 tJacobi mnBA = nBA * nBA.transpose() * K;
1437 for (auto &ME : A[iN])
1438 if (ME.j == jN)
1439 {
1440 A[iN][0].m += mnBA;
1441 ME.m -= mnBA;
1442 }
1443 }
1444 if (coordsElevDisp[iN](0) != largeReal)
1445 {
1446 nFix++;
1447 real nDiag = 0;
1448 for (auto &ME : A[iN])
1449 if (ME.j == iN)
1450 nDiag = ME.m.array().abs().maxCoeff();
1451 bO2[iN] = coordsElevDisp[iN] * nDiag;
1452 // bO2[iN].setZero();
1453 // if (coords[iN].norm() < 2)
1454 // bO2[iN](1) = sqr(coords[iN](0)) * 0.1 * nDiag;
1455 for (auto &ME : A[iN])
1456 if (ME.j == iN)
1457 ME.m = tJacobi::Identity() * nDiag;
1458 else
1459 ME.m.setZero();
1460 }
1461 if (coordsElevDisp[iN](2) == 2 * largeReal || coordsElevDisp[iN](2) == 3 * largeReal)
1462 {
1463 nFixB++;
1464 real nDiag = 0;
1465 for (auto &ME : A[iN])
1466 if (ME.j == iN)
1467 nDiag = ME.m.array().abs().maxCoeff();
1468 bO2[iN].setZero();
1469 // bO2[iN].setZero();
1470 // if (coords[iN].norm() < 2)
1471 // bO2[iN](1) = sqr(coords[iN](0)) * 0.1 * nDiag;
1472 for (auto &ME : A[iN])
1473 if (ME.j == iN)
1474 ME.m = tJacobi::Identity() * nDiag;
1475 else
1476 ME.m.setZero();
1477 }
1478 if (dim != 3)
1479 {
1480 real nDiag = 0;
1481 for (auto &ME : A[iN])
1482 if (ME.j == iN)
1483 nDiag = ME.m.array().abs().maxCoeff();
1484 for (auto &ME : A[iN])
1485 if (ME.j == iN)
1486 ME.m(2, 2) += nDiag; // the z-direction need a constraint
1487 }
1488 }
1489
1490 bO2.trans.startPersistentPull();
1491 bO2.trans.waitPersistentPull();
1492 // for (index iN = 0; iN < A.size(); iN++)
1493 // for (auto &ME : A[iN])
1494 // std::cout << "ij: " << iN << " " << ME.j << "\n"
1495 // << ME.m << std::endl;
1496 MPI::AllreduceOneIndex(nFix, MPI_SUM, mpi);
1497 MPI::AllreduceOneIndex(nFixB, MPI_SUM, mpi);
1498 if (mpi.rank == mRank)
1499 {
1500 log() << fmt::format("UnstructuredMesh === ElevatedNodesSolveInternalSmooth(): Matrix Assembled, nFix {}, nFixB {} ",
1501 nFix, nFixB)
1502 << std::endl;
1503 }
1504 };
1505 dispO2.setConstant(0.0);
1506 AssembleRHSMat(dispO2);
1507
1508 auto LimitDisp = [&](tCoordPair &uCur, tCoordPair &duCur, tCoordPair &duCurLim)
1509 {
1510 index nFix{0};
1511 real minLim{1};
1512 if (true) // bisect fem method
1513 for (index iCell = 0; iCell < cell2node.Size(); iCell++)
1514 {
1515 auto c2n = cell2node[iCell];
1516 rowsize nnLoc = c2n.size();
1517 SmallCoordsAsVector coordsC, coordsCu, coordsCuinc;
1518 GetCoordsOnCell(iCell, coordsC);
1519 GetCoordsOnCell(iCell, coordsCu, uCur);
1520 GetCoordsOnCell(iCell, coordsCuinc, duCur);
1521 auto localNodeIdx = Eigen::ArrayXi::LinSpaced(c2n.size(), 0, c2n.size() - 1);
1522
1523 auto eCell = GetCellElement(iCell);
1524 int nBi = eCell.GetO2NumBisect();
1525 int iBiVariant = GetO2ElemBisectVariant(eCell, coordsC);
1526
1527 auto FJDetRatio = [&](real alpha) -> real
1528 {
1529 real JDetRatio = 1;
1530 for (int iBi = 0; iBi < nBi; iBi++)
1531 {
1532 auto eCellSub = eCell.ObtainO2BisectElem(iBi);
1533 Eigen::ArrayXi c2nSubLocal;
1534 c2nSubLocal.resize(eCellSub.GetNumNodes());
1535 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, localNodeIdx, c2nSubLocal);
1536 SmallCoordsAsVector coordsCSub, coordsCuSub, coordsCuincSub;
1537 coordsCSub.resize(3, eCellSub.GetNumNodes());
1538 coordsCuSub.resize(3, eCellSub.GetNumNodes());
1539 coordsCuincSub.resize(3, eCellSub.GetNumNodes());
1540 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsC, coordsCSub);
1541 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsCu, coordsCuSub);
1542 eCell.ExtractO2BisectElemNodes(iBi, iBiVariant, coordsCuinc, coordsCuincSub);
1543 auto qCellSub = Elem::Quadrature{eCellSub, 2}; // for O1 FEM
1544 Elem::SummationNoOp NoOp;
1545 qCellSub.Integration(
1546 NoOp,
1547 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
1548 {
1549 // J = dx_i/dxii_j, Jinv = dxii_i/dx_j
1550 auto getJ = [&](auto coos)
1551 {
1552 tJacobi J = Elem::ShapeJacobianCoordD01Nj(coos, DiNj);
1553 tJacobi JInv = tJacobi::Identity();
1554 real JDet;
1555 if (dim == 2)
1556 JDet = J(EigenAll, 0).cross(J(EigenAll, 1))(2); // ! note: for 2D, incompatible with curve surface situation
1557 else
1558 JDet = J.determinant();
1559 if (dim == 2)
1560 JInv({0, 1}, {0, 1}) = J({0, 1}, {0, 1}).inverse().eval();
1561 else
1562 JInv = J.inverse().eval();
1563 return std::make_tuple(J, JInv, JDet);
1564 };
1565 auto [J, JInv, JDet] = getJ(coordsCSub + coordsCuSub);
1566 DNDS_assert(JDet > 0);
1567 auto [Jn, JInvn, JDetn] = getJ(coordsCSub + coordsCuSub + coordsCuincSub * alpha);
1568 JDetRatio = std::min(JDetRatio, JDetn / JDet);
1569 });
1570 }
1571 return JDetRatio;
1572 };
1573 real aLim = 1;
1574 real decay = 0.5;
1575 real minJDetRatio = 0.1;
1576 for (int i = 0; i < 20; i++)
1577 if (FJDetRatio(aLim) < minJDetRatio)
1578 aLim *= decay;
1579 for (auto iN : c2n)
1580 if (iN < duCurLim.father->Size())
1581 duCurLim[iN] = duCur[iN].array().sign() * duCurLim[iN].array().abs().min(duCur[iN].array().abs() * aLim);
1582 minLim = std::min(minLim, aLim);
1583 if (aLim < 1)
1584 nFix++;
1585 }
1586
1587 duCurLim.trans.startPersistentPull();
1588 duCurLim.trans.waitPersistentPull();
1589 // for (index iN = 0; iN < A.size(); iN++)
1590 // for (auto &ME : A[iN])
1591 // std::cout << "ij: " << iN << " " << ME.j << "\n"
1592 // << ME.m << std::endl;
1593 MPI::AllreduceOneIndex(nFix, MPI_SUM, mpi);
1594 MPI::AllreduceOneReal(minLim, MPI_MIN, mpi);
1595 if (mpi.rank == mRank)
1596 {
1597 log() << fmt::format("UnstructuredMesh === ElevatedNodesSolveInternalSmooth(): Disp Limited, nLim {}, minLim {:.3e} ",
1598 nFix, minLim)
1599 << std::endl;
1600 }
1601 return minLim;
1602 };
1603
1604 auto MatVec = [&](tCoordPair &x, tCoordPair &Ax)
1605 {
1606 for (index iN = 0; iN < x.father->Size(); iN++)
1607 {
1608 Ax[iN].setZero();
1609 for (auto &ME : A[iN])
1610 Ax[iN] += ME.m * x[ME.j];
1611 }
1612 Ax.trans.startPersistentPull();
1613 Ax.trans.waitPersistentPull();
1614 };
1615
1616 struct InitForEigenSparse
1617 {
1618 using value_type = index;
1619 std::vector<std::vector<MatElem>> &A;
1620 index operator[](index i) const
1621 {
1622 return A[i / 3].size() * 3;
1623 }
1624 } initForEigenSparse{A};
1625
1626 Eigen::SparseMatrix<real> ADiag(dispO2.father->Size() * 3, dispO2.father->Size() * 3);
1627 ADiag.reserve(initForEigenSparse);
1628 // std::cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" << std::endl;
1629 for (index iN = 0; iN < dispO2.father->Size(); iN++)
1630 {
1631 for (auto &ME : A[iN])
1632 if (ME.j < dispO2.father->Size())
1633 for (int iii = 0; iii < 3; iii++)
1634 for (int jjj = 0; jjj < 3; jjj++)
1635 ADiag.insert(iN * 3 + iii, ME.j * 3 + jjj) =
1636 ME.m(iii, jjj);
1637 }
1638 if (mpi.rank == mRank)
1639 {
1640 log() << "UnstructuredMesh === ElevatedNodesSolveInternalSmooth(): Eigen Matrix Filled" << std::endl;
1641 }
1642 // std::cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" << std::endl;
1643 Eigen::SparseLU<Eigen::SparseMatrix<real>, Eigen::COLAMDOrdering<int>> LUSolver;
1644 if (dispO2.father->Size())
1645 {
1646 // LUSolver.analyzePattern(ADiag);
1647 // LUSolver.factorize(ADiag); // full LU perconditioner
1648 // std::cout << "det " << LUSolver.absDeterminant() << std::endl;
1649 }
1650 if (mpi.rank == mRank)
1651 {
1652 log() << "UnstructuredMesh === ElevatedNodesSolveInternalSmooth(): LU Factorized" << std::endl;
1653 }
1654
1655 auto DiagSolveVec = [&](tCoordPair &x, tCoordPair &ADIx)
1656 {
1657 Eigen::Vector<real, Eigen::Dynamic> v;
1658 v.resize(x.father->Size() * 3);
1659 for (index iN = 0; iN < x.father->Size(); iN++)
1660 {
1661 v[iN * 3 + 0] = x[iN](0);
1662 v[iN * 3 + 1] = x[iN](1);
1663 v[iN * 3 + 2] = x[iN](2);
1664 }
1665 Eigen::Vector<real, Eigen::Dynamic> sol;
1666 if (dispO2.father->Size())
1667 {
1668 sol = LUSolver.solve(v);
1669 }
1670 for (index iN = 0; iN < x.father->Size(); iN++)
1671 {
1672 ADIx[iN](0) = sol[iN * 3 + 0];
1673 ADIx[iN](1) = sol[iN * 3 + 1];
1674 ADIx[iN](2) = sol[iN * 3 + 2];
1675 }
1676 ADIx.trans.startPersistentPull();
1677 ADIx.trans.waitPersistentPull();
1678 };
1679 auto JacobiPrec = [&](tCoordPair &x, tCoordPair &ADIx)
1680 {
1681 for (index iN = 0; iN < x.father->Size(); iN++)
1682 {
1683 ADIx[iN] = A[iN][0].m.fullPivLu().solve(x[iN]);
1684 }
1685 ADIx.trans.startPersistentPull();
1686 ADIx.trans.waitPersistentPull();
1687 };
1688 // for(int iIter = 1; iIter <= 200; iIter ++)
1689 // {
1690 // real incC = UpdateDispO2Jacobi(dispO2, dispO2, false); // JOR iteration
1691 // dispO2.trans.startPersistentPull();
1692 // dispO2.trans.waitPersistentPull();
1693 // if(mpi.rank == mRank)
1694 // log() << fmt::format("iIter [{}] Jacobi Increment [{:3e}]", iIter,incC) << std::endl;
1695 // }
1696
1697 /**********************************************/
1698 // Linear GMRES
1699
1700 Linear::GMRES_LeftPreconditioned<CoordPairDOF> gmres(5, initDOF);
1701 gmres.solve(
1702 [&](CoordPairDOF &x, CoordPairDOF &Ax)
1703 {
1704 MatVec(x, Ax);
1705 },
1706 [&](CoordPairDOF &x, CoordPairDOF &MLx)
1707 {
1708 JacobiPrec(x, MLx);
1709 },
1710 [&](CoordPairDOF &a, CoordPairDOF &b) -> real
1711 {
1712 return a.dot(b);
1713 },
1714 bO2,
1715 dispO2,
1716 elevationInfo.nIter,
1717 [&](int iRestart, real res, real resB)
1718 {
1719 // std::cout << mpi.rank << "gmres here" << std::endl;
1720 if (mpi.rank == mRank)
1721 log() << fmt::format("iRestart [{}] res [{:3e}] -> [{:3e}]", iRestart, resB, res) << std::endl;
1722 return res < resB * 1e-10;
1723 });
1724
1725 /**********************************************/
1726 // Nonlin
1727
1728 // Linear::GMRES_LeftPreconditioned<CoordPairDOF> gmres(5, initDOF);
1729 // real incNorm0{0}, rhsNorm0{0};
1730 // dispO2.setConstant(0.0);
1731 // for (int iIter = 1; iIter <= elevationInfo.nIter; iIter++)
1732 // {
1733 // bO2.setConstant(0.0); // remember this
1734 // AssembleRHSMat(dispO2);
1735 // MatVec(dispO2, dispO2PrecRHS); // "dispO2PrecRHS" = Ax
1736 // bO2.addTo(dispO2PrecRHS, -1.0); // RHS
1737 // real rhsNorm = bO2.norm2();
1738 // dispO2Inc.setConstant(0.0);
1739
1740 // gmres.solve(
1741 // [&](CoordPairDOF &x, CoordPairDOF &Ax)
1742 // {
1743 // MatVec(x, Ax);
1744 // },
1745 // [&](CoordPairDOF &x, CoordPairDOF &MLx)
1746 // {
1747 // JacobiPrec(x, MLx);
1748 // },
1749 // bO2,
1750 // dispO2Inc,
1751 // 2,
1752 // [&](int iRestart, real res, real resB)
1753 // {
1754 // // std::cout << mpi.rank << "gmres here" << std::endl;
1755 // if (mpi.rank == mRank)
1756 // log() << fmt::format("step [{}] iRestart [{}] res [{:3e}] -> [{:3e}]", iIter, iRestart, resB, res) << std::endl;
1757 // return res < resB * 1e-10;
1758 // });
1759 // dispO2IncLimited = dispO2Inc;
1760 // real minLim = LimitDisp(dispO2, dispO2Inc, dispO2IncLimited);
1761 // dispO2IncLimited = dispO2Inc;
1762 // dispO2IncLimited *= minLim;
1763 // dispO2.addTo(dispO2IncLimited, 1.0);
1764 // real incNorm = dispO2IncLimited.norm2();
1765 // if (iIter == 1)
1766 // incNorm0 = incNorm, rhsNorm0 = rhsNorm;
1767 // if (mpi.rank == mRank)
1768 // log() << fmt::format("step [{}] increment [{:3e}] -> [{:3e}] rhs [{:3e}] -> [{:3e}]", iIter, incNorm0, incNorm, rhsNorm0, rhsNorm) << std::endl;
1769 // }
1770 /**********************************************/
1771 // dispO2 = bO2;
1772
1773 for (index iN = 0; iN < coords.father->Size(); iN++)
1774 {
1775 if (dim == 2)
1776 dispO2[iN](2) = 0;
1777 // if(dispO2[iN].norm() != 0)
1778 // std::cout << dispO2[iN].transpose() << std::endl;
1779 coords[iN] += dispO2[iN];
1780 }
1781 coords.trans.pullOnce();
1782 }
1783}
ArrayEigenUniMatrixBatch.hpp
Batch of uniform-sized Eigen matrices per row, with variable batch count.
EigenPCH.hpp
Pre-compiled-header style shim that includes the heavy Eigen headers under DNDSR's warning suppressio...
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_MPI_InsertCheck
#define DNDS_MPI_InsertCheck(mpi, info)
Debug helper: barrier + print "CHECK <info> RANK r @ fn @ file:line".
Definition MPI.hpp:320
DNDS::Linear::GMRES_LeftPreconditioned
Definition Linear.hpp:11
DNDS::Linear::GMRES_LeftPreconditioned::solve
bool solve(TFA &&FA, TFML &&FML, TFDot &&fDot, TDATA &b, TDATA &x, uint32_t nRestart, TFstop &&FStop)
Definition Linear.hpp:45
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::RBF::RBFCPC2
MatrixXR RBFCPC2(Tcent cent, Tx xs, real R, RBFKernelType kernel=Gaussian)
Definition RadialBasisFunction.hpp:99
DNDS::Geom::RBF::FRBFBasis
MatrixXR FRBFBasis(TIn RiXj, RBFKernelType kernel)
Definition RadialBasisFunction.hpp:57
DNDS::Geom::RBF::RBFInterpolateSolveCoefsNoPoly
MatrixXR RBFInterpolateSolveCoefsNoPoly(const tSmallCoords &xs, const TF fs, real R, RBFKernelType kernel=Gaussian)
Definition RadialBasisFunction.hpp:136
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::Geom::tCoord
decltype(tCoordPair::father) tCoord
Definition Mesh_DeviceView.hpp:79
DNDS::MPI::AllreduceOneReal
void AllreduceOneReal(real &v, MPI_Op op, const MPIInfo &mpi)
Single-scalar Allreduce helper for reals (in-place, count = 1).
Definition MPI.hpp:673
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::AllreduceOneIndex
void AllreduceOneIndex(index &v, MPI_Op op, const MPIInfo &mpi)
Single-scalar Allreduce helper for indices (in-place, count = 1).
Definition MPI.hpp:679
DNDS::rowsize
int32_t rowsize
Row-width / per-row element-count type (signed 32-bit).
Definition Defines.hpp:109
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:205
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::log
std::ostream & log()
Return the current DNDSR log stream (either std::cout or the installed file).
Definition Defines.cpp:49
DNDS::DNDS_MPI_REAL
const MPI_Datatype DNDS_MPI_REAL
MPI datatype matching real (= MPI_REAL8).
Definition MPI.hpp:92
DNDS::ArrayPairDeviceView_Base::Size
DNDS_DEVICE_CALLABLE index Size() const
Combined father + son row count.
Definition ArrayPair.hpp:33
DNDS::ArrayPairDeviceView
Mutable device view onto an ArrayPair (for CUDA kernels).
Definition ArrayPair.hpp:105
DNDS::ArrayPair< DNDS::ArrayEigenVector< 3 > >::t_deviceView
ArrayPairDeviceView< B, DNDS::ArrayEigenVector< 3 > > t_deviceView
Device-view template alias: t_deviceView<DeviceBackend::CUDA> gives the mutable CUDA view type for th...
Definition ArrayPair.hpp:653
DNDS::ArrayPair< DNDS::ArrayEigenVector< 3 > >::father
ssp< DNDS::ArrayEigenVector< 3 > > father
Owned-side array (must be resized before ghost setup).
Definition ArrayPair.hpp:170
DNDS::ArrayPair< DNDS::ArrayEigenVector< 3 > >::Size
index Size() const
Combined row count (father->Size() + son->Size()).
Definition ArrayPair.hpp:292
DNDS::ArrayPair::trans
TTrans trans
Ghost-communication engine bound to father and son.
Definition ArrayPair.hpp:175
DNDS::Geom::CoordPairDOF::addTo
void addTo(CoordPairDOF &R, real alpha)
Definition Mesh_Elevation_SmoothSolver.cpp:39
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::PointCloudKDTreeCoordPair::kdtree_get_pt
real kdtree_get_pt(const size_t idx, const size_t dim) const
Definition Mesh_Elevation_SmoothSolver.cpp:81
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::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::GetCoordsOnCell
void GetCoordsOnCell(index iCell, tSmallCoords &cs)
Definition Mesh.hpp:566
DNDS::Geom::UnstructuredMesh::coords
tCoordPair coords
reader
Definition Mesh.hpp:54
DNDS::Geom::UnstructuredMesh::adjPrimaryState
MeshAdjState adjPrimaryState
Definition Mesh.hpp:41
DNDS::Geom::UnstructuredMesh::elevationInfo
struct DNDS::Geom::UnstructuredMesh::ElevationInfo elevationInfo
DNDS::Geom::UnstructuredMesh::adjFacialState
MeshAdjState adjFacialState
Definition Mesh.hpp:43
DNDS::Geom::UnstructuredMesh::elevState
MeshElevationState elevState
Definition Mesh.hpp:52
DNDS::MPIInfo::rank
int rank
This rank's 0-based index within comm (-1 until initialised).
Definition MPI.hpp:219
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
Ax
tVec Ax(NCells)
r
tVec r(NCells)
x
tVec x(NCells)
b
tVec b(NCells)
alpha
real alpha
Definition test_GasThermo.cpp:454
res
auto res
Definition test_ODE.cpp:196