STYLE: medialAxisMover: moved routines into meshRefinement

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C

@@ -51,7 +51,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-
 Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
 (
     const pointVectorField& fld
@@ -118,63 +117,6 @@ Foam::medialAxisMeshMover::getPatch
 }
 
 
-void Foam::medialAxisMeshMover::calculateEdgeWeights
-(
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    scalarField& edgeWeights,
-    scalarField& invSumWeight
-) const
-{
-    const pointField& pts = mesh().points();
-
-    // Calculate edgeWeights and inverse sum of edge weights
-    edgeWeights.setSize(meshEdges.size());
-    invSumWeight.setSize(meshPoints.size());
-
-    forAll(edges, edgeI)
-    {
-        const edge& e = edges[edgeI];
-        scalar eMag = max
-        (
-            VSMALL,
-            mag
-            (
-                pts[meshPoints[e[1]]]
-              - pts[meshPoints[e[0]]]
-            )
-        );
-        edgeWeights[edgeI] = 1.0/eMag;
-    }
-
-    // Sum per point all edge weights
-    weightedSum
-    (
-        mesh(),
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        edgeWeights,
-        scalarField(meshPoints.size(), 1.0),  // data
-        invSumWeight
-    );
-
-    // Inplace invert
-    forAll(invSumWeight, pointI)
-    {
-        scalar w = invSumWeight[pointI];
-
-        if (w > 0.0)
-        {
-            invSumWeight[pointI] = 1.0/w;
-        }
-    }
-}
-
-
 void Foam::medialAxisMeshMover::smoothPatchNormals
 (
     const label nSmoothDisp,
@@ -193,8 +135,9 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
 
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -207,7 +150,7 @@ void Foam::medialAxisMeshMover::smoothPatchNormals
     vectorField average;
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -284,8 +227,9 @@ void Foam::medialAxisMeshMover::smoothNormals
 
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -297,7 +241,7 @@ void Foam::medialAxisMeshMover::smoothNormals
     vectorField average;
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -1059,8 +1003,9 @@ void Foam::medialAxisMeshMover::minSmoothField
 
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -1075,7 +1020,7 @@ void Foam::medialAxisMeshMover::minSmoothField
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
         scalarField average(pp.nPoints());
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -1534,8 +1479,9 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
     // Calculate inverse sum of weights
     scalarField edgeWeights(meshEdges.size());
     scalarField invSumWeight(meshPoints.size());
-    calculateEdgeWeights
+    meshRefinement::calculateEdgeWeights
     (
+        mesh(),
         isMasterEdge,
         meshEdges,
         meshPoints,
@@ -1554,7 +1500,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
 
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,
@@ -1575,7 +1521,7 @@ void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
             }
         }
 
-        weightedSum
+        meshRefinement::weightedSum
         (
             mesh(),
             isMasterEdge,

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H

@@ -105,31 +105,6 @@ class medialAxisMeshMover
             const labelList&
         );
 
-        //- Weighted sum (over all subset of mesh points) by
-        //  summing contribution from (master) edges
-        template<class Type>
-        static void weightedSum
-        (
-            const polyMesh& mesh,
-            const PackedBoolList& isMasterEdge,
-            const labelList& meshEdges,
-            const labelList& meshPoints,
-            const edgeList& edges,
-            const scalarField& edgeWeights,
-            const Field<Type>& data,
-            Field<Type>& sum
-        );
-
-        void calculateEdgeWeights
-        (
-            const PackedBoolList& isMasterEdge,
-            const labelList& meshEdges,
-            const labelList& meshPoints,
-            const edgeList& edges,
-            scalarField& edgeWeights,
-            scalarField& invSumWeight
-        ) const;
-
 
         // Calculation of medial axis information
 
@@ -314,12 +289,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "medialAxisMeshMoverTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C

@@ -1,93 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "medialAxisMeshMover.H"
-#include "syncTools.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::medialAxisMeshMover::weightedSum
-(
-    const polyMesh& mesh,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    const scalarField& edgeWeights,
-    const Field<Type>& pointData,
-    Field<Type>& sum
-)
-{
-    if
-    (
-        mesh.nEdges() != isMasterEdge.size()
-     || edges.size() != meshEdges.size()
-     || edges.size() != edgeWeights.size()
-     || meshPoints.size() != pointData.size()
-    )
-    {
-        FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
-            << "Inconsistent sizes for edge or point data:"
-            << " meshEdges:" << meshEdges.size()
-            << " isMasterEdge:" << isMasterEdge.size()
-            << " edgeWeights:" << edgeWeights.size()
-            << " edges:" << edges.size()
-            << " pointData:" << pointData.size()
-            << " meshPoints:" << meshPoints.size()
-            << abort(FatalError);
-    }
-
-    sum.setSize(meshPoints.size());
-    sum = pTraits<Type>::zero;
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-
-            scalar eWeight = edgeWeights[edgeI];
-
-            label v0 = e[0];
-            label v1 = e[1];
-
-            sum[v0] += eWeight*pointData[v1];
-            sum[v1] += eWeight*pointData[v0];
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        meshPoints,
-        sum,
-        plusEqOp<Type>(),
-        pTraits<Type>::zero     // null value
-    );
-}
-
-
-// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -1949,6 +1949,64 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
 }
 
 
