46 averageNeighbourFvGeometryScheme,
56 const scalar minRatio,
69 for (label facei = 0; facei <
mesh_.
nFaces(); facei++)
74 const vector& fcCorr = faceCorrection[facei];
82 const vector& fn = faceNormals[facei];
83 const point& fc = faceCentres[facei];
99 faceCorrection[facei] = vector::zero;
117 faceCorrection[facei] = vector::zero;
140 ownHeight.setSize(mesh_.nFaces());
141 neiHeight.setSize(mesh_.nInternalFaces());
143 const labelList& own = mesh_.faceOwner();
144 const labelList& nei = mesh_.faceNeighbour();
146 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
153 ownHeight[facei] = ((fc-ownCc)&
n);
154 neiHeight[facei] = ((neiCc-fc)&
n);
157 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
163 ownHeight[facei] = ((fc-ownCc)&
n);
183 const labelList& own = mesh_.faceOwner();
184 const labelList& nei = mesh_.faceNeighbour();
187 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
193 const vector&
n = faceNormals[facei];
194 const point& fc = faceCentres[facei];
197 const label ownCelli = own[facei];
201 const scalar ownHeight = ((fc-ownCc)&
n);
202 if (ownHeight < minOwnHeight[facei])
214 const label neiCelli = nei[facei];
218 const scalar neiHeight = ((neiCc-fc)&
n);
219 if (neiHeight < minNeiHeight[facei])
232 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
238 const vector&
n = faceNormals[facei];
239 const point& fc = faceCentres[facei];
242 const label ownCelli = own[facei];
246 const scalar ownHeight = ((fc-ownCc)&
n);
247 if (ownHeight < minOwnHeight[facei])
271 const labelList& own = mesh_.faceOwner();
272 const labelList& nei = mesh_.faceNeighbour();
278 Field<solveScalar> cellWeights(mesh_.nCells(),
Zero);
281 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
286 const vector&
n = faceNormals[facei];
288 const point& ownCc = cellCentres[own[facei]];
289 const point& neiCc = cellCentres[nei[facei]];
303 const scalar w = 0.5*faceWeights[facei];
304 cc[own[facei]] +=
point(w*d);
305 cellWeights[own[facei]] += w;
307 cc[nei[facei]] -=
point(w*d);
308 cellWeights[nei[facei]] += w;
316 const polyBoundaryMesh&
pbm = mesh_.boundaryMesh();
317 for (
const polyPatch&
pp :
pbm)
325 const label meshFacei =
pp.start()+i;
326 const label bFacei = meshFacei-mesh_.nInternalFaces();
331 const vector&
n = faceNormals[meshFacei];
334 const point& ownCc = cellCentres[fc[i]];
335 const point& neiCc = neiCellCentres[bFacei];
348 const scalar w = 0.5*faceWeights[meshFacei];
349 cc[fc[i]] +=
point(w*d);
350 cellWeights[fc[i]] += w;
358 if (cellWeights[celli] > VSMALL)
360 cc[celli] = cellCentres[celli] + cc[celli]/cellWeights[celli];
364 cc[celli] = cellCentres[celli];
381 const labelList& own = mesh_.faceOwner();
382 const labelList& nei = mesh_.faceNeighbour();
389 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
395 const vector&
n = faceNormals[facei];
396 const point& oldFc = faceCentres[facei];
416 const solveScalar
f = (deltaFc&
n)/(deltaCc&
n);
428 newFc[facei] = oldFc + d;
436 const polyBoundaryMesh&
pbm = mesh_.boundaryMesh();
437 for (
const polyPatch&
pp :
pbm)
446 const label facei =
pp.start()+i;
447 const label bFacei = facei-mesh_.nInternalFaces();
453 const vector&
n = faceNormals[facei];
454 const point& oldFc = faceCentres[facei];
464 const solveScalar
f = (deltaFc&
n)/(deltaCc&
n);
472 newFc[facei] = oldFc + d;
480 const label facei =
pp.start()+i;
486 const vector&
n = faceNormals[facei];
487 const point& oldFc = faceCentres[facei];
497 newFc[facei] = oldFc+d;
523 faceWeights.setSize(cosAngles.size());
530 const scalar cosAngle = cosAngles[facei];
531 if (cosAngle < minCos)
533 faceWeights[facei] = 1.0;
535 else if (cosAngle > maxCos)
537 faceWeights[facei] = 0.0;
542 1.0-(cosAngle-minCos)/(maxCos-minCos);
551 Foam::averageNeighbourFvGeometryScheme::averageNeighbourFvGeometryScheme
554 const dictionary&
dict 557 highAspectRatioFvGeometryScheme(
mesh,
dict),
560 dict.getCheckOrDefault<label>
564 [](label val) {
return val >= 0; }
572 [](scalar val) { return val > 0 && val <= 1; }
581 [](scalar val) { return val >= 0 && val <= 1; }
587 Pout<<
"averageNeighbourFvGeometryScheme :" 588 <<
" nIters:" << nIters_
589 <<
" relax:" << relax_
590 <<
" minRatio:" << minRatio_ <<
endl;
604 Pout<<
"averageNeighbourFvGeometryScheme::movePoints() : " 605 <<
"recalculating primitiveMesh centres" <<
endl;
623 const vectorField faceNormals(faceAreas/magFaceAreas);
630 calcAspectRatioWeights(cellWeight, faceWeight);
633 cellWeight *= relax_;
650 minOwnHeight *= minRatio_;
651 minNeiHeight *= minRatio_;
655 autoPtr<OBJstream> osPtr;
656 autoPtr<surfaceWriter> writerPtr;
663 mesh_.time().timePath()
664 /
"cellCentre_correction.obj" 667 Pout<<
"averageNeighbourFvGeometryScheme::movePoints() :" 668 <<
" writing cell centre path to " << osPtr().
name() <<
endl;
674 mesh_.time().globalPath()
676 / mesh_.pointsInstance()
686 writerPtr->useTimeDir(
true);
688 writerPtr->beginTime(mesh_.time());
694 (outputDir /
"cosAngle"),
698 writerPtr->endTime();
707 for (label iter = 0; iter < nIters_; iter++)
740 writerPtr->beginTime(instant(scalar(iter)));
741 writerPtr->write(
"cosAngles", cosAngles);
742 writerPtr->endTime();
751 Pout<<
" face:" << facei
752 <<
" at:" << mesh_.faceCentres()[facei]
753 <<
" cosAngle:" << cosAngles[facei]
754 <<
" faceWeight:" << faceWeights[facei]
760 tcc = averageNeighbourCentres
774 const label nClipped = clipPyramids
790 forAll(cellCentres, celli)
792 const point& cc = cellCentres[celli];
817 Pout<<
" iter:" << iter
818 <<
" nClipped:" << nClipped
819 <<
" average displacement:" <<
gAverage(magCorrection)
820 <<
" non-ortho angle : average:" <<
gAverage(nonOrthoAngle)
821 <<
" max:" <<
gMax(nonOrthoAngle) <<
endl;
834 vectorField faceCorrection(faceWeight*(tfc-mesh_.faceCentres()));
844 vectorField faceCentres(mesh_.faceCentres() + faceCorrection);
848 Pout<<
"averageNeighbourFvGeometryScheme::movePoints() :" 849 <<
" averageNeighbour weight" 850 <<
" max:" <<
gMax(cellWeight) <<
" min:" <<
gMin(cellWeight)
854 const fileName tp(mesh_.time().timePath());
856 OBJstream str(tp/
"averageNeighbourCellCentres.obj");
857 Pout<<
"Writing lines from old to new cell centre to " << str.
name()
859 forAll(mesh_.cellCentres(), celli)
861 const point& oldCc = mesh_.cellCentres()[celli];
862 const point& newCc = cellCentres[celli];
863 str.writeLine(oldCc, newCc);
869 const fileName tp(mesh_.time().timePath());
870 OBJstream str(tp/
"averageFaceCentres.obj");
871 Pout<<
"Writing lines from old to new face centre to " << str.
name()
873 forAll(mesh_.faceCentres(), facei)
875 const point& oldFc = mesh_.faceCentres()[facei];
876 const point& newFc = faceCentres[facei];
877 str.writeLine(oldFc, newFc);
887 std::move(faceCentres),
888 std::move(faceAreas),
889 std::move(cellCentres),
890 std::move(cellVolumes)
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix by subtracting the matrix multiplied by the current fi...
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
const polyBoundaryMesh & pbm
dimensionedScalar acos(const dimensionedScalar &ds)
const fvMesh & mesh_
Hold reference to mesh.
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Type gMin(const FieldField< Field, Type > &f)
tmp< pointField > averageCentres(const pointField &cellCentres, const pointField &faceCentres, const vectorField &faceNormals) const
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
Unit conversion functions.
T getCheck(const word &keyword, const Predicate &pred, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T with additional checking FatalIOError if not found, or if the number of tokens is...
Ostream & endl(Ostream &os)
Add newline and flush stream.
label nFaces() const noexcept
Number of mesh faces.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Macros for easy insertion into run-time selection tables.
UList< label > labelUList
A UList of labels.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
Vector< solveScalar > solveVector
void resetGeometry(pointField &&faceCentres, pointField &&faceAreas, pointField &&cellCentres, scalarField &&cellVolumes)
Reset the local geometry.
vectorField pointField
pointField is a vectorField.
virtual void movePoints()
Do what is necessary if the mesh has moved.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
dimensionedScalar cos(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
T getCheckOrDefault(const word &keyword, const T &deflt, const Predicate &pred, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
virtual const faceList & faces() const
Return raw faces.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
tmp< pointField > averageNeighbourCentres(const pointField &cellCentres, const vectorField &faceNormals, const scalarField &faceWeights) const
Average neighbouring cell centres to minimise non-ortho.
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross-product operators.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
label clipFaceTet(const scalar minRatio, const vectorField &faceCentres, const vectorField &faceNormals, vectorField &faceCorrection) const
Clip face-centre correction vector if new triangle area would get below min. Return number of clipped...
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
void makeNonOrthoWeights(const pointField &cellCentres, const vectorField &faceNormals, scalarField &cosAngles, scalarField &faceWeights) const
Make weights based on non-orthogonality.
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
vector point
Point is a vector.
Type gAverage(const FieldField< Field, Type > &f)
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
label clipPyramids(const pointField &cellCentres, const vectorField &faceCentres, const vectorField &faceNormals, const scalarField &minOwnHeight, const scalarField &minNeiHeight, vectorField &correction) const
Clip correction vector if new pyramid height would get below min. Return number of clipped cells...
static word outputPrefix
Directory prefix.
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< label > labelList
A List of labels.
A class for managing temporary objects.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
virtual void movePoints()
Do what is necessary if the mesh has moved.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
static autoPtr< surfaceWriter > New(const word &writeType)
Return a reference to the selected surfaceWriter.
static const Vector< Cmpt > zero
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
dimensionSet clamp(const dimensionSet &a, const dimensionSet &range)
void makePyrHeights(const pointField &cellCentres, const vectorField &faceCentres, const vectorField &faceNormals, scalarField &ownHeight, scalarField &neiHeight) const
Calculate pyramid heights.
static constexpr const zero Zero
Global zero (0)