cyclic ordering fix

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -60,26 +60,34 @@ const NamedEnum<cyclicPolyPatch::transformType, 3>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::cyclicPolyPatch::calcTransforms()
+Foam::label Foam::cyclicPolyPatch::findMaxArea
+(
+    const pointField& points,
+    const faceList& faces
+)
 {
-    if (size() > 0)
-    {
-        const pointField& points = this->points();
+    label maxI = -1;
+    scalar maxAreaSqr = -GREAT;
 
-        maxI_ = -1;
-        scalar maxAreaSqr = -GREAT;
+    forAll(faces, faceI)
+    {
+        scalar areaSqr = magSqr(faces[faceI].normal(points));
 
-        for (label faceI = 0; faceI < size()/2; faceI++)
+        if (areaSqr > maxAreaSqr)
         {
-            const face& f = operator[](faceI);
-            scalar areaSqr = magSqr(f.normal(points));
-
-            if (areaSqr > maxAreaSqr)
-            {
-                maxAreaSqr = areaSqr;
-                maxI_ = faceI;
-            }
+            maxAreaSqr = areaSqr;
+            maxI = faceI;
         }
+    }
+    return maxI;
+}
+
+
+void Foam::cyclicPolyPatch::calcTransforms()
+{
+    if (size() > 0)
+    {
+        const pointField& points = this->points();
 
         primitivePatch half0
         (
@@ -369,6 +377,9 @@ bool Foam::cyclicPolyPatch::getGeometricHalves
 }
 
 
+// Given a split of faces into left and right half calculate the centres
+// and anchor points. Transform the left points so they align with the
+// right ones.
 void Foam::cyclicPolyPatch::getCentresAndAnchors
 (
     const primitivePatch& pp,
@@ -403,11 +414,17 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
         }
         default:
         {
-            // Assumes that cyclic is correctly ordered, so that face[maxI_]
-            // on each side is equivalent.
-            n0 = half0Faces[maxI_].normal(pp.points());
+            // Assumes that cyclic is planar. This is also the initial
+            // condition for patches without faces.
+
+            // Determine the face with max area on both halves. These
+            // two faces are used to determine the transformation tensors
+            label max0I = findMaxArea(pp.points(), half0Faces);
+            n0 = half0Faces[max0I].normal(pp.points());
             n0 /= mag(n0) + VSMALL;
-            n1 = half1Faces[maxI_].normal(pp.points());
+
+            label max1I = findMaxArea(pp.points(), half1Faces);
+            n1 = half1Faces[max1I].normal(pp.points());
             n1 /= mag(n1) + VSMALL;
         }
     }
@@ -596,7 +613,6 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     coupledPointsPtr_(NULL),
     coupledEdgesPtr_(NULL),
     featureCos_(0.9),
-    maxI_(-1),
     transform_(UNKNOWN),
     rotationAxis_(vector::zero),
     rotationCentre_(point::zero)
@@ -617,7 +633,6 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     coupledPointsPtr_(NULL),
     coupledEdgesPtr_(NULL),
     featureCos_(0.9),
-    maxI_(-1),
     transform_(UNKNOWN),
     rotationAxis_(vector::zero),
     rotationCentre_(point::zero)
@@ -658,7 +673,6 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     coupledPointsPtr_(NULL),
     coupledEdgesPtr_(NULL),
     featureCos_(pp.featureCos_),
-    maxI_(pp.maxI_),
     transform_(pp.transform_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_)
@@ -680,7 +694,6 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     coupledPointsPtr_(NULL),
     coupledEdgesPtr_(NULL),
     featureCos_(pp.featureCos_),
-    maxI_(pp.maxI_),
     transform_(pp.transform_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_)

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H

@@ -98,9 +98,6 @@ private:
         //  Used to split cyclic into halves.
         scalar featureCos_;
 
-        //- Index of largest cell face
-        label maxI_;
-
         //- Type of transformation - rotational or translational
         transformType transform_;
 
@@ -113,8 +110,12 @@ private:
 
     // Private member functions
 
+        //- Find amongst selected faces the one with the largest area
+        static label findMaxArea(const pointField&, const faceList&);
+
         void calcTransforms();
 
+
         // Face ordering
 
             //- Find the two parts of the faces of pp using feature edges.
@@ -169,16 +170,16 @@ protected:
     // Protected Member functions
 
         //- Initialise the calculation of the patch geometry
-        void initGeometry();
+        virtual void initGeometry();
 
         //- Calculate the patch geometry
-        void calcGeometry();
+        virtual void calcGeometry();
 
         //- Initialise the patches for moving points
-        void initMovePoints(const pointField&);
+        virtual void initMovePoints(const pointField&);
 
         //- Correct patches after moving points
-        void movePoints(const pointField&);
+        virtual void movePoints(const pointField&);
 
         //- Initialise the update of the patch topology
         virtual void initUpdateMesh();