46 averageNeighbourFvGeometryScheme,
56 const scalar minRatio,
71 for (label facei = 0; facei < mesh_.nFaces(); facei++)
74 const solveVector fcCorr(faceCorrection[facei]);
76 const vector& fcCorr = faceCorrection[facei];
78 if (fcCorr != solveVector::zero)
82 const solveVector fc(faceCentres[facei]);
85 const point& fc = faceCentres[facei];
87 const face&
f = mesh_.faces()[facei];
91 const solveVector thisPt(
p[
f[fp]]);
92 const solveVector nextPt(
p[
f.fcValue(fp)]);
93 const solveVector d(nextPt-thisPt);
96 const solveVector
nCorr(d^(fc+fcCorr - thisPt));
101 faceCorrection[facei] = vector::zero;
108 const solveVector
n(d^(fc - thisPt));
119 faceCorrection[facei] = vector::zero;
142 ownHeight.setSize(mesh_.nFaces());
143 neiHeight.setSize(mesh_.nInternalFaces());
147 const labelList& own = mesh_.faceOwner();
148 const labelList& nei = mesh_.faceNeighbour();
150 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
153 const solveVector fc = faceCentres[facei];
154 const solveVector ownCc = cellCentres[own[facei]];
155 const solveVector neiCc = cellCentres[nei[facei]];
157 ownHeight[facei] = ((fc-ownCc)&
n);
158 neiHeight[facei] = ((neiCc-fc)&
n);
161 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
164 const solveVector fc = faceCentres[facei];
165 const solveVector ownCc = cellCentres[own[facei]];
167 ownHeight[facei] = ((fc-ownCc)&
n);
189 const labelList& own = mesh_.faceOwner();
190 const labelList& nei = mesh_.faceNeighbour();
193 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
197 const solveVector fc(faceCentres[facei]);
200 const point& fc = faceCentres[facei];
203 const label ownCelli = own[facei];
206 const solveVector ownCc(cellCentres[ownCelli]+
correction[ownCelli]);
207 const scalar ownHeight = ((fc-ownCc)&
n);
208 if (ownHeight < minOwnHeight[facei])
220 const label neiCelli = nei[facei];
223 const solveVector neiCc(cellCentres[neiCelli]+
correction[neiCelli]);
224 const scalar neiHeight = ((neiCc-fc)&
n);
225 if (neiHeight < minNeiHeight[facei])
238 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
242 const solveVector fc(faceCentres[facei]);
245 const point& fc = faceCentres[facei];
248 const label ownCelli = own[facei];
251 const solveVector ownCc(cellCentres[ownCelli]+
correction[ownCelli]);
252 const scalar ownHeight = ((fc-ownCc)&
n);
253 if (ownHeight < minOwnHeight[facei])
279 const labelList& own = mesh_.faceOwner();
280 const labelList& nei = mesh_.faceNeighbour();
289 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
296 const point& ownCc = cellCentres[own[facei]];
297 const point& neiCc = cellCentres[nei[facei]];
299 solveVector d(neiCc-ownCc);
311 const scalar w = 0.5*faceWeights[facei];
312 cc[own[facei]] +=
point(w*d);
313 cellWeights[own[facei]] += w;
315 cc[nei[facei]] -=
point(w*d);
316 cellWeights[nei[facei]] += w;
322 syncTools::swapBoundaryCellPositions(mesh_, cellCentres, neiCellCentres);
333 const label meshFacei = pp.start()+i;
334 const label bFacei = meshFacei-mesh_.nInternalFaces();
342 const point& ownCc = cellCentres[fc[i]];
343 const point& neiCc = neiCellCentres[bFacei];
345 solveVector d(neiCc-ownCc);
356 const scalar w = 0.5*faceWeights[meshFacei];
357 cc[fc[i]] +=
point(w*d);
358 cellWeights[fc[i]] += w;
366 if (cellWeights[celli] > VSMALL)
368 cc[celli] = cellCentres[celli] + cc[celli]/cellWeights[celli];
372 cc[celli] = cellCentres[celli];
391 const labelList& own = mesh_.faceOwner();
392 const labelList& nei = mesh_.faceNeighbour();
399 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
403 const solveVector oldFc(faceCentres[facei]);
406 const point& oldFc = faceCentres[facei];
409 const solveVector ownCc(cellCentres[own[facei]]);
410 const solveVector neiCc(cellCentres[nei[facei]]);
412 solveVector deltaCc(neiCc-ownCc);
413 solveVector deltaFc(oldFc-ownCc);
426 const solveScalar
f = (deltaFc&
n)/(deltaCc&
n);
427 const solveVector avgCc((1.0-
f)*ownCc +
f*neiCc);
429 solveVector d(avgCc-oldFc);
438 newFc[facei] = oldFc + d;
444 syncTools::swapBoundaryCellPositions(mesh_, cellCentres, neiCellCentres);
456 const label facei = pp.start()+i;
457 const label bFacei = facei-mesh_.nInternalFaces();
461 const solveVector oldFc(faceCentres[facei]);
464 const point& oldFc = faceCentres[facei];
467 const solveVector ownCc(cellCentres[fc[i]]);
468 const solveVector neiCc(neiCellCentres[bFacei]);
470 solveVector deltaCc(neiCc-ownCc);
471 solveVector deltaFc(oldFc-ownCc);
474 const solveScalar
f = (deltaFc&
n)/(deltaCc&
n);
475 const solveVector avgCc((1.0-
f)*ownCc +
f*neiCc);
477 solveVector d(avgCc-oldFc);
482 newFc[facei] = oldFc + d;
490 const label facei = pp.start()+i;
494 const solveVector oldFc(faceCentres[facei]);
497 const point& oldFc = faceCentres[facei];
500 const solveVector ownCc(cellCentres[fc[i]]);
502 solveVector d(ownCc-oldFc);
507 newFc[facei] = oldFc+d;
532 polyMeshTools::faceOrthogonality
544 faceWeights.setSize(cosAngles.size());
551 const scalar cosAngle = cosAngles[facei];
552 if (cosAngle < minCos)
554 faceWeights[facei] = 1.0;
556 else if (cosAngle > maxCos)
558 faceWeights[facei] = 0.0;
563 1.0-(cosAngle-minCos)/(maxCos-minCos);
572 Foam::averageNeighbourFvGeometryScheme::averageNeighbourFvGeometryScheme
581 dict.getCheckOrDefault<label>
585 [&](
const label& nIters)
593 dict.getCheck<scalar>
596 [&](
const scalar&
relax)
598 return relax > 0 && relax <= 1;
604 dict.getCheckOrDefault<scalar>
608 [&](
const scalar& minRatio)
610 return minRatio >= 0 && minRatio <= 1;
617 Pout<<
"averageNeighbourFvGeometryScheme :"
618 <<
" nIters:" << nIters_
619 <<
" relax:" << relax_
620 <<
" minRatio:" << minRatio_ <<
endl;
634 Pout<<
"averageNeighbourFvGeometryScheme::movePoints() : "
635 <<
"recalculating primitiveMesh centres" <<
endl;
694 /
"cellCentre_correction.obj"
697 Pout<<
"averageNeighbourFvGeometryScheme::movePoints() :"
698 <<
" writing cell centre path to " << osPtr().name() <<
endl;
716 writerPtr->useTimeDir(
true);
724 (outputDir /
"cosAngle"),
728 writerPtr->endTime();
737 for (label iter = 0; iter <
nIters_; iter++)
770 writerPtr->beginTime(
instant(scalar(iter)));
771 writerPtr->write(
"cosAngles", cosAngles);
772 writerPtr->endTime();
781 Pout<<
" face:" << facei
783 <<
" cosAngle:" << cosAngles[facei]
784 <<
" faceWeight:" << faceWeights[facei]
820 forAll(cellCentres, celli)
822 const point& cc = cellCentres[celli];
847 Pout<<
" iter:" << iter
848 <<
" nClipped:" << nClipped
849 <<
" average displacement:" <<
gAverage(magCorrection)
850 <<
" non-ortho angle : average:" <<
gAverage(nonOrthoAngle)
851 <<
" max:" <<
gMax(nonOrthoAngle) <<
endl;
878 Pout<<
"averageNeighbourFvGeometryScheme::movePoints() :"
879 <<
" averageNeighbour weight"
880 <<
" max:" <<
gMax(cellWeight) <<
" min:" <<
gMin(cellWeight)
886 OBJstream str(tp/
"averageNeighbourCellCentres.obj");
887 Pout<<
"Writing lines from old to new cell centre to " << str.
name()
892 const point& newCc = cellCentres[celli];
900 OBJstream str(tp/
"averageFaceCentres.obj");
901 Pout<<
"Writing lines from old to new face centre to " << str.
name()
906 const point& newFc = faceCentres[facei];
917 std::move(faceCentres),
918 std::move(faceAreas),
919 std::move(cellCentres),
920 std::move(cellVolumes)