Skip to content

GitLab

Commit 43c92bb0 authored by mattijs's avatar mattijs
Browse files

BUG: polyMeshGeometry: fixed updating of local geometric properties

parent b4f78778
Hide whitespace changes
Inline Side-by-side
src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
View file @ 43c92bb0
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -332,6 +332,17 @@ public:
const dictionary& paramDict() const;
//- Return reference to the point motion displacement field
pointVectorField& pointDisplacement()
{
return displacement_;
}
//- Return const reference to the point motion displacement field
const pointVectorField& pointDisplacement() const
{
return displacement_;
}
// Edit
......@@ -339,7 +350,6 @@ public:
//- Take over existing mesh position.
void correct();
//- Set patch fields on patchIDs to be consistent with
// all other boundary conditions
static void setDisplacementPatchFields
......@@ -471,6 +481,7 @@ public:
const bool report,
const dictionary& dict,
const polyMeshGeometry&,
const pointField&,
const labelList& checkFaces,
labelHashSet& wrongFaces
);
......@@ -483,6 +494,7 @@ public:
const bool report,
const dictionary& dict,
const polyMeshGeometry&,
const pointField&,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
......
src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
View file @ 43c92bb0
......@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -434,6 +434,7 @@ bool Foam::motionSmootherAlgo::checkMesh
const bool report,
const dictionary& dict,
const polyMeshGeometry& meshGeom,
const pointField& points,
const labelList& checkFaces,
labelHashSet& wrongFaces
)
......@@ -445,6 +446,7 @@ bool Foam::motionSmootherAlgo::checkMesh
report,
dict,
meshGeom,
points,
checkFaces,
emptyBaffles,
wrongFaces
......@@ -457,6 +459,7 @@ bool Foam::motionSmootherAlgo::checkMesh
const bool report,
const dictionary& dict,
const polyMeshGeometry& meshGeom,
const pointField& points,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
......@@ -482,14 +485,14 @@ bool Foam::motionSmootherAlgo::checkMesh
(
readScalar(dict.lookup("minArea", true))
);
//const scalar maxIntSkew
//(
// readScalar(dict.lookup("maxInternalSkewness", true))
//);
//const scalar maxBounSkew
//(
// readScalar(dict.lookup("maxBoundarySkewness", true))
//);
const scalar maxIntSkew
(
readScalar(dict.lookup("maxInternalSkewness", true))
);
const scalar maxBounSkew
(
readScalar(dict.lookup("maxBoundarySkewness", true))
);
const scalar minWeight
(
readScalar(dict.lookup("minFaceWeight", true))
......@@ -512,7 +515,6 @@ bool Foam::motionSmootherAlgo::checkMesh
(
readScalar(dict.lookup("minDeterminant", true))
);
label nWrongFaces = 0;
Info<< "Checking faces in error :" << endl;
......@@ -544,8 +546,8 @@ bool Foam::motionSmootherAlgo::checkMesh
meshGeom.checkFacePyramids
(
report,
minVol,
meshGeom.mesh().points(),
minTetQuality,
points,
checkFaces,
baffles,
&wrongFaces
......@@ -566,7 +568,7 @@ bool Foam::motionSmootherAlgo::checkMesh
(
report,
minTetQuality,
meshGeom.mesh().points(),
points,
checkFaces,
baffles,
&wrongFaces
......@@ -575,7 +577,7 @@ bool Foam::motionSmootherAlgo::checkMesh
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with face-decomposition tet quality < "
<< setw(5) << minTetQuality << " : "
<< setw(5) << minTetQuality << " : "
<< nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
......@@ -587,7 +589,7 @@ bool Foam::motionSmootherAlgo::checkMesh
(
report,
maxConcave,
meshGeom.mesh().points(),
points,
checkFaces,
&wrongFaces
);
......@@ -604,7 +606,13 @@ bool Foam::motionSmootherAlgo::checkMesh
if (minArea > -SMALL)
{
meshGeom.checkFaceArea(report, minArea, checkFaces, &wrongFaces);
meshGeom.checkFaceArea
(
report,
minArea,
checkFaces,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
......@@ -616,30 +624,32 @@ bool Foam::motionSmootherAlgo::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (maxIntSkew > 0 || maxBounSkew > 0)
{
polyMeshGeometry::checkFaceSkewness
(
report,
maxIntSkew,
maxBounSkew,
meshGeom.mesh(),
points,
meshGeom.cellCentres(),
meshGeom.faceCentres(),
meshGeom.faceAreas(),
checkFaces,
baffles,
&wrongFaces
);
//- Note: cannot check the skewness without the points and don't want
// to store them on polyMeshGeometry.
//if (maxIntSkew > 0 || maxBounSkew > 0)
//{
// meshGeom.checkFaceSkewness
// (
// report,
// maxIntSkew,
// maxBounSkew,
// checkFaces,
// baffles,
// &wrongFaces
// );
//
// label nNewWrongFaces = returnReduce(wrongFaces.size(),sumOp<label>());
//
// Info<< " faces with skewness > "
// << setw(3) << maxIntSkew
// << " (internal) or " << setw(3) << maxBounSkew
// << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
//
// nWrongFaces = nNewWrongFaces;
//}
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
Info<< " faces with skewness > "
<< setw(3) << maxIntSkew
<< " (internal) or " << setw(3) << maxBounSkew
<< " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
nWrongFaces = nNewWrongFaces;
}
if (minWeight >= 0 && minWeight < 1)
{
......@@ -691,7 +701,7 @@ bool Foam::motionSmootherAlgo::checkMesh
(
report,
minTwist,
meshGeom.mesh().points(),
points,
checkFaces,
&wrongFaces
);
......@@ -714,7 +724,7 @@ bool Foam::motionSmootherAlgo::checkMesh
(
report,
minTriangleTwist,
meshGeom.mesh().points(),
points,
checkFaces,
&wrongFaces
);
......@@ -729,13 +739,13 @@ bool Foam::motionSmootherAlgo::checkMesh
nWrongFaces = nNewWrongFaces;
}
if (minFaceFlatness > -1)
if (minFaceFlatness > -SMALL)
{
meshGeom.checkFaceFlatness
(
report,
minFaceFlatness,
meshGeom.mesh().points(),
points,
checkFaces,
&wrongFaces
);
......@@ -757,7 +767,7 @@ bool Foam::motionSmootherAlgo::checkMesh
report,
minDet,
checkFaces,
meshGeom.affectedCells(meshGeom.mesh(), checkFaces),
polyMeshGeometry::affectedCells(meshGeom.mesh(), checkFaces),
&wrongFaces
);
......
src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
View file @ 43c92bb0
......@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -105,86 +105,82 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
const labelList& changedFaces
)
{
const labelList& own = mesh().faceOwner();
const cellList& cells = mesh().cells();
// Clear the fields for accumulation
UIndirectList<vector>(cellCentres_, changedCells) = vector::zero;
UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
// first estimate the approximate cell centre as the average of face centres
// Re-calculate the changed cell centres and volumes
forAll(changedCells, changedCellI)
{
const label cellI(changedCells[changedCellI]);
vectorField cEst(mesh_.nCells());
UIndirectList<vector>(cEst, changedCells) = vector::zero;
scalarField nCellFaces(mesh_.nCells());
UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
const labelList& cFaces(cells[cellI]);
forAll(changedFaces, i)
{
label faceI = changedFaces[i];
cEst[own[faceI]] += faceCentres_[faceI];
nCellFaces[own[faceI]] += 1;
if (mesh_.isInternalFace(faceI))
// Estimate the cell centre and bounding box using the face centres
vector cEst = vector::zero;
boundBox bb(boundBox::invertedBox);
forAll(cFaces, cFaceI)
{
cEst[nei[faceI]] += faceCentres_[faceI];
nCellFaces[nei[faceI]] += 1;
const point& fc = faceCentres_[cFaces[cFaceI]];
cEst += fc;
bb.max() = max(bb.max(), fc);
bb.min() = min(bb.min(), fc);
}
}
cEst /= cFaces.size();
forAll(changedCells, i)
{
label cellI = changedCells[i];
cEst[cellI] /= nCellFaces[cellI];
}
forAll(changedFaces, i)
{
label faceI = changedFaces[i];
// Sum up the face-pyramid contributions
forAll(cFaces, cFaceI)
{
const label faceI(cFaces[cFaceI]);
// Calculate 3*face-pyramid volume
scalar pyr3Vol = max
(
faceAreas_[faceI] & (faceCentres_[faceI] - cEst[own[faceI]]),
VSMALL
);
// Calculate 3* the face-pyramid volume
scalar pyr3Vol = faceAreas_[faceI] & (faceCentres_[faceI] - cEst);
// Calculate face-pyramid centre
vector pc = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst[own[faceI]];
if (own[faceI] != cellI)
{
pyr3Vol = -pyr3Vol;
}
// Accumulate volume-weighted face-pyramid centre
cellCentres_[own[faceI]] += pyr3Vol*pc;
// Accumulate face-pyramid volume
cellVolumes_[cellI] += pyr3Vol;
// Accumulate face-pyramid volume
cellVolumes_[own[faceI]] += pyr3Vol;
// Calculate the face-pyramid centre
const vector pCtr = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst;
if (mesh_.isInternalFace(faceI))
{
// Calculate 3*face-pyramid volume
scalar pyr3Vol = max
(
faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
VSMALL
);
// Accumulate volume-weighted face-pyramid centre
cellCentres_[cellI] += pyr3Vol*pCtr;
}
// Calculate face-pyramid centre
vector pc =
(3.0/4.0)*faceCentres_[faceI]
+ (1.0/4.0)*cEst[nei[faceI]];
// Average the accumulated quantities
// Accumulate volume-weighted face-pyramid centre
cellCentres_[nei[faceI]] += pyr3Vol*pc;
if (mag(cellVolumes_[cellI]) > VSMALL)
{
point cc = cellCentres_[cellI] / cellVolumes_[cellI];
// Accumulate face-pyramid volume
cellVolumes_[nei[faceI]] += pyr3Vol;
// Do additional check for collapsed cells since some volumes
// (e.g. 1e-33) do not trigger above but do return completely
// wrong cell centre
if (bb.contains(cc))
{
cellCentres_[cellI] = cc;
}
else
{
cellCentres_[cellI] = cEst;
}
}
else
{
cellCentres_[cellI] = cEst;
}
}
forAll(changedCells, i)
{
label cellI = changedCells[i];
cellCentres_[cellI] /= cellVolumes_[cellI] + VSMALL;
cellVolumes_[cellI] *= (1.0/3.0);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment