From fac508c56d8ca347c60c03b7c08f38d9804c6a5c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 11 Nov 2011 13:47:26 +0000
Subject: [PATCH] ENH: cyclicAMI: consistent parallel transformations

---
 .../constraint/cyclic/cyclicPolyPatch.C       |  42 ++---
 .../constraint/cyclic/cyclicPolyPatch.H       |   5 +-
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.C   | 151 +++++++++++-------
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.H   |   9 +-
 4 files changed, 123 insertions(+), 84 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 5709f70b733..5fb4658db33 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -251,6 +251,7 @@ void Foam::cyclicPolyPatch::calcTransforms
         if (debug)
         {
             Pout<< "cyclicPolyPatch::calcTransforms :"
+                << " patch:" << name()
                 << " Max area error:" << 100*maxAreaDiff << "% at face:"
                 << maxAreaFacei << " at:" << half0Ctrs[maxAreaFacei]
                 << " coupled face at:" << half1Ctrs[maxAreaFacei]
@@ -264,17 +265,15 @@ void Foam::cyclicPolyPatch::calcTransforms
         {
             // Calculate using the given rotation axis and centre. Do not
             // use calculated normals.
-            label face0 = getConsistentRotationFace(half0Ctrs);
-            label face1 = face0;
-
-            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
-            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
+            vector n0 = findFaceMaxRadius(half0Ctrs);
+            vector n1 = findFaceMaxRadius(half1Ctrs);
             n0 /= mag(n0) + VSMALL;
             n1 /= mag(n1) + VSMALL;
 
             if (debug)
             {
                 Pout<< "cyclicPolyPatch::calcTransforms :"
+                    << " patch:" << name()
                     << " Specified rotation :"
                     << " n0:" << n0 << " n1:" << n1 << endl;
             }
@@ -330,6 +329,7 @@ void Foam::cyclicPolyPatch::calcTransforms
                 if (debug)
                 {
                     Pout<< "cyclicPolyPatch::calcTransforms :"
+                        << " patch:" << name()
                         << " Specified separation vector : "
                         << separationVector_ << endl;
                 }
@@ -426,17 +426,15 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
     {
         case ROTATIONAL:
         {
-            label face0 = getConsistentRotationFace(half0Ctrs);
-            label face1 = getConsistentRotationFace(half1Ctrs);
-
-            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
-            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
+            vector n0 = findFaceMaxRadius(half0Ctrs);
+            vector n1 = findFaceMaxRadius(half1Ctrs);
             n0 /= mag(n0) + VSMALL;
             n1 /= mag(n1) + VSMALL;
 
             if (debug)
             {
                 Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                    << " patch:" << name()
                     << " Specified rotation :"
                     << " n0:" << n0 << " n1:" << n1 << endl;
             }
@@ -485,6 +483,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
             if (debug)
             {
                 Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                    << " patch:" << name()
                     << "Specified translation : " << separationVector_
                     << endl;
             }
@@ -516,6 +515,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
                 if (debug)
                 {
                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                        << " patch:" << name()
                         << " Detected rotation :"
                         << " n0:" << n0 << " n1:" << n1 << endl;
                 }
@@ -548,6 +548,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
                 if (debug)
                 {
                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                        << " patch:" << name()
                         << " Detected translation :"
                         << " n0:" << n0 << " n1:" << n1
                         << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
@@ -566,30 +567,29 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
 }
 
 
-Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
+Foam::vector Foam::cyclicPolyPatch::findFaceMaxRadius
 (
     const pointField& faceCentres
 ) const
 {
     // Determine a face furthest away from the axis
 
-    const scalarField magRadSqr
-    (
-        magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
-    );
+    const vectorField n((faceCentres - rotationCentre_) ^ rotationAxis_);
+
+    const scalarField magRadSqr(magSqr(n));
 
-    label rotFace = findMax(magRadSqr);
+    label faceI = findMax(magRadSqr);
 
     if (debug)
     {
-        Info<< "getConsistentRotationFace(const pointField&)" << nl
-            << "    rotFace  = " << rotFace << nl
-            << "    point    =  " << faceCentres[rotFace] << nl
-            << "    distance = " << Foam::sqrt(magRadSqr[rotFace])
+        Info<< "findFaceMaxRadius(const pointField&) : patch: " << name() << nl
+            << "    rotFace  = " << faceI << nl
+            << "    point    = " << faceCentres[faceI] << nl
+            << "    distance = " << Foam::sqrt(magRadSqr[faceI])
             << endl;
     }
 
-    return rotFace;
+    return n[faceI];
 }
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
index b170f267edb..daffc183b3e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
@@ -130,9 +130,8 @@ class cyclicPolyPatch
                 scalarField& tols
             ) const;
 
-            //- For rotational cases, try to find a unique face on each side
-            //  of the cyclic.
-            label getConsistentRotationFace(const pointField&) const;
+            //- Return normal of face at max distance from rotation axis
+            vector findFaceMaxRadius(const pointField& faceCentres) const;
 
 
 protected:
diff --git a/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
index 5309057d0b4..b6769849f7d 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
@@ -30,6 +30,7 @@ License
 #include "Time.H"
 #include "addToRunTimeSelectionTable.H"
 #include "faceAreaIntersect.H"
+#include "ops.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,63 +45,68 @@ namespace Foam
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-Foam::label Foam::cyclicAMIPolyPatch::findFaceMaxRadius
+Foam::vector Foam::cyclicAMIPolyPatch::findFaceMaxRadius
 (
     const pointField& faceCentres
 ) const
 {
     // Determine a face furthest away from the axis
 
-    const scalarField magRadSqr
-    (
-        magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
-    );
+    const vectorField n((faceCentres - rotationCentre_) ^ rotationAxis_);
+
+    const scalarField magRadSqr(magSqr(n));
 
     label faceI = findMax(magRadSqr);
 
     if (debug)
     {
-        Info<< "findFaceMaxRadius(const pointField&)" << nl
+        Info<< "findFaceMaxRadius(const pointField&) : patch: " << name() << nl
             << "    rotFace  = " << faceI << nl
-            << "    point    =  " << faceCentres[faceI] << nl
+            << "    point    = " << faceCentres[faceI] << nl
             << "    distance = " << Foam::sqrt(magRadSqr[faceI])
             << endl;
     }
 
-    return faceI;
+    return n[faceI];
 }
 
 
-// * * * * * * * * * * * Protecetd Member Functions  * * * * * * * * * * * * //
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 void Foam::cyclicAMIPolyPatch::calcTransforms()
 {
-    if (size())
+    // Half0
+    const cyclicAMIPolyPatch& half0 = *this;
+    vectorField half0Areas(half0.size());
+    forAll(half0, facei)
     {
-        // Half0
-        const cyclicAMIPolyPatch& half0 = *this;
-        vectorField half0Areas(half0.size());
-        forAll(half0, facei)
-        {
-            half0Areas[facei] = half0[facei].normal(half0.points());
-        }
+        half0Areas[facei] = half0[facei].normal(half0.points());
+    }
 
-        // Half1
-        const cyclicAMIPolyPatch& half1 = neighbPatch();
-        vectorField half1Areas(half1.size());
-        forAll(half1, facei)
-        {
-            half1Areas[facei] = half1[facei].normal(half1.points());
-        }
+    // Half1
+    const cyclicAMIPolyPatch& half1 = neighbPatch();
+    vectorField half1Areas(half1.size());
+    forAll(half1, facei)
+    {
+        half1Areas[facei] = half1[facei].normal(half1.points());
+    }
 
-        calcTransforms
-        (
-            half0,
-            half0.faceCentres(),
-            half0Areas,
-            half1.faceCentres(),
-            half1Areas
-        );
+    calcTransforms
+    (
+        half0,
+        half0.faceCentres(),
+        half0Areas,
+        half1.faceCentres(),
+        half1Areas
+    );
+
+    if (debug)
+    {
+        Pout<< "calcTransforms() : patch: " << name() << nl
+            << "    forwardT = " << forwardT() << nl
+            << "    reverseT = " << reverseT() << nl
+            << "    separation = " << separation() << nl
+            << "    collocated = " << collocated() << nl << endl;
     }
 }
 
@@ -128,26 +134,30 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
 
     // Calculate transformation tensors
 
-    if ((half0Ctrs.size() <= 0) || (half1Ctrs.size() <= 0))
-    {
-        return;
-    }
-
     switch (transform_)
     {
         case ROTATIONAL:
         {
-            label face0 = findFaceMaxRadius(half0Ctrs);
-            label face1 = findFaceMaxRadius(half1Ctrs);
+            point n0 = vector::zero;
+            if (half0Ctrs.size())
+            {
+                n0 = findFaceMaxRadius(half0Ctrs);
+            }
+            point n1 = vector::zero;
+            if (half1Ctrs.size())
+            {
+                n1 = -findFaceMaxRadius(half1Ctrs);
+            }
+
+            reduce(n0, maxMagSqrOp<point>());
+            reduce(n1, maxMagSqrOp<point>());
 
-            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
-            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
             n0 /= mag(n0) + VSMALL;
             n1 /= mag(n1) + VSMALL;
 
             if (debug)
             {
-                Pout<< "cyclicAMIPolyPatch::calcTransforms :"
+                Pout<< "cyclicAMIPolyPatch::calcTransforms : patch:" << name()
                     << " Specified rotation :"
                     << " n0:" << n0 << " n1:" << n1 << endl;
             }
@@ -178,8 +188,8 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
         {
             if (debug)
             {
-                Pout<< "cyclicAMIPolyPatch::calcTransforms :"
-                    << "Specified translation : " << separationVector_
+                Pout<< "cyclicAMIPolyPatch::calcTransforms : patch:" << name()
+                    << " Specified translation : " << separationVector_
                     << endl;
             }
 
@@ -198,7 +208,8 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
         {
             if (debug)
             {
-                Pout<< "Assuming cyclic AMI pairs are colocated" << endl;
+                Pout<< " patch:" << name()
+                    << " Assuming cyclic AMI pairs are colocated" << endl;
             }
 
             const_cast<tensorField&>(forwardT()).clear();
@@ -232,7 +243,8 @@ void Foam::cyclicAMIPolyPatch::resetAMI() const
 
         if (debug)
         {
-            OFstream os(name() + "_neighbourPatch-org.obj");
+            const Time& t = boundaryMesh().mesh().time();
+            OFstream os(t.path()/name() + "_neighbourPatch-org.obj");
             meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
         }
 
@@ -250,10 +262,11 @@ void Foam::cyclicAMIPolyPatch::resetAMI() const
 
         if (debug)
         {
-            OFstream osN(name() + "_neighbourPatch-trans.obj");
+            const Time& t = boundaryMesh().mesh().time();
+            OFstream osN(t.path()/name() + "_neighbourPatch-trans.obj");
             meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
 
-            OFstream osO(name() + "_ownerPatch.obj");
+            OFstream osO(t.path()/name() + "_ownerPatch.obj");
             meshTools::writeOBJ(osO, this->localFaces(), localPoints());
         }
 
@@ -451,10 +464,10 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     coupledPolyPatch(pp, bm),
     nbrPatchName_(pp.nbrPatchName_),
     nbrPatchID_(-1),
-    transform_(UNKNOWN),
-    rotationAxis_(vector::zero),
-    rotationCentre_(point::zero),
-    separationVector_(vector::zero),
+    transform_(pp.transform_),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_),
     AMIPtr_(NULL),
     surfPtr_(NULL),
     surfDict_(dictionary::null)
@@ -477,10 +490,10 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     coupledPolyPatch(pp, bm, index, newSize, newStart),
     nbrPatchName_(nbrPatchName),
     nbrPatchID_(-1),
-    transform_(UNKNOWN),
-    rotationAxis_(vector::zero),
-    rotationCentre_(point::zero),
-    separationVector_(vector::zero),
+    transform_(pp.transform_),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_),
     AMIPtr_(NULL),
     surfPtr_(NULL),
     surfDict_(dictionary::null)
@@ -535,6 +548,23 @@ Foam::cyclicAMIPolyPatch::~cyclicAMIPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::cyclicAMIPolyPatch::coupled() const
+{
+    if
+    (
+        Pstream::parRun()
+     || (size() && neighbPatch().size())
+    )
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::label Foam::cyclicAMIPolyPatch::neighbPatchID() const
 {
     if (nbrPatchID_ == -1)
@@ -769,6 +799,7 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
     coupledPolyPatch::write(os);
     os.writeKeyword("neighbourPatch") << nbrPatchName_
         << token::END_STATEMENT << nl;
+
     switch (transform_)
     {
         case ROTATIONAL:
@@ -801,8 +832,12 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
         }
     }
 
-    os.writeKeyword(surfDict_.dictName());
-    os  << surfDict_;
+
+    if (surfDict_ != dictionary::null)
+    {
+        os.writeKeyword(surfDict_.dictName());
+        os  << surfDict_;
+    }
 }
 
 
diff --git a/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H b/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
index fa447ec6fe5..ced7968171c 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
+++ b/src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
@@ -97,8 +97,8 @@ private:
 
     // Private Member Functions
 
-        //- Return face ID of face at max distance from rotation axis
-        label findFaceMaxRadius(const pointField& faceCentres) const;
+        //- Return normal of face at max distance from rotation axis
+        vector findFaceMaxRadius(const pointField& faceCentres) const;
 
         void calcTransforms
         (
@@ -254,6 +254,11 @@ public:
 
         // Access
 
+            //- Return true only if is coupled. Note that for non-parallel
+            //  operation of a decomposed case this can return the wrong
+            //  result
+            virtual bool coupled() const;
+
             //- Neighbour patch name
             inline const word& neighbPatchName() const;
 
-- 
GitLab