From 8e0f93d344c9e3e41d00168a9197392af764b5a0 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 16 Dec 2013 11:53:21 +0000
Subject: [PATCH] ENH: primitiveMesh: apply stabilisation to cell volume, not
 to individual pyr vols

---
 .../primitiveMeshCellCentresAndVols.C         | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C
index 41e5366c9ae..bd58fe39a99 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,7 +111,7 @@ void Foam::primitiveMesh::makeCellCentresAndVols
     {
         // Calculate 3*face-pyramid volume
         scalar pyr3Vol =
-            max(fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]), VSMALL);
+            fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]);
 
         // Calculate face-pyramid centre
         vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[own[facei]];
@@ -127,7 +127,7 @@ void Foam::primitiveMesh::makeCellCentresAndVols
     {
         // Calculate 3*face-pyramid volume
         scalar pyr3Vol =
-            max(fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]), VSMALL);
+            fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei]);
 
         // Calculate face-pyramid centre
         vector pc = (3.0/4.0)*fCtrs[facei] + (1.0/4.0)*cEst[nei[facei]];
@@ -139,7 +139,18 @@ void Foam::primitiveMesh::makeCellCentresAndVols
         cellVols[nei[facei]] += pyr3Vol;
     }
 
-    cellCtrs /= cellVols;
+    forAll(cellCtrs, celli)
+    {
+        if (cellVols[celli] > VSMALL)
+        {
+            cellCtrs[celli] /= cellVols[celli];
+        }
+        else
+        {
+            cellCtrs[celli] = cEst[celli];
+        }
+    }
+
     cellVols *= (1.0/3.0);
 }
 
-- 
GitLab