BUG: polyMeshGeometry: fixed updating of local geometric properties

This commit is contained in:
mattijs
2015-10-13 16:22:03 +01:00
committed by Andrew Heather
parent d292a909aa
commit 431042dc23
3 changed files with 128 additions and 110 deletions

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -332,6 +332,17 @@ public:
const dictionary& paramDict() const; 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 // Edit
@ -339,7 +350,6 @@ public:
//- Take over existing mesh position. //- Take over existing mesh position.
void correct(); void correct();
//- Set patch fields on patchIDs to be consistent with //- Set patch fields on patchIDs to be consistent with
// all other boundary conditions // all other boundary conditions
static void setDisplacementPatchFields static void setDisplacementPatchFields
@ -471,6 +481,7 @@ public:
const bool report, const bool report,
const dictionary& dict, const dictionary& dict,
const polyMeshGeometry&, const polyMeshGeometry&,
const pointField&,
const labelList& checkFaces, const labelList& checkFaces,
labelHashSet& wrongFaces labelHashSet& wrongFaces
); );
@ -483,6 +494,7 @@ public:
const bool report, const bool report,
const dictionary& dict, const dictionary& dict,
const polyMeshGeometry&, const polyMeshGeometry&,
const pointField&,
const labelList& checkFaces, const labelList& checkFaces,
const List<labelPair>& baffles, const List<labelPair>& baffles,
labelHashSet& wrongFaces labelHashSet& wrongFaces

View File