+void Foam::meshRefinement::calculateEdgeWeights
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    scalarField& edgeWeights,
+    scalarField& invSumWeight
+)
+{
+    const pointField& pts = mesh.points();
+
+    // Calculate edgeWeights and inverse sum of edge weights
+    edgeWeights.setSize(meshEdges.size());
+    invSumWeight.setSize(meshPoints.size());
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        scalar eMag = max
+        (
+            VSMALL,
+            mag
+            (
+                pts[meshPoints[e[1]]]
+              - pts[meshPoints[e[0]]]
+            )
+        );
+        edgeWeights[edgeI] = 1.0/eMag;
+    }
+
+    // Sum per point all edge weights
+    weightedSum
+    (
+        mesh,
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        scalarField(meshPoints.size(), 1.0),  // data
+        invSumWeight
+    );
+
+    // Inplace invert
+    forAll(invSumWeight, pointI)
+    {
+        scalar w = invSumWeight[pointI];
+
+        if (w > 0.0)
+        {
+            invSumWeight[pointI] = 1.0/w;
+        }
+    }
+}
+
+
 Foam::label Foam::meshRefinement::appendPatch
 (
     fvMesh& mesh,
@@ -2167,9 +2225,7 @@ Foam::label Foam::meshRefinement::addMeshedPatch
 //        );
 
         // Store
-        label sz = meshedPatches_.size();
+        meshedPatches_.append(name);
-        meshedPatches_.setSize(sz+1);
-        meshedPatches_[sz] = name;
 
         return patchI;
     }

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -729,6 +729,33 @@ public:
             //- Helper function: check that face zones are synced
             static void checkCoupledFaceZones(const polyMesh&);
 
+            //- Helper: calculate edge weights (1/length)
+            static void calculateEdgeWeights
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const labelList& meshPoints,
+                const edgeList& edges,
+                scalarField& edgeWeights,
+                scalarField& invSumWeight
+            );
+
+            //- Helper: weighted sum (over all subset of mesh points) by
+            //  summing contribution from (master) edges
+            template<class Type>
+            static void weightedSum
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const labelList& meshPoints,
+                const edgeList& edges,
+                const scalarField& edgeWeights,
+                const Field<Type>& data,
+                Field<Type>& sum
+            );
+
 
         // Refinement
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -458,10 +458,18 @@ void Foam::meshRefinement::markFeatureCellLevel
     // Database to pass into trackedParticle::move
     trackedParticle::trackingData td(startPointCloud, maxFeatureLevel);
 
+
     // Track all particles to their end position (= starting feature point)
     // Note that the particle might have started on a different processor
     // so this will transfer across nicely until we can start tracking proper.
     scalar maxTrackLen = 2.0*mesh_.bounds().mag();
+
+    if (debug&meshRefinement::FEATURESEEDS)
+    {
+        Pout<< "Tracking " << startPointCloud.size()
+            << " particles over distance " << maxTrackLen
+            << " to find the starting cell" << endl;
+    }
     startPointCloud.move(td, maxTrackLen);
 
 
@@ -485,6 +493,11 @@ void Foam::meshRefinement::markFeatureCellLevel
         IDLList<trackedParticle>()
     );
 
