From 431042dc232f8f0eff5af80b71d8ad5be55b8915 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 13 Oct 2015 16:22:03 +0100
Subject: [PATCH] BUG: polyMeshGeometry: fixed updating of local geometric
 properties

---
 .../motionSmoother/motionSmootherAlgo.H       |  16 ++-
 .../motionSmoother/motionSmootherAlgoCheck.C  | 100 +++++++-------
 .../polyMeshGeometry/polyMeshGeometry.C       | 122 +++++++++---------
 3 files changed, 128 insertions(+), 110 deletions(-)

diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
index d4c1054860..724658012e 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
@@ -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
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
index 07859539a5..8928dc29bb 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
@@ -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
         );
 
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 57725dc57c..e6426ca9da 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -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
-
-    vectorField cEst(mesh_.nCells());
-    UIndirectList<vector>(cEst, changedCells) = vector::zero;
-    scalarField nCellFaces(mesh_.nCells());
-    UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
-
-    forAll(changedFaces, i)
+    // Re-calculate the changed cell centres and volumes
+    forAll(changedCells, changedCellI)
     {
-        label faceI = changedFaces[i];
-        cEst[own[faceI]] += faceCentres_[faceI];
-        nCellFaces[own[faceI]] += 1;
+        const label cellI(changedCells[changedCellI]);
 
-        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];
-            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];
-
-        // 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
-        cellVolumes_[own[faceI]] += pyr3Vol;
-
-        if (mesh_.isInternalFace(faceI))
+        // Sum up the face-pyramid contributions
+        forAll(cFaces, cFaceI)
         {
-            // Calculate 3*face-pyramid volume
-            scalar pyr3Vol = max
-            (
-                faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
-                VSMALL
-            );
+            const label faceI(cFaces[cFaceI]);
 
-            // Calculate face-pyramid centre
-            vector pc =
-                (3.0/4.0)*faceCentres_[faceI]
-              + (1.0/4.0)*cEst[nei[faceI]];
+            // Calculate 3* the face-pyramid volume
+            scalar pyr3Vol = faceAreas_[faceI] & (faceCentres_[faceI] - cEst);
 
-            // Accumulate volume-weighted face-pyramid centre
-            cellCentres_[nei[faceI]] += pyr3Vol*pc;
+            if (own[faceI] != cellI)
+            {
+                pyr3Vol = -pyr3Vol;
+            }
 
             // Accumulate face-pyramid volume
-            cellVolumes_[nei[faceI]] += pyr3Vol;
+            cellVolumes_[cellI] += pyr3Vol;
+
+            // Calculate the face-pyramid centre
+            const vector pCtr = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst;
+
+            // Accumulate volume-weighted face-pyramid centre
+            cellCentres_[cellI] += pyr3Vol*pCtr;
         }
-    }
 
-    forAll(changedCells, i)
-    {
-        label cellI = changedCells[i];
+        // Average the accumulated quantities
+
+        if (mag(cellVolumes_[cellI]) > VSMALL)
+        {
+            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);
     }
 }