@ -2,8 +2,8 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -434,6 +434,7 @@ bool Foam::motionSmootherAlgo::checkMesh
const bool report, const bool report,
const dictionary& dict, const dictionary& dict,
const polyMeshGeometry& meshGeom, const polyMeshGeometry& meshGeom,
const pointField& points,
const labelList& checkFaces, const labelList& checkFaces,
labelHashSet& wrongFaces labelHashSet& wrongFaces
) )
@ -445,6 +446,7 @@ bool Foam::motionSmootherAlgo::checkMesh
report, report,
dict, dict,
meshGeom, meshGeom,
points,
checkFaces, checkFaces,
emptyBaffles, emptyBaffles,
wrongFaces wrongFaces
@ -457,6 +459,7 @@ bool Foam::motionSmootherAlgo::checkMesh
const bool report, const bool report,
const dictionary& dict, const dictionary& dict,
const polyMeshGeometry& meshGeom, const polyMeshGeometry& meshGeom,
const pointField& points,
const labelList& checkFaces, const labelList& checkFaces,
const List<labelPair>& baffles, const List<labelPair>& baffles,
labelHashSet& wrongFaces labelHashSet& wrongFaces
@ -482,14 +485,14 @@ bool Foam::motionSmootherAlgo::checkMesh
( (
readScalar(dict.lookup("minArea", true)) readScalar(dict.lookup("minArea", true))
); );
//const scalar maxIntSkew const scalar maxIntSkew
//( (
// readScalar(dict.lookup("maxInternalSkewness", true)) readScalar(dict.lookup("maxInternalSkewness", true))
//); );
//const scalar maxBounSkew const scalar maxBounSkew
//( (
// readScalar(dict.lookup("maxBoundarySkewness", true)) readScalar(dict.lookup("maxBoundarySkewness", true))
//); );
const scalar minWeight const scalar minWeight
( (
readScalar(dict.lookup("minFaceWeight", true)) readScalar(dict.lookup("minFaceWeight", true))
@ -512,7 +515,6 @@ bool Foam::motionSmootherAlgo::checkMesh
( (
readScalar(dict.lookup("minDeterminant", true)) readScalar(dict.lookup("minDeterminant", true))
); );
label nWrongFaces = 0; label nWrongFaces = 0;
Info<< "Checking faces in error :" << endl; Info<< "Checking faces in error :" << endl;
@ -544,8 +546,8 @@ bool Foam::motionSmootherAlgo::checkMesh
meshGeom.checkFacePyramids meshGeom.checkFacePyramids
( (
report, report,
minVol, minTetQuality,
meshGeom.mesh().points(), points,
checkFaces, checkFaces,
baffles, baffles,
&wrongFaces &wrongFaces
@ -566,7 +568,7 @@ bool Foam::motionSmootherAlgo::checkMesh
( (
report, report,
minTetQuality, minTetQuality,
meshGeom.mesh().points(), points,
checkFaces, checkFaces,
baffles, baffles,
&wrongFaces &wrongFaces
@ -587,7 +589,7 @@ bool Foam::motionSmootherAlgo::checkMesh
( (
report, report,
maxConcave, maxConcave,
meshGeom.mesh().points(), points,
checkFaces, checkFaces,
&wrongFaces &wrongFaces
); );
@ -604,7 +606,13 @@ bool Foam::motionSmootherAlgo::checkMesh
if (minArea > -SMALL) if (minArea > -SMALL)
{ {
meshGeom.checkFaceArea(report, minArea, checkFaces, &wrongFaces); meshGeom.checkFaceArea
(
report,
minArea,
checkFaces,
&wrongFaces
);
label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>()); label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
@ -616,30 +624,32 @@ bool Foam::motionSmootherAlgo::checkMesh
nWrongFaces = nNewWrongFaces; 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 label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
// to store them on polyMeshGeometry.
//if (maxIntSkew > 0 || maxBounSkew > 0) Info<< " faces with skewness > "
//{ << setw(3) << maxIntSkew
// meshGeom.checkFaceSkewness << " (internal) or " << setw(3) << maxBounSkew
// ( << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
// report,
// maxIntSkew, nWrongFaces = nNewWrongFaces;
// 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;
//}
if (minWeight >= 0 && minWeight < 1) if (minWeight >= 0 && minWeight < 1)
{ {
@ -691,7 +701,7 @@ bool Foam::motionSmootherAlgo::checkMesh
( (
report, report,
minTwist, minTwist,
meshGeom.mesh().points(), points,
checkFaces, checkFaces,
&wrongFaces &wrongFaces
); );
@ -714,7 +724,7 @@ bool Foam::motionSmootherAlgo::checkMesh
( (
report, report,
minTriangleTwist, minTriangleTwist,
meshGeom.mesh().points(), points,
checkFaces, checkFaces,
&wrongFaces &wrongFaces
); );
@ -729,13 +739,13 @@ bool Foam::motionSmootherAlgo::checkMesh
nWrongFaces = nNewWrongFaces; nWrongFaces = nNewWrongFaces;
} }
if (minFaceFlatness > -1) if (minFaceFlatness > -SMALL)
{ {
meshGeom.checkFaceFlatness meshGeom.checkFaceFlatness
( (
report, report,
minFaceFlatness, minFaceFlatness,
meshGeom.mesh().points(), points,
checkFaces, checkFaces,
&wrongFaces &wrongFaces
); );
@ -757,7 +767,7 @@ bool Foam::motionSmootherAlgo::checkMesh
report, report,
minDet, minDet,
checkFaces, checkFaces,
meshGeom.affectedCells(meshGeom.mesh(), checkFaces), polyMeshGeometry::affectedCells(meshGeom.mesh(), checkFaces),
&wrongFaces &wrongFaces
); );

View File

@ -2,8 +2,8 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -105,86 +105,82 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
const labelList& changedFaces const labelList& changedFaces
) )
{ {
const labelList& own = mesh().faceOwner();
const cellList& cells = mesh().cells();
// Clear the fields for accumulation // Clear the fields for accumulation
UIndirectList<vector>(cellCentres_, changedCells) = vector::zero; UIndirectList<vector>(cellCentres_, changedCells) = vector::zero;
UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0; 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)
vectorField cEst(mesh_.nCells());
UIndirectList<vector>(cEst, changedCells) = vector::zero;
scalarField nCellFaces(mesh_.nCells());
UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
forAll(changedFaces, i)
{ {
label faceI = changedFaces[i]; const label cellI(changedCells[changedCellI]);
cEst[own[faceI]] += faceCentres_[faceI];
nCellFaces[own[faceI]] += 1;
if (mesh_.isInternalFace(faceI)) const labelList& cFaces(cells[cellI]);
// 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]; const point& fc = faceCentres_[cFaces[cFaceI]];
nCellFaces[nei[faceI]] += 1; cEst += fc;
} bb.max() = max(bb.max(), fc);
bb.min() = min(bb.min(), fc);
} }
cEst /= cFaces.size();
forAll(changedCells, i)
// Sum up the face-pyramid contributions
forAll(cFaces, cFaceI)
{ {
label cellI = changedCells[i]; const label faceI(cFaces[cFaceI]);
cEst[cellI] /= nCellFaces[cellI];
}
forAll(changedFaces, i) // Calculate 3* the face-pyramid volume
scalar pyr3Vol = faceAreas_[faceI] & (faceCentres_[faceI] - cEst);
if (own[faceI] != cellI)
{ {
label faceI = changedFaces[i]; pyr3Vol = -pyr3Vol;
}
// Calculate 3*face-pyramid volume
scalar pyr3Vol = max
(
faceAreas_[faceI] & (faceCentres_[faceI] - cEst[own[faceI]]),
VSMALL
);
// Calculate face-pyramid centre
vector pc = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst[own[faceI]];
// Accumulate volume-weighted face-pyramid centre
cellCentres_[own[faceI]] += pyr3Vol*pc;
// Accumulate face-pyramid volume // Accumulate face-pyramid volume
cellVolumes_[own[faceI]] += pyr3Vol; cellVolumes_[cellI] += pyr3Vol;
if (mesh_.isInternalFace(faceI)) // Calculate the face-pyramid centre
{ const vector pCtr = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst;
// Calculate 3*face-pyramid volume
scalar pyr3Vol = max
(
faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
VSMALL
);
// Calculate face-pyramid centre
vector pc =
(3.0/4.0)*faceCentres_[faceI]
+ (1.0/4.0)*cEst[nei[faceI]];
// Accumulate volume-weighted face-pyramid centre // Accumulate volume-weighted face-pyramid centre
cellCentres_[nei[faceI]] += pyr3Vol*pc; cellCentres_[cellI] += pyr3Vol*pCtr;
// Accumulate face-pyramid volume
cellVolumes_[nei[faceI]] += pyr3Vol;
}
} }
forAll(changedCells, i) // Average the accumulated quantities
if (mag(cellVolumes_[cellI]) > VSMALL)
{ {
label cellI = changedCells[i]; point cc = cellCentres_[cellI] / cellVolumes_[cellI];
// 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;
}
cellCentres_[cellI] /= cellVolumes_[cellI] + VSMALL;
cellVolumes_[cellI] *= (1.0/3.0); cellVolumes_[cellI] *= (1.0/3.0);
} }
} }