DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
VariationalReconstruction.cpp
Go to the documentation of this file.
1
2#include "VariationalReconstruction.hpp"
3#include "DNDS/HardEigen.hpp"
4
5namespace DNDS::CFV
6{
7 template <int dim>
8 void
9 VariationalReconstruction<dim>::
10 ConstructMetrics()
11 {
12 using namespace Geom;
13 using namespace Geom::Elem;
14
15 /***************************************/
16 // get volumes
17 this->SetCellAtrBasic();
18 this->ConstructCellVolume();
19 this->ConstructCellBary();
20 this->ConstructCellCent();
21 this->ConstructCellIntJacobiDet();
22 this->ConstructCellIntPPhysics();
23 this->ConstructCellAlignedHBox();
24 this->ConstructCellMajorHBoxCoordInertia();
25
26 if (mpi.rank == mRank)
27 log()
28 << fmt::format(
29 "VariationalReconstruction<dim>::ConstructMetrics() === \n ===Sum/Min/Max Volume [{:.10g}; {:.5g}, {:.5g}]",
30 sumVolume, minVolume, maxVolume)
31 << std::endl;
32
33 /***************************************/
34 // get areas
35 this->SetFaceAtrBasic();
36 this->ConstructFaceCent();
37 this->ConstructFaceArea();
38 this->ConstructFaceIntJacobiDet();
39 this->ConstructFaceIntPPhysics();
40 this->ConstructFaceUnitNorm();
41 this->ConstructFaceMeanNorm();
42
43 this->ConstructCellSmoothScale();
44 }
45
46 template void
47 VariationalReconstruction<2>::
48 ConstructMetrics();
49 template void
50 VariationalReconstruction<3>::
51 ConstructMetrics();
52
53 template <int dim>
54 void
55 VariationalReconstruction<dim>::
56 ConstructBaseAndWeight(const tFGetBoundaryWeight &id2faceDircWeight)
57 {
58 using namespace Geom;
59 using namespace Geom::Elem;
60 using namespace Geom::Base;
61 int maxNDOF = GetNDof<dim>(settings.maxOrder);
62 // for polynomial: ndiff = ndof
63 int maxNDIFF = maxNDOF;
64 /******************************/
65 // *cell's moment and cache
66 this->MakePairDefaultOnCell(baseWeight_.cellBaseMoment, "VR::cellBaseMoment", maxNDOF);
67 if (settings.cacheDiffBase)
68 {
69 this->MakePairDefaultOnCell(
70 baseWeight_.cellDiffBaseCache, "VR::cellDiffBaseCache",
71 std::min(maxNDIFF, static_cast<int>(settings.cacheDiffBaseSize)), maxNDOF);
72 this->MakePairDefaultOnCell(
73 baseWeight_.cellDiffBaseCacheCent, "VR::cellDiffBaseCacheCent",
74 std::min(maxNDIFF, static_cast<int>(settings.cacheDiffBaseSize)), maxNDOF);
75 }
76 if (coefficients_.needVolIntCholeskyL)
77 coefficients_.volIntCholeskyL.resize(mesh->NumCellProc());
78#ifdef DNDS_USE_OMP
79# pragma omp parallel for
80#endif
81 for (index iCell = 0; iCell < mesh->NumCellProc(); iCell++)
82 {
83 GetCellAtr(iCell).relax = settings.jacobiRelax;
84 auto qCell = this->GetCellQuad(iCell);
85 auto qCellO1 = this->GetCellQuadO1(iCell);
86 auto qCellMax = Quadrature{mesh->GetCellElement(iCell), INT_ORDER_MAX};
87 if (settings.cacheDiffBase)
88 {
89 baseWeight_.cellDiffBaseCache.ResizeRow(iCell, qCell.GetNumPoints());
90 }
91 tSmallCoords coordsCell;
92 mesh->GetCoordsOnCell(iCell, coordsCell);
93
94 Eigen::RowVector<real, Eigen::Dynamic> m;
95 m.setZero(GetCellAtr(iCell).NDOF);
96 qCellMax.Integration(
97 m,
98 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
99 {
100 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coordsCell, DiNj);
101 real JDet = CellJacobianDet<dim>(coordsCell, DiNj) * (axisSymmetric ? pPhy(1) : 1.);
102 Eigen::RowVector<real, Eigen::Dynamic> vv;
103 vv.resizeLike(m);
104 this->FDiffBaseValue(vv, pPhy, iCell, -1, -2, 1); // un-dispatched call
105 DNDS_assert(vv(0) == 1); // must have 0th base
106 vInc = vv * JDet;
107 });
108 baseWeight_.cellBaseMoment[iCell] = m.transpose() / m(0); // 0th must be good
109 SummationNoOp noOp;
110 qCell.Integration(
111 noOp,
112 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
113 {
114 if (settings.cacheDiffBase)
115 {
116 MatrixXR dbv;
117 dbv.resize(
118 std::min(GetCellAtr(iCell).NDIFF, settings.cacheDiffBaseSize),
119 GetCellAtr(iCell).NDOF);
120 this->FDiffBaseValue(dbv, this->GetCellQuadraturePPhys(iCell, iG), iCell, -1, iG, 0);
121 baseWeight_.cellDiffBaseCache(iCell, iG) = dbv;
122 }
123 });
124 qCellO1.Integration(
125 noOp,
126 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
127 {
128 if (settings.cacheDiffBase)
129 {
130 MatrixXR dbv;
131 dbv.resize(
132 std::min(GetCellAtr(iCell).NDIFF, settings.cacheDiffBaseSize),
133 GetCellAtr(iCell).NDOF);
134 this->FDiffBaseValue(dbv, this->GetCellQuadraturePPhys(iCell, -1), iCell, -1, -1, 0);
135 baseWeight_.cellDiffBaseCacheCent[iCell] = dbv;
136 }
137 });
138
139 //****** Get Orthogonization Coefs
140 MatrixXR MBiBj;
141 MBiBj.setZero(GetCellAtr(iCell).NDOF - 1, GetCellAtr(iCell).NDOF - 1);
142 qCell.Integration(
143 MBiBj,
144 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
145 {
146 MatrixXR D0Bj =
147 this->GetIntPointDiffBaseValue(iCell, -1, -1, iG, std::array<int, 1>{0}, 1);
148 vInc = (D0Bj.transpose() * D0Bj);
149 vInc *= this->GetCellJacobiDet(iCell, iG);
150 });
151 if (coefficients_.needVolIntCholeskyL)
152 coefficients_.volIntCholeskyL.at(iCell).resizeLike(MBiBj),
153 coefficients_.volIntCholeskyL.at(iCell) = MBiBj.llt().matrixL();
154 }
155
156 /******************************/
157 // *face's weight and cache
158
159 this->MakePairDefaultOnFace(baseWeight_.faceWeight, "VR::faceWeight", settings.maxOrder + 1);
160 this->MakePairDefaultOnFace(baseWeight_.faceAlignedScales, "VR::faceAlignedScales");
161 this->MakePairDefaultOnFace(baseWeight_.faceMajorCoordScale, "VR::faceMajorCoordScale", 3, 3);
162 if (settings.cacheDiffBase)
163 {
164 this->MakePairDefaultOnFace(
165 baseWeight_.faceDiffBaseCache, "VR::faceDiffBaseCache",
166 std::min(maxNDIFF, static_cast<int>(settings.cacheDiffBaseSize)), maxNDOF);
167 this->MakePairDefaultOnFace(
168 baseWeight_.faceDiffBaseCacheCent, "VR::faceDiffBaseCacheCent",
169 std::min(maxNDIFF, static_cast<int>(settings.cacheDiffBaseSize)), maxNDOF * 2);
170 }
171#ifdef DNDS_USE_OMP
172# pragma omp parallel for
173#endif
174 for (index iFace = 0; iFace < mesh->NumFaceProc(); iFace++)
175 {
176 auto qFace = this->GetFaceQuad(iFace);
177 auto qFaceO1 = this->GetFaceQuadO1(iFace);
178 if (settings.cacheDiffBase)
179 baseWeight_.faceDiffBaseCache.ResizeRow(iFace, 2 * qFace.GetNumPoints());
180
181 tSmallCoords coords, coordsC;
182 mesh->GetCoordsOnFace(iFace, coords);
183 coordsC = coords;
184 coordsC.colwise() -= faceCent[iFace];
185
186 // * get bdv cache
187 SummationNoOp noOp;
188 for (int if2c = 0; if2c < 2; if2c++)
189 {
190 index iCell = mesh->face2cell(iFace, if2c);
191 if (iCell == UnInitIndex)
192 continue;
193 if (FaceIDIsExternalBC(mesh->GetFaceZone(iFace)))
194 DNDS_assert(if2c == 0);
195 else if (FaceIDIsPeriodic(mesh->GetFaceZone(iFace)))
196 {
197 // TODO: handle the case with periodic
198 // DNDS_assert(false); //! do nothing?
199 }
200 qFace.Integration(
201 noOp,
202 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
203 {
204 if (settings.cacheDiffBase)
205 {
206 MatrixXR dbv;
207 dbv.resize(
208 std::min(GetFaceAtr(iFace).NDIFF, settings.cacheDiffBaseSize),
209 GetCellAtr(iCell).NDOF);
210 this->FDiffBaseValue(dbv, this->GetFaceQuadraturePPhysFromCell(iFace, iCell, if2c, iG), iCell, iFace, iG, 0);
211 baseWeight_.faceDiffBaseCache(iFace, if2c * qFace.GetNumPoints() + iG) = dbv;
212 }
213 });
214 qFaceO1.Integration(
215 noOp,
216 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
217 {
218 if (settings.cacheDiffBase)
219 {
220 MatrixXR dbv;
221 dbv.resize(
222 std::min(GetFaceAtr(iFace).NDIFF, settings.cacheDiffBaseSize),
223 GetCellAtr(iCell).NDOF);
224 this->FDiffBaseValue(dbv, this->GetFaceQuadraturePPhysFromCell(iFace, iCell, if2c, -1), iCell, iFace, -1, 0);
225 int maxNDOF = GetNDof<dim>(settings.maxOrder);
226 baseWeight_.faceDiffBaseCacheCent[iFace](
227 EigenAll,
228 Eigen::seq(if2c * maxNDOF,
229 if2c * maxNDOF + maxNDOF - 1)) = dbv;
230 }
231 });
232 }
233
234 // *get face (derivatives) scale: cell average AlignedHBox mode
235 //! note: if use cell-2-cell distance, note periodic, use wrapped call here
236 tPoint faceScale{0, 0, 0};
237 int nF2C{0};
238 for (int if2c = 0; if2c < 2; if2c++)
239 {
240 index iCell = mesh->face2cell(iFace, if2c);
241 if (iCell == UnInitIndex)
242 continue;
243 if (FaceIDIsExternalBC(mesh->GetFaceZone(iFace)))
244 DNDS_assert(if2c == 0);
245 else if (FaceIDIsPeriodic(mesh->GetFaceZone(iFace)))
246 {
247 // TODO: handle the case with periodic
248 // DNDS_assert(false); // !do nothing for now
249 }
250 faceScale += cellAlignedHBox[iCell];
251 nF2C++;
252 }
253 DNDS_assert(nF2C > 0);
254 faceScale /= nF2C;
255 baseWeight_.faceAlignedScales[iFace] = faceScale;
256 Geom::tPoint faceBaryDiffV;
257 if (mesh->face2cell(iFace, 1) != UnInitIndex)
258 {
259 faceBaryDiffV =
260 this->GetOtherCellBaryFromCell(mesh->face2cell(iFace, 0),
261 mesh->face2cell(iFace, 1), iFace) -
262 this->GetCellBary(mesh->face2cell(iFace, 0));
263 }
264 else
265 {
266 Geom::tPoint vCB = this->GetFaceQuadraturePPhysFromCell(
267 iFace,
268 mesh->face2cell(iFace, 0), 0, -1) -
269 this->GetCellBary(mesh->face2cell(iFace, 0));
270 auto uNorm = this->GetFaceNormFromCell(
271 iFace,
272 mesh->face2cell(iFace, 0), 0, -1);
273 faceBaryDiffV =
274 2 * vCB.dot(uNorm) *
275 this->GetFaceNorm(iFace, -1);
276 }
277 Geom::tPoint faceBaryDiffL, faceBaryDiffR;
278 faceBaryDiffL =
279 this->GetCellBary(mesh->face2cell(iFace, 0)) -
280 this->GetFaceQuadraturePPhysFromCell(iFace, mesh->face2cell(iFace, 0), 0, -1);
281 faceBaryDiffR = faceBaryDiffV + faceBaryDiffL;
282 real volL, volR;
283 volL = volR = std::pow(GetCellVol(mesh->face2cell(iFace, 0)) + verySmallReal, settings.functionalSettings.inertiaWeightPower);
284 Geom::tGPoint cellInertiaL = cellInertia[mesh->face2cell(iFace, 0)] * volL;
285 Geom::tGPoint cellInertiaR = cellInertiaL;
286 if (mesh->face2cell(iFace, 1) != UnInitIndex)
287 {
288 volR = std::pow(GetCellVol(mesh->face2cell(iFace, 1)) + verySmallReal, settings.functionalSettings.inertiaWeightPower);
289 cellInertiaR = this->GetOtherCellInertiaFromCell(
290 mesh->face2cell(iFace, 0),
291 mesh->face2cell(iFace, 1), iFace) *
292 volR;
293 }
294 Geom::tGPoint faceInertiaC =
295 (cellInertiaL + cellInertiaR) / (volL + volR);
296
297 Geom::tGPoint faceCoord;
298 faceCoord.setIdentity();
299 if constexpr (dim == 3)
300 faceCoord = HardEigen::Eigen3x3RealSymEigenDecomposition(faceInertiaC);
301 else
302 faceCoord({0, 1}, {0, 1}) = HardEigen::Eigen2x2RealSymEigenDecomposition(faceInertiaC({0, 1}, {0, 1}));
303 baseWeight_.faceMajorCoordScale[iFace] = faceCoord;
304 if (settings.functionalSettings.anisotropicType == VRSettings::FunctionalSettings::AnisotropicType::InertiaCoordBB ||
305 settings.functionalSettings.anisotropicType == VRSettings::FunctionalSettings::AnisotropicType::InertiaCoordBBNorm ||
306 settings.functionalSettings.anisotropicType == VRSettings::FunctionalSettings::AnisotropicType::InertiaCoordBBSym)
307 {
308 Geom::tGPoint faceCoordNorm = faceCoord.colwise().normalized();
309 tSmallCoords coordsL, coordsR;
310 coordsL = coordsC.colwise() - faceBaryDiffL;
311 coordsR = coordsC.colwise() - faceBaryDiffR;
312 coordsL = faceCoordNorm.transpose() * coordsL;
313 coordsR = faceCoordNorm.transpose() * coordsR;
314 Geom::tPoint faceCoordBB =
315 coordsL.array().abs().rowwise().maxCoeff().max(
316 coordsR.array().abs().rowwise().maxCoeff());
317
318 baseWeight_.faceMajorCoordScale[iFace] = faceCoordNorm * faceCoordBB.asDiagonal();
319 }
320 real AR2 = baseWeight_.faceMajorCoordScale[iFace].col(0).norm() / baseWeight_.faceMajorCoordScale[iFace].col(1).norm();
321 if (settings.functionalSettings.scaleType == VRSettings::FunctionalSettings::ScaleType::BaryDiff)
322 {
323 baseWeight_.faceAlignedScales[iFace] = faceBaryDiffV;
324 if constexpr (dim == 2)
325 baseWeight_.faceAlignedScales[iFace](2) = 1;
326 }
327
328 // *get geom weight ic2f
329 real wg = 1;
330 if (settings.functionalSettings.geomWeightScheme == VRSettings::FunctionalSettings::GeomWeightScheme::HQM_SD)
331 wg = settings.functionalSettings.geomWeightBias +
332 std::pow(std::pow(this->GetFaceArea(iFace), 1. / real(dim - 1)) / faceBaryDiffV.norm(), settings.functionalSettings.geomWeightPower * 0.5);
333 if (settings.functionalSettings.geomWeightScheme == VRSettings::FunctionalSettings::GeomWeightScheme::SD_Power)
334 wg = settings.functionalSettings.geomWeightBias +
335 std::pow(this->GetFaceArea(iFace), settings.functionalSettings.geomWeightPower1 * 0.5) *
336 std::pow(faceBaryDiffV.norm(), settings.functionalSettings.geomWeightPower2 * 0.5);
337
338 // *get dir weight
339 Eigen::Vector<real, Eigen::Dynamic> wd;
340 wd.resize(settings.maxOrder + 1);
341 switch (settings.functionalSettings.dirWeightScheme)
342 {
343 case VRSettings::FunctionalSettings::DirWeightScheme::Factorial:
344 for (int p = 0; p < wd.size(); p++)
345 wd[p] = 1. / factorials[p];
346 break;
347 case VRSettings::FunctionalSettings::DirWeightScheme::HQM_OPT:
348 switch (settings.maxOrder)
349 {
350 case 1:
351 wd[0] = wd[1] = 1;
352 break;
353 case 2:
354 wd[0] = 1, wd[1] = 0.4643, wd[2] = 0.1559;
355 break;
356 case 3:
357 wd[0] = 1, wd[1] = .5295, wd[2] = wd[3] = .2117;
358 break;
359 default:
360 DNDS_assert(false);
361 break;
362 }
363 break;
364 case VRSettings::FunctionalSettings::DirWeightScheme::ManualDirWeight:
365 switch (settings.maxOrder)
366 {
367 case 3:
368 wd[3] = settings.functionalSettings.manualDirWeights[3];
369 [[fallthrough]];
370 case 2:
371 wd[2] = settings.functionalSettings.manualDirWeights[2];
372 [[fallthrough]];
373 case 1:
374 wd[1] = settings.functionalSettings.manualDirWeights[1];
375 wd[0] = settings.functionalSettings.manualDirWeights[0];
376 break;
377 default:
378 DNDS_assert(false);
379 break;
380 }
381 break;
382 case VRSettings::FunctionalSettings::DirWeightScheme::TEST_OPT:
383 {
384 if (settings.maxOrder != 3) // use manual value
385 {
386 switch (settings.maxOrder)
387 {
388 case 3:
389 wd[3] = settings.functionalSettings.manualDirWeights[3];
390 [[fallthrough]];
391 case 2:
392 wd[2] = settings.functionalSettings.manualDirWeights[2];
393 [[fallthrough]];
394 case 1:
395 wd[1] = settings.functionalSettings.manualDirWeights[1];
396 wd[0] = settings.functionalSettings.manualDirWeights[0];
397 break;
398 default:
399 DNDS_assert(false);
400 break;
401 }
402 }
403 else if (dim == 3) // current scheme of a/d b/d geom weighting
404 {
405 auto faceSpan = Geom::Elem::GetElemNodeMajorSpan(coords);
406 real a = faceSpan(0);
407 real b = faceSpan(1);
408 real d = faceBaryDiffV.norm();
409 real r1 = a / (d + verySmallReal);
410 real r2 = b / (d + verySmallReal);
411 real rr = std::sqrt(sqr(r1) + sqr(r2));
412 real lrr = std::log10(rr);
413
414 wd[0] = std::sqrt(std::pow(1.0031 + 0.86 * std::pow(r1 * r2, 0.5327), 0.98132));
415 wd[1] = wd[0] * std::sqrt(0.3604 - 0.0774 * std::exp(-0.5 * sqr((lrr - 1.8e-10) / 0.1325)));
416 wd[2] = wd[0] * std::sqrt(0.0676 - 0.0547 * std::exp(-0.5 * sqr((lrr + 2.74e-10) / 0.1376)));
417 wd[3] = wd[0] * std::sqrt(5.299e-2 - 4.63e-2 * std::exp(-0.5 * sqr((lrr - 9.028e-12) / 0.175)));
418 wg = 1; // force no wg
419 }
420 else if (dim == 2)
421 {
422 real r = this->GetFaceArea(iFace) / faceBaryDiffV.norm();
423 real lr = std::log10(r);
424 wd[0] = std::sqrt(std::pow(1.00745 + r, 1.004));
425 wd[1] = wd[0] * std::sqrt(0.359 - 0.0674 * std::exp(-sqr(lr - 5.18e-11) / (2 * sqr(0.124))));
426 wd[2] = wd[0] * std::sqrt(0.0676 - 0.0507 * std::exp(-sqr(lr - 3.15e-10) / (2 * sqr(0.137))));
427 wd[3] = wd[0] * std::sqrt(0.0529 - 0.0487 * std::exp(-sqr(lr - 1.33e-10) / (2 * sqr(0.167))));
428 wg = 1; // force no wg
429 }
430 }
431 break;
432 default:
433 DNDS_assert(false);
434 break;
435 }
436
437 if (FaceIDIsExternalBC(mesh->GetFaceZone(iFace)))
438 {
439 for (int iOrder = 0; iOrder <= settings.maxOrder; iOrder++)
440 wd(iOrder) *= id2faceDircWeight(mesh->GetFaceZone(iFace), iOrder) * settings.bcWeight;
441 }
442 baseWeight_.faceWeight[iFace] = wd * wg;
443 }
444
445 if (settings.cacheDiffBase)
446 {
447 baseWeight_.faceDiffBaseCache.CompressBoth();
448 baseWeight_.cellDiffBaseCache.CompressBoth();
449 }
450
451 if (!settings.intOrderVRBCIsSame())
452 {
453 for (index iCell = 0; iCell < mesh->NumCell(); iCell++)
454 {
455 auto c2f = mesh->cell2face[iCell];
456 for (int ic2f = 0; ic2f < c2f.size(); ic2f++)
457 {
458 index iFace = c2f[ic2f];
459 index iCellOther = CellFaceOther(iCell, iFace);
460 auto faceID = mesh->GetFaceZone(iFace);
461 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
462 if (iCellOther == UnInitIndex)
463 {
464 DNDS_assert(FaceIDIsExternalBC(faceID));
465 DNDS_assert(baseWeight_.bndVRCaches.count(iFace) == 0);
466 auto qFace = Quadrature(mesh->GetFaceElement(iFace), settings.intOrderVRBCValue());
467 baseWeight_.bndVRCaches[iFace].reserve(qFace.GetNumPoints());
468 tSmallCoords coords;
469 mesh->GetCoordsOnFace(iFace, coords);
470 SummationNoOp noOp;
471 qFace.Integration(
472 noOp,
473 [&](auto &vInc, int __xxx_iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
474 {
475 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coords, DiNj);
476 tJacobi J = Elem::ShapeJacobianCoordD01Nj(coords, DiNj);
477 real JDet = JacobiDetFace<dim>(J) * (axisSymmetric ? pPhy(1) : 1.);
478 tPoint np = FacialJacobianToNormVec<dim>(J);
479 RowVectorXR dbv, dbvD;
480 dbvD.resize(1, GetCellAtr(iCell).NDOF);
481 this->FDiffBaseValue(dbvD, this->GetFacePointFromCell(iFace, iCell, -1, pPhy), iCell, iFace, -2, 0);
482 dbv = dbvD(0, Eigen::seq(Eigen::fix<1>, EigenLast));
483 BndVRPointCache cacheEntry;
484 cacheEntry.D0Bj = dbv;
485 cacheEntry.JDet = JDet;
486 cacheEntry.norm = np;
487 cacheEntry.PPhy = pPhy;
488 baseWeight_.bndVRCaches.at(iFace).emplace_back(std::move(cacheEntry));
489 });
490 }
491 }
492 }
493 }
494
495 if (mpi.rank == mRank)
496 log()
497 << "VariationalReconstruction<dim>::ConstructBaseAndWeight() done" << std::endl;
498 }
499
500 template void
501 VariationalReconstruction<2>::
502 ConstructBaseAndWeight(const tFGetBoundaryWeight &id2faceDircWeight);
503 template void
504 VariationalReconstruction<3>::
505 ConstructBaseAndWeight(const tFGetBoundaryWeight &id2faceDircWeight);
506
507 template <int dim>
508 void
509 VariationalReconstruction<dim>::ConstructRecCoeff()
510 {
511 auto &c_ = coefficients_;
512 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
513 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
514 using namespace Geom;
515 using namespace Geom::Elem;
516 using namespace Geom::Base;
517 int maxNDOF = GetNDof<dim>(settings.maxOrder);
518 if (c_.needOriginalMatrix)
519 this->MakePairDefaultOnCell(c_.matrixAB, "VR::matrixAB", maxNDOF - 1, maxNDOF - 1);
520 this->MakePairDefaultOnCell(c_.matrixAAInvB, "VR::matrixAAInvB", maxNDOF - 1, maxNDOF - 1);
521 if (c_.needOriginalMatrix)
522 this->MakePairDefaultOnCell(c_.vectorB, "VR::vectorB", maxNDOF - 1, 1);
523 this->MakePairDefaultOnCell(c_.vectorAInvB, "VR::vectorAInvB", maxNDOF - 1, 1);
524 if (settings.functionalSettings.greenGauss1Weight != 0)
525 this->MakePairDefaultOnCell(c_.matrixAHalf_GG, "VR::matrixAHalf_GG", dim, (maxNDOF - 1));
526 this->MakePairDefaultOnCell(c_.matrixA, "VR::matrixA", (maxNDOF - 1), (maxNDOF - 1));
527 if (c_.needMatrixACholeskyL)
528 c_.matrixACholeskyL.resize(mesh->NumCellProc());
529 real maxCond = 0.0;
530#ifdef DNDS_USE_OMP
531# pragma omp parallel for
532#endif
533 for (index iCell = 0; iCell < mesh->NumCell(); iCell++) // only non-ghost
534 {
535 if (c_.needOriginalMatrix)
536 c_.matrixAB.ResizeRow(iCell, mesh->cell2face.RowSize(iCell) + 1);
537 c_.matrixAAInvB.ResizeRow(iCell, mesh->cell2face.RowSize(iCell) + 1);
538 if (c_.needOriginalMatrix)
539 c_.vectorB.ResizeRow(iCell, mesh->cell2face.RowSize(iCell));
540 c_.vectorAInvB.ResizeRow(iCell, mesh->cell2face.RowSize(iCell));
541 std::vector<MatrixXR> local_Bs;
542 std::vector<Eigen::Vector<real, Eigen::Dynamic>> local_bs;
543 local_Bs.reserve(mesh->cell2face.RowSize(iCell));
544 local_bs.reserve(mesh->cell2face.RowSize(iCell));
545 //*get A
546 MatrixXR A;
547 A.setZero(GetCellAtr(iCell).NDOF - 1, GetCellAtr(iCell).NDOF - 1);
548 for (int ic2f = 0; ic2f < mesh->cell2face.RowSize(iCell); ic2f++)
549 {
550 index iFace = mesh->cell2face(iCell, ic2f);
551 index iCellOther = CellFaceOther(iCell, iFace);
552 auto qFace = this->GetFaceQuad(iFace);
553 auto qFaceVR = Quadrature(mesh->GetFaceElement(iFace), iCellOther == UnInitIndex ? settings.intOrderVRBCValue()
554 : settings.intOrderVRValue());
555 bool cur_face_intOrderVRIsSame = qFace.int_order == qFaceVR.int_order;
556 tSmallCoords coords;
557 if (!cur_face_intOrderVRIsSame)
558 mesh->GetCoordsOnFace(iFace, coords);
559 qFaceVR.Integration(
560 A,
561 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
562 {
563 decltype(A) DiffI;
564 real JDet{0};
565 if (cur_face_intOrderVRIsSame)
566 {
567 DiffI = this->GetIntPointDiffBaseValue(iCell, iFace, -1, iG, EigenAll);
568 JDet = this->GetFaceJacobiDet(iFace, iG);
569 }
570 else
571 {
572 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coords, DiNj);
573 JDet = FaceJacobianDet(dim, coords, DiNj) * (axisSymmetric ? pPhy(1) : 1.);
574 MatrixXR dbv;
575 dbv.resize(GetFaceAtr(iFace).NDIFF, GetCellAtr(iCell).NDOF);
576 this->FDiffBaseValue(dbv, this->GetFacePointFromCell(iFace, iCell, -1, pPhy), iCell, iFace, -2, 0);
577 DiffI = dbv(EigenAll, Eigen::seq(Eigen::fix<1>, EigenLast));
578 }
579 vInc = this->FFaceFunctional(DiffI, DiffI, iFace, iCell, iCell);
580 vInc *= JDet;
581 });
582 }
583 Eigen::Matrix<real, dim, Eigen::Dynamic> AHalf_GG;
584 if (settings.functionalSettings.greenGauss1Weight != 0)
585 {
586 //* reduced functional on compact stencil
587 AHalf_GG.resize(Eigen::NoChange, A.cols());
588 AHalf_GG.setZero();
589 // AHalf's central part:
590 auto qCell = this->GetCellQuad(iCell);
591 qCell.IntegrationSimple(
592 AHalf_GG,
593 [&](decltype(AHalf_GG) &inc, int iG)
594 {
595 inc = this->GetIntPointDiffBaseValue(
596 iCell, -1, -1, iG,
597 Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), dim + 1);
598 // using no scaling here
599 inc *= this->GetCellJacobiDet(iCell, iG);
600 });
601 for (int ic2f = 0; ic2f < mesh->cell2face.RowSize(iCell); ic2f++)
602 {
603 index iFace = mesh->cell2face(iCell, ic2f);
604 index iCellOther = this->CellFaceOther(iCell, iFace);
605 auto qFace = this->GetFaceQuad(iFace);
606 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
607 qFace.IntegrationSimple(
608 AHalf_GG,
609 [&](decltype(AHalf_GG) &inc, int iG)
610 {
611 tPoint normOut = this->GetFaceNormFromCell(iFace, iCell, if2c, iG) *
612 (if2c ? -1 : 1);
613 inc = normOut(Seq012) *
614 this->GetIntPointDiffBaseValue(
615 iCell, iFace, -1, iG,
616 0, 1);
617 inc *= (1 - settings.functionalSettings.greenGauss1Bias);
618 if (iCellOther != UnInitIndex)
619 {
620 real dLR = (GetOtherCellBaryFromCell(iCell, iCellOther, iFace) - GetCellBary(iCell)).norm();
621 inc -= normOut(Seq012) *
622 (normOut(Seq012).transpose() *
623 this->GetIntPointDiffBaseValue(iCell, iFace, -1, iG, Seq123, dim + 1)) *
624 dLR *
625 settings.functionalSettings.greenGauss1Penalty;
626 }
627
628 inc *= (-1) * this->GetFaceJacobiDet(iFace, iG);
629 });
630 }
631 A += (AHalf_GG.transpose() * AHalf_GG) * this->GetGreenGauss1WeightOnCell(iCell);
632 c_.matrixAHalf_GG[iCell] = AHalf_GG;
633 }
634
635 //*get B
636 for (int ic2f = 0; ic2f < mesh->cell2face.RowSize(iCell); ic2f++)
637 {
638 index iFace = mesh->cell2face(iCell, ic2f);
639 index iCellOther = this->CellFaceOther(iCell, iFace);
640 auto qFace = this->GetFaceQuad(iFace);
641 auto qFaceVR = Quadrature(mesh->GetFaceElement(iFace),
642 iCellOther == UnInitIndex ? settings.intOrderVRBCValue()
643 : settings.intOrderVRValue());
644 bool cur_face_intOrderVRIsSame = qFace.int_order == qFaceVR.int_order;
645 tSmallCoords coords;
646 if (!cur_face_intOrderVRIsSame)
647 mesh->GetCoordsOnFace(iFace, coords);
648 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
649 auto faceID = mesh->GetFaceZone(iFace);
650 if (FaceIDIsExternalBC(mesh->GetFaceZone(iFace)))
651 {
652 DNDS_assert(iCellOther == UnInitIndex);
653 local_Bs.emplace_back();
654 continue;
655 // if is periodic, then use internal //TODO!
656 }
657 MatrixXR B;
658 B.setZero(GetCellAtr(iCell).NDOF - 1, GetCellAtr(iCellOther).NDOF - 1);
659 qFaceVR.Integration(
660 B,
661 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
662 {
663 decltype(B) DiffI;
664 decltype(B) DiffJ;
665 real JDet{0};
666 if (cur_face_intOrderVRIsSame)
667 {
668 JDet = this->GetFaceJacobiDet(iFace, iG);
669 DiffI = this->GetIntPointDiffBaseValue(iCell, iFace, -1, iG, EigenAll);
670 DiffJ = this->GetIntPointDiffBaseValue(iCellOther, iFace, -1, iG, EigenAll);
671 }
672 else
673 {
674 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coords, DiNj);
675 JDet = FaceJacobianDet(dim, coords, DiNj) * (axisSymmetric ? pPhy(1) : 1.);
676 MatrixXR dbvI, dbvJ;
677 dbvI.resize(GetFaceAtr(iFace).NDIFF, GetCellAtr(iCell).NDOF);
678 dbvJ.resize(GetFaceAtr(iFace).NDIFF, GetCellAtr(iCellOther).NDOF);
679 this->FDiffBaseValue(dbvI, this->GetFacePointFromCell(iFace, iCell, -1, pPhy), iCell, iFace, -2, 0);
680 this->FDiffBaseValue(dbvJ, this->GetFacePointFromCell(iFace, iCellOther, -1, pPhy), iCellOther, iFace, -2, 0);
681
682 DiffI = dbvI(EigenAll, Eigen::seq(Eigen::fix<1>, EigenLast));
683 DiffJ = dbvJ(EigenAll, Eigen::seq(Eigen::fix<1>, EigenLast));
684 }
685 if (mesh->isPeriodic)
686 {
687 if ((if2c == 1 && Geom::FaceIDIsPeriodicMain(faceID)) ||
688 (if2c == 0 && Geom::FaceIDIsPeriodicDonor(faceID))) // I am donor
689 {
690 periodicity.TransDiValueInplace<dim>(DiffJ, faceID);
691 }
692 if ((if2c == 1 && Geom::FaceIDIsPeriodicDonor(faceID)) ||
693 (if2c == 0 && Geom::FaceIDIsPeriodicMain(faceID))) // I am main
694 {
695 periodicity.TransDiValueBackInplace<dim>(DiffJ, faceID);
696 }
697 }
698
699 vInc = this->FFaceFunctional(DiffI, DiffJ, iFace, iCell, iCellOther);
700 vInc *= JDet;
701 });
702 if (settings.functionalSettings.greenGauss1Weight != 0)
703 {
704 decltype(AHalf_GG) BHalf_GG;
705 BHalf_GG.resize(Eigen::NoChange, B.cols());
706 BHalf_GG.setZero();
707 qFace.IntegrationSimple(
708 BHalf_GG,
709 [&](decltype(BHalf_GG) &inc, int iG)
710 {
711 tPoint normOut = this->GetFaceNormFromCell(iFace, iCell, if2c, iG) *
712 (if2c ? -1 : 1);
713 inc = normOut(Seq012) *
714 this->GetIntPointDiffBaseValue(iCellOther, iFace, 1 - if2c, iG, 0, 1); // need 1-if2c!!if2c is for iCell!
715 inc *= settings.functionalSettings.greenGauss1Bias;
716 real dLR = (GetOtherCellBaryFromCell(iCell, iCellOther, iFace) - GetCellBary(iCell)).norm();
717 inc += normOut(Seq012) *
718 (normOut(Seq012).transpose() *
719 this->GetIntPointDiffBaseValue(iCellOther, iFace, 1 - if2c, iG, Seq123, dim + 1)) *
720 dLR *
721 settings.functionalSettings.greenGauss1Penalty;
722 inc *= this->GetFaceJacobiDet(iFace, iG);
723 });
724 B += AHalf_GG.transpose() * BHalf_GG * this->GetGreenGauss1WeightOnCell(iCell);
725 }
726 if (iCellOther == iCell) //* coincide periodic
727 A -= B, B *= 0;
728 if (c_.needOriginalMatrix)
729 c_.matrixAB(iCell, 1 + ic2f) = B;
730 local_Bs.emplace_back(std::move(B));
731 }
732
733 //*get b
734 for (int ic2f = 0; ic2f < mesh->cell2face.RowSize(iCell); ic2f++)
735 {
736 index iFace = mesh->cell2face(iCell, ic2f);
737 index iCellOther = this->CellFaceOther(iCell, iFace);
738 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
739 auto qFace = this->GetFaceQuad(iFace);
740 auto qFaceVR = Quadrature(mesh->GetFaceElement(iFace),
741 iCellOther == UnInitIndex ? settings.intOrderVRBCValue()
742 : settings.intOrderVRValue());
743 bool cur_face_intOrderVRIsSame = qFace.int_order == qFaceVR.int_order;
744 tSmallCoords coords;
745 if (!cur_face_intOrderVRIsSame)
746 mesh->GetCoordsOnFace(iFace, coords);
747 VectorXR b;
748 b.setZero(GetCellAtr(iCell).NDOF - 1, 1);
749 qFaceVR.Integration(
750 b,
751 [&](auto &vInc, int iG, const tPoint &pParam, const Elem::tD01Nj &DiNj)
752 {
753 Eigen::RowVector<real, Eigen::Dynamic> DiffI;
754 real JDet{0};
755 if (cur_face_intOrderVRIsSame)
756 {
757 JDet = this->GetFaceJacobiDet(iFace, iG);
758 DiffI = this->GetIntPointDiffBaseValue(iCell, iFace, -1, iG, std::array<int, 1>{0}, 1);
759 }
760 else
761 {
762 tPoint pPhy = Elem::PPhysicsCoordD01Nj(coords, DiNj);
763 JDet = FaceJacobianDet(dim, coords, DiNj) * (axisSymmetric ? pPhy(1) : 1.);
764 MatrixXR dbv;
765 dbv.resize(1, GetCellAtr(iCell).NDOF);
766 this->FDiffBaseValue(dbv, this->GetFacePointFromCell(iFace, iCell, -1, pPhy), iCell, iFace, -2, 0);
767 DiffI = dbv(EigenAll, Eigen::seq(Eigen::fix<1>, EigenLast));
768 }
769 vInc = this->FFaceFunctional(DiffI, Eigen::MatrixXd::Ones(1, 1), iFace, iCell, iCell);
770 vInc *= JDet;
771 });
772 if (settings.functionalSettings.greenGauss1Weight != 0)
773 {
774 Eigen::Matrix<real, dim, 1> bHalf_GG;
775 bHalf_GG.setZero();
776 qFace.IntegrationSimple(
777 bHalf_GG,
778 [&](auto &inc, int iG)
779 {
780 inc = this->GetFaceNormFromCell(iFace, iCell, -1, iG)(Seq012) *
781 (if2c ? -1 : 1);
782 inc *= settings.functionalSettings.greenGauss1Bias * this->GetFaceJacobiDet(iFace, iG);
783 });
784 b += AHalf_GG.transpose() * bHalf_GG * this->GetGreenGauss1WeightOnCell(iCell);
785 }
786 if (c_.needOriginalMatrix)
787 c_.vectorB(iCell, ic2f) = b;
788 local_bs.emplace_back(std::move(b));
789 }
790
791 // * get AInv and fill A, AInv
792
793 decltype(A) AInv;
794 real aCond{0};
795 if (settings.svdTolerance)
796 aCond = HardEigen::EigenLeastSquareInverse_Filtered(A, AInv, settings.svdTolerance, 1);
797 else
798 aCond = HardEigen::EigenLeastSquareInverse(A, AInv, settings.svdTolerance);
799 if (c_.needOriginalMatrix)
800 c_.matrixAB(iCell, 0) = A;
801 c_.matrixAAInvB(iCell, 0) = AInv;
802 c_.matrixA[iCell] = A;
803
804 maxCond = std::max(aCond, maxCond);
805 if (c_.needMatrixACholeskyL)
806 c_.matrixACholeskyL.at(iCell) = A.llt().matrixL();
807 DNDS_assert(AInv.allFinite());
808
809 // * get AInvB and AInvb
810 for (int ic2f = 0; ic2f < mesh->cell2face.RowSize(iCell); ic2f++)
811 {
812 index iFace = mesh->cell2face(iCell, ic2f);
813 auto qFace = this->GetFaceQuad(iFace);
814 index iCellOther = CellFaceOther(iCell, iFace);
815 int if2c = CellIsFaceBack(iCell, iFace) ? 0 : 1;
816 auto faceID = mesh->GetFaceZone(iFace);
817 if (FaceIDIsExternalBC(mesh->GetFaceZone(iFace)))
818 {
819 DNDS_assert(iCellOther == UnInitIndex);
820 continue;
821 }
822 c_.matrixAAInvB(iCell, 1 + ic2f) = AInv * local_Bs.at(ic2f);
823 c_.vectorAInvB(iCell, ic2f) = AInv * local_bs.at(ic2f);
824 }
825 }
826 if (c_.needOriginalMatrix)
827 c_.matrixAB.CompressBoth();
828 c_.matrixAAInvB.CompressBoth();
829 c_.vectorAInvB.CompressBoth();
830 if (c_.needOriginalMatrix)
831 c_.vectorB.CompressBoth();
832 c_.matrixAHalf_GG.CompressBoth();
833
834 // standard comm for c_.matrixA
835 c_.matrixA.TransAttach();
836 c_.matrixA.trans.BorrowGGIndexing(mesh->cell2cell.trans);
837 c_.matrixA.trans.createMPITypes();
838 c_.matrixA.trans.pullOnce();
839
840 // Get Secondary matrices
841 this->MakePairDefaultOnFace(c_.matrixSecondary, "VR::matrixSecondary", maxNDOF - 1, (maxNDOF - 1) * 2);
842 for (index iFace = 0; iFace < mesh->NumFaceProc(); iFace++)
843 {
844 index iCellL = mesh->face2cell(iFace, 0);
845 index iCellR = mesh->face2cell(iFace, 1);
846 if (iCellR == UnInitIndex) // only for faces with two
847 continue;
848 // get dbv from 1st derivative to last
849 MatrixXR DiffI = this->GetIntPointDiffBaseValue(iCellL, iFace, 0, -1, Eigen::seq(1, EigenLast));
850 MatrixXR DiffJ = this->GetIntPointDiffBaseValue(iCellR, iFace, 1, -1, Eigen::seq(1, EigenLast));
851 MatrixXR M2_L2R, M2_R2L;
852 HardEigen::EigenLeastSquareSolve(DiffI, DiffJ, M2_R2L);
853 HardEigen::EigenLeastSquareSolve(DiffJ, DiffI, M2_L2R);
854 // TODO: cleanse M2s' lower triangle
855 int maxNDOFM1 = c_.matrixSecondary[iFace].cols() / 2;
856 c_.matrixSecondary[iFace](
857 EigenAll, Eigen::seq(
858 0 * maxNDOFM1 + 0,
859 0 * maxNDOFM1 + maxNDOFM1 - 1)) = M2_R2L;
860 c_.matrixSecondary[iFace](
861 EigenAll, Eigen::seq(
862 1 * maxNDOFM1 + 0,
863 1 * maxNDOFM1 + maxNDOFM1 - 1)) = M2_L2R;
864 }
865 c_.matrixSecondary.CompressBoth();
866
867 real maxCondAll = 0;
868 MPI::Allreduce(&maxCond, &maxCondAll, 1, DNDS_MPI_REAL, MPI_MAX, mesh->getMPI().comm);
869 if (mesh->getMPI().rank == 0)
870 log() << std::scientific << std::setprecision(3)
871 << "VariationalReconstruction<dim>::ConstructRecCoeff() === A cond Max: [ " << maxCondAll << "] " << std::endl;
872 }
873
874 template void
875 VariationalReconstruction<2>::ConstructRecCoeff();
876
877 template void
878 VariationalReconstruction<3>::ConstructRecCoeff();
879}
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
HardEigen.hpp
Robust linear-algebra primitives for small / moderately-sized Eigen matrices that are more numericall...
DNDS::CFV::VariationalReconstruction::ConstructBaseAndWeight
void ConstructBaseAndWeight(const tFGetBoundaryWeight &id2faceDircWeight=[](Geom::t_index id, int iOrder) { return 1.0;})
Definition VariationalReconstruction.cpp:56
DNDS::CFV::VariationalReconstruction::tFGetBoundaryWeight
std::function< real(Geom::t_index, int)> tFGetBoundaryWeight
Definition VariationalReconstruction.hpp:284
DNDS::Geom::Elem::GetElemNodeMajorSpan
tPoint GetElemNodeMajorSpan(const tSmallCoords &coords)
Definition Elements.hpp:513
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::Geom::FaceIDIsPeriodicDonor
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicDonor(t_index id)
Definition BoundaryCondition.hpp:130
DNDS::Geom::FaceIDIsPeriodicMain
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicMain(t_index id)
Definition BoundaryCondition.hpp:123
DNDS::Geom::tGPoint
Eigen::Matrix3d tGPoint
Definition Geometric.hpp:11
DNDS::HardEigen::EigenLeastSquareSolve
Eigen::Index EigenLeastSquareSolve(const Eigen::MatrixXd &A, const Eigen::MatrixXd &B, Eigen::MatrixXd &AIB)
Least-squares solve A * AIB ~= B via a rank-revealing QR-style decomposition; returns the computed ra...
Definition HardEigen.cpp:160
DNDS::HardEigen::EigenLeastSquareInverse_Filtered
real EigenLeastSquareInverse_Filtered(const Eigen::MatrixXd &A, Eigen::MatrixXd &AI, real svdTol, int mode)
Pseudoinverse with a choice of singular-value filter.
Definition HardEigen.cpp:25
DNDS::HardEigen::Eigen3x3RealSymEigenDecomposition
Eigen::Matrix3d Eigen3x3RealSymEigenDecomposition(const Eigen::Matrix3d &A)
Analytic eigen-decomposition of a 3x3 real symmetric matrix. Returns the eigenvector matrix (columns ...
Definition HardEigen.cpp:72
DNDS::HardEigen::Eigen2x2RealSymEigenDecomposition
Eigen::Matrix2d Eigen2x2RealSymEigenDecomposition(const Eigen::Matrix2d &A)
Analytic 2x2 analogue of Eigen3x3RealSymEigenDecomposition.
Definition HardEigen.cpp:88
DNDS::HardEigen::EigenLeastSquareInverse
real EigenLeastSquareInverse(const Eigen::MatrixXd &A, Eigen::MatrixXd &AI, real svdTol)
Moore-Penrose pseudoinverse via SVD, dropping singular values below svdTol (relative to the largest)....
Definition HardEigen.cpp:13
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::UnInitIndex
DNDS_CONSTANT const index UnInitIndex
Sentinel "not initialised" index value (= INT64_MIN).
Definition Defines.hpp:176
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:205
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::RowVectorXR
Eigen::RowVector< real, Eigen::Dynamic > RowVectorXR
Dynamic row-vector of reals.
Definition Defines.hpp:209
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::VectorXR
Eigen::Vector< real, Eigen::Dynamic > VectorXR
Dynamic Eigen vector of reals.
Definition Defines.hpp:207
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
r
tVec r(NCells)
b
tVec b(NCells)