ENH: periodicAMI: added debug geometry writing

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -326,9 +326,30 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
 
         // So now we have agglomeration of the target side in
         // allRestrict:
-        //  0..size-1 : local agglomeration (= targetRestrictAddressing)
+        //  0..size-1 : local agglomeration (= targetRestrictAddressing
+        //              (but potentially permutated))
         //  size..    : agglomeration data from other processors
 
+
+        // The trickiness in this algorithm is finding out the compaction
+        // of the remote data (i.e. allocation of the coarse 'slots'). We could
+        // either send across the slot compaction maps or just make sure
+        // that we allocate the slots in exactly the same order on both sending
+        // and receiving side (e.g. if the submap is set up to send 4 items,
+        // the constructMap is also set up to receive 4 items.
+
+
+        // Short note about the various types of indices:
+        // - face indices : indices into the geometry.
+        // - coarse face indices : how the faces get agglomerated
+        // - transferred data : how mapDistribute sends/receives data,
+        // - slots : indices into data after distribution (e.g. stencil,
+        //           srcAddress/tgtAddress). Note: for fully local addressing
+        //           the slots are equal to face indices.
+        // A mapDistribute has:
+        // - a subMap : these are face indices
+        // - a constructMap : these are from 'transferred-date' to slots
+
         labelListList tgtSubMap(Pstream::nProcs());
 
         // Local subMap is just identity
@@ -340,10 +361,14 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         {
             if (procI != Pstream::myProcNo())
             {
-                // Combine entries that point to the same coarse element. All
+                // Combine entries that point to the same coarse element.
-                // the elements refer to local data so index into
+                // The important bit is to loop over the data (and hand out
-                // targetRestrictAddressing or allRestrict (since the same
+                // compact indices ) in 'transferred data' order. This
-                // for local data).
+                // guarantees that we're doing exactly the
+                // same on sending and receiving side - e.g. the fourth element
+                // in the subMap is the fourth element received in the
+                // constructMap
+
                 const labelList& elems = map.subMap()[procI];
                 const labelList& elemsMap =
                     map.constructMap()[Pstream::myProcNo()];
@@ -383,8 +408,13 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         labelList tgtCompactMap(map.constructSize());
 
         {
-            const labelList& elemsMap =
+            // Note that in special cases (e.g. 'appending' two AMIs) the
-                map.constructMap()[Pstream::myProcNo()];
+            // local size after distributing can be longer than the number
+            // of faces. I.e. it duplicates elements.
+            // Since we don't know this size instead we loop over all
+            // reachable elements (using the local constructMap)
+
+            const labelList& elemsMap = map.constructMap()[Pstream::myProcNo()];
             forAll(elemsMap, i)
             {
                 label fineElem = elemsMap[i];
@@ -503,7 +533,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         forAll(fineSrcAddress, faceI)
         {
             // All the elements contributing to faceI. Are slots in
-            // mapDistribute'd data.
+            // target data.
             const labelList& elems = fineSrcAddress[faceI];
             const scalarList& weights = fineSrcWeights[faceI];
             const scalar fineArea = fineSrcMagSf[faceI];

src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C

@@ -25,6 +25,8 @@ License
 
 #include "cyclicPeriodicAMIPolyPatch.H"
 #include "addToRunTimeSelectionTable.H"
+#include "OBJstream.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -263,6 +265,25 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
             nbrPoints
         );
 
+
+        autoPtr<OBJstream> thisStr;
+        autoPtr<OBJstream> nbrStr;
+        if (debug)
+        {
+            const fileName dir(boundaryMesh().mesh().time().timePath());
+
+            mkDir(dir);
+
+            thisStr.reset(new OBJstream(dir/name()+"_this.obj"));
+            nbrStr.reset(new OBJstream(dir/name()+"_nbr.obj"));
+            Info<< name() << " : dumping local duplicated geometry to "
+                << thisStr().name() << " and " << nbrStr().name() << endl;
+
+            thisStr().write(thisPatch0, thisPatch0.points());
+            nbrStr().write(nbrPatch0, nbrPatch0.points());
+        }
+
+
         // Construct a new AMI interpolation between the initial patch locations
         AMIPtr_.reset
         (
@@ -287,12 +308,16 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
         scalar srcSum(gAverage(AMIPtr_->srcWeightsSum()));
         scalar tgtSum(gAverage(AMIPtr_->tgtWeightsSum()));
 
-        // Direction of geometry replication
+        // Direction (or rather side of AMI : this or nbr patch) of 
+        // geometry replication
         bool direction = nTransforms_ >= 0;
 
         // Increase in the source weight sum for the last iteration in the
         // opposite direction. If the current increase is less than this, the
-        // direction is reversed.
+        // direction (= side of AMI to transform) is reversed.
+        // We switch the side to replicate instead of reversing the transform
+        // since at the coupledPolyPatch level there is no
+        // 'reverseTransformPosition' functionality.
         scalar srcSumDiff = 0;
 
         // Loop, replicating the geometry
@@ -311,6 +336,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
 
                 thisPatch.movePoints(thisPoints);
 
+                if (thisStr.valid())
+                {
+                    thisStr().write(thisPatch, thisPatch.points());
+                }
+
                 AMIPtr_->append(thisPatch, nbrPatch0);
             }
             else
@@ -319,6 +349,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
 
                 nbrPatch.movePoints(nbrPoints);
 
+                if (nbrStr.valid())
+                {
+                    nbrStr().write(nbrPatch, nbrPatch.points());
+                }
+
                 AMIPtr_->append(thisPatch0, nbrPatch);
             }
 
@@ -346,7 +381,7 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
         // Check that the match is complete
         if (iter == maxIter_)
         {
-            FatalErrorIn
+            SeriousErrorIn
             (
                 "void Foam::cyclicPeriodicAMIPolyPatch::resetPeriodicAMI"
                 "("
@@ -357,7 +392,12 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
                 << " do not couple to within a tolerance of "
                 << matchTolerance()
                 << " when transformed according to the periodic patch "
-                << periodicPatch.name() << "." << exit(FatalError);
+                << periodicPatch.name() << "." << nl
+                << "This is only acceptable during post-processing"
+                << "; not during running. Improve your mesh or increase"
+                << " the 'matchTolerance' setting in the patch specification."
+                << endl;
+                // << exit(FatalError);
         }
 
         // Check that both patches have replicated an integer number of times
@@ -367,7 +407,7 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
          || mag(tgtSum - floor(tgtSum + 0.5)) > tgtSum*matchTolerance()
         )
         {
-            FatalErrorIn
+            SeriousErrorIn
             (
                 "void Foam::cyclicPeriodicAMIPolyPatch::resetPeriodicAMI"
                 "("
@@ -377,7 +417,11 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
                 << "Patches " << name() << " and " << neighbPatch().name()
                 << " do not overlap an integer number of times when transformed"
                 << " according to the periodic patch "
-                << periodicPatch.name() << "." << exit(FatalError);
+                << periodicPatch.name() << "." << nl
+                << "This is only acceptable during post-processing"
+                << "; not during running. Improve your mesh or increase"
+                << " the 'matchTolerance' setting in the patch specification."
+                << endl;    //<< exit(FatalError);
         }
 
         // Normalise the weights