+    if (debug&meshRefinement::FEATURESEEDS)
+    {
+        Pout<< "Constructing cloud for cell marking" << endl;
+    }
+
     forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
     {
         const trackedParticle& startTp = iter();
@@ -532,11 +545,14 @@ void Foam::meshRefinement::markFeatureCellLevel
     while (true)
     {
         // Track all particles to their end position.
+        if (debug&meshRefinement::FEATURESEEDS)
+        {
+            Pout<< "Tracking " << cloud.size()
+                << " particles over distance " << maxTrackLen
+                << " to mark cells" << endl;
+        }
         cloud.move(td, maxTrackLen);
 
-
-        label nParticles = 0;
-
         // Make particle follow edge.
         forAllIter(Cloud<trackedParticle>, cloud, iter)
         {
@@ -578,15 +594,15 @@ void Foam::meshRefinement::markFeatureCellLevel
                 // seeded. Delete particle.
                 cloud.deleteParticle(tp);
             }
-            else
-            {
-                // Keep particle
-                nParticles++;
-            }
         }
 
-        reduce(nParticles, sumOp<label>());
+
-        if (nParticles == 0)
+        if (debug&meshRefinement::FEATURESEEDS)
+        {
+            Pout<< "Remaining particles " << cloud.size() << endl;
+        }
+
+        if (returnReduce(cloud.size(), sumOp<label>()) == 0)
         {
             break;
         }

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C

@@ -26,16 +26,12 @@ License
 #include "meshRefinement.H"
 #include "fvMesh.H"
 #include "globalIndex.H"
-
+#include "syncTools.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // Add a T entry
-template<class T> void meshRefinement::updateList
+template<class T> void Foam::meshRefinement::updateList
 (
     const labelList& newToOld,
     const T& nullValue,
@@ -59,7 +55,7 @@ template<class T> void meshRefinement::updateList
 
 
 template<class T>
-T meshRefinement::gAverage
+T Foam::meshRefinement::gAverage
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterElem,
@@ -109,7 +105,7 @@ T meshRefinement::gAverage
 
 
 template<class T>
-T meshRefinement::gAverage
+T Foam::meshRefinement::gAverage
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterElem,
@@ -162,7 +158,7 @@ T meshRefinement::gAverage
 
 // Compare two lists over all boundary faces
 template<class T>
-void meshRefinement::testSyncBoundaryFaceList
+void Foam::meshRefinement::testSyncBoundaryFaceList
 (
     const scalar tol,
     const string& msg,
@@ -222,7 +218,7 @@ void meshRefinement::testSyncBoundaryFaceList
 // Print list sorted by coordinates. Used for comparing non-parallel v.s.
 // parallel operation
 template<class T>
-void meshRefinement::collectAndPrint
+void Foam::meshRefinement::collectAndPrint
 (
     const UList<point>& points,
     const UList<T>& data
@@ -268,7 +264,11 @@ void meshRefinement::collectAndPrint
 
 //template<class T, class Mesh>
 template<class GeoField>
-void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
+void Foam::meshRefinement::addPatchFields
+(
+    fvMesh& mesh,
+    const word& patchFieldType
+)
 {
     HashTable<GeoField*> flds
     (
@@ -298,7 +298,11 @@ void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
 
 // Reorder patch field
 template<class GeoField>
-void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
+void Foam::meshRefinement::reorderPatchFields
+(
+    fvMesh& mesh,
+    const labelList& oldToNew
+)
 {
     HashTable<GeoField*> flds
     (
@@ -316,7 +320,11 @@ void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
 
 
 template<class Enum>
-int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
+int Foam::meshRefinement::readFlags
+(
+    const Enum& namedEnum,
+    const wordList& words
+)
 {
     int flags = 0;
 
@@ -330,8 +338,66 @@ int meshRefinement::readFlags(const Enum& namedEnum, const wordList& words)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class Type>
+void Foam::meshRefinement::weightedSum
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    const scalarField& edgeWeights,
+    const Field<Type>& pointData,
+    Field<Type>& sum
+)
+{
+    if
+    (
+        mesh.nEdges() != isMasterEdge.size()
+     || edges.size() != meshEdges.size()
+     || edges.size() != edgeWeights.size()
+     || meshPoints.size() != pointData.size()
+    )
+    {
+        FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
+            << "Inconsistent sizes for edge or point data:"
+            << " meshEdges:" << meshEdges.size()
+            << " isMasterEdge:" << isMasterEdge.size()
+            << " edgeWeights:" << edgeWeights.size()
+            << " edges:" << edges.size()
+            << " pointData:" << pointData.size()
+            << " meshPoints:" << meshPoints.size()
+            << abort(FatalError);
+    }
+
+    sum.setSize(meshPoints.size());
+    sum = pTraits<Type>::zero;
+
+    forAll(edges, edgeI)
+    {
+        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        {
+            const edge& e = edges[edgeI];
+
+            scalar eWeight = edgeWeights[edgeI];
+
+            label v0 = e[0];
+            label v1 = e[1];
+
+            sum[v0] += eWeight*pointData[v1];
+            sum[v1] += eWeight*pointData[v0];
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        meshPoints,
+        sum,
+        plusEqOp<Type>(),
+        pTraits<Type>::zero     // null value
+    );
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //