Initial implementation of the whole methodology

meshLibrary/utilities/smoothers/geometry/meshOptimizer/symmetryPlaneOptimisation/symmetryPlaneOptimisation.C

@@ -27,20 +27,10 @@ Description
 
 #include "demandDrivenData.H"
 #include "symmetryPlaneOptimisation.H"
-#include "partTetMesh.H"
-#include "meshSurfaceEngine.H"
-#include "meshSurfaceEngineModifier.H"
-
-#include "partTetMeshSimplex.H"
-#include "meshUntangler.H"
-#include "volumeOptimizer.H"
-#include "knuppMetric.H"
-
-#include <map>
-
-# ifdef USE_OMP
-#include <omp.h>
-# endif
+#include "polyMeshGenAddressing.H"
+#include "helperFunctions.H"
+#include "polyMeshGenChecks.H"
+#include "meshOptimizer.H"
 
 // #define DEBUGSearch
 
@@ -64,6 +54,7 @@ void symmetryPlaneOptimisation::detectSymmetryPlanes()
 
     forAll(boundaries, patchI)
     {
+        Info << "Patch " << patchI << " is of type " << boundaries[patchI].type() << endl;
         if( boundaries[patchI].type() == "symmetryPlane" )
         {
             std::pair<vector, label>& cs = centreSum[patchI];
@@ -108,14 +99,82 @@ void symmetryPlaneOptimisation::detectSymmetryPlanes()
         const std::pair<vector, label>& ns = normalSum[it->first];
         const point n = ns.first / ns.second;
 
-        const word pName = boundaries[it->first].patchName();
-        symmetryPlanes_.insert(std::make_pair(pName, plane(c, n)));
+        symmetryPlanes_.insert(std::make_pair(it->first, plane(c, n)));
     }
 }
 
-void symmetryPlaneOptimisation::pointInPlanes(VRWGraph&) const
+void symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
 {
+    const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
+    const pointFieldPMG& points = mesh_.points();
+    const faceListPMG& faces = mesh_.faces();
+
+    pointInPlanes.clear();
+    pointInPlanes.setSize(points.size());
+
+    forAll(boundaries, patchI)
+    {
+        if( boundaries[patchI].type() == "symmetryPlane" )
+        {
+            const label start = boundaries[patchI].patchStart();
+            const label end = start + boundaries[patchI].patchSize();
+            for(label faceI=start;faceI<end;++faceI)
+            {
+                const face& f = faces[faceI];
+
+                forAll(f, pI)
+                    pointInPlanes.appendIfNotIn(f[pI], patchI);
+            }
+        }
+    }
+
+    if( Pstream::parRun() )
+    {
+        const Map<label>& globalToLocal =
+            mesh_.addressingData().globalToLocalPointAddressing();
+        const VRWGraph& pointAtProcs = mesh_.addressingData().pointAtProcs();
+        const DynList<label>& neiProcs =
+            mesh_.addressingData().pointNeiProcs();
 
+        std::map<label, labelLongList> exchangeData;
+        forAll(neiProcs, i)
+            exchangeData[neiProcs[i]].clear();
+
+        forAllConstIter(Map<label>, globalToLocal, it)
+        {
+            const label pointI = it();
+
+            if( pointInPlanes.sizeOfRow(pointI) == 0 )
+                continue;
+
+            forAllRow(pointAtProcs, pointI, i)
+            {
+                const label neiProc = pointAtProcs(pointI, i);
+
+                if( neiProc == Pstream::myProcNo() )
+                    continue;
+
+                labelLongList& dataToSend = exchangeData[neiProc];
+
+                dataToSend.append(it.key());
+                dataToSend.append(pointInPlanes.sizeOfRow(pointI));
+                forAllRow(pointInPlanes, pointI, pipI)
+                    dataToSend.append(pointInPlanes(pointI, pipI));
+            }
+        }
+
+        labelLongList receivedData;
+        help::exchangeMap(exchangeData, receivedData);
+
+        for(label counter=0;counter<receivedData.size();)
+        {
+            const label pointI = globalToLocal[receivedData[counter++]];
+
+            const label size = receivedData[counter++];
+            for(label i=0;i<size;++i)
+                pointInPlanes.appendIfNotIn(pointI, receivedData[counter++]);
+        }
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -137,7 +196,55 @@ symmetryPlaneOptimisation::~symmetryPlaneOptimisation()
 
 void symmetryPlaneOptimisation::optimizeSymmetryPlanes()
 {
+    pointFieldPMG& points = mesh_.points();
 
+    Info << "Detected " << label(symmetryPlanes_.size()) << " symmetry planes in the mesh" << endl;
+
+    VRWGraph pointInPlane;
+    pointInPlanes(pointInPlane);
+
+    forAll(pointInPlane, pointI)
+    {
+        const label nPlanes = pointInPlane.sizeOfRow(pointI);
+
+        if( nPlanes > 3 )
+        {
+            WarningIn
+            (
+                "void symmetryPlaneOptimisation::optimizeSymmetryPlanes()"
+            ) << "Point " << pointI << " is in more than three symmetry"
+              << " planes. Cannot move it" << endl;
+        }
+
+        point& p = points[pointI];
+        vector disp(vector::zero);
+        for(label plI=0;plI<nPlanes;++plI)
+        {
+            //- point is in a plane
+            std::map<label, plane>::const_iterator it =
+                symmetryPlanes_.find(pointInPlane(pointI, 0));
+
+            const point newP = it->second.nearestPoint(points[pointI]);
+            disp += newP - p;
+        }
+
+        p += disp;
+    }
+
+    labelHashSet badFaces;
+    polyMeshGenChecks::checkFacePyramids(mesh_, false, VSMALL, &badFaces);
+
+    if( badFaces.size() )
+    {
+        WarningIn
+        (
+            "void symmetryPlaneOptimisation::optimizeSymmetryPlanes()"
+        ) << "Bad quality or inverted faces found in the mesh" << endl;
+
+        const label badFacesId = mesh_.addFaceSubset("invalidFaces");
+        forAllConstIter(labelHashSet, badFaces, it)
+            mesh_.addFaceToSubset(badFacesId, it.key());
+    }
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

meshLibrary/utilities/smoothers/geometry/meshOptimizer/symmetryPlaneOptimisation/symmetryPlaneOptimisation.H

@@ -61,7 +61,7 @@ class symmetryPlaneOptimisation
         polyMeshGen& mesh_;
 
         //- symmetry planes in the mesh
-        std::map<word, plane> symmetryPlanes_;
+        std::map<label, plane> symmetryPlanes_;
 
     // Private member functions
         //- detect symmetry planes

utilities/improveSymmetryPlanes/improveSymmetryPlanes.C

@@ -44,8 +44,10 @@ int main(int argc, char *argv[])
     polyMeshGen pmg(runTime);
     pmg.read();
 
+    Info << "Starting optimisation of symmetry planes" << endl;
     symmetryPlaneOptimisation(pmg).optimizeSymmetryPlanes();
 
+    Info << "Writing mesh" << endl;
     pmg.write();
 
     Info << "End\n" << endl;