BUG: cyclicPeriodicAMI: fixed map distribution on decomposed periodic AMI patches

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -345,6 +345,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
                 // targetRestrictAddressing or allRestrict (since the same
                 // for local data).
                 const labelList& elems = map.subMap()[procI];
+                const labelList& elemsMap =
+                    map.constructMap()[Pstream::myProcNo()];
                 labelList& newSubMap = tgtSubMap[procI];
                 newSubMap.setSize(elems.size());
 
@@ -353,7 +355,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
 
                 forAll(elems, i)
                 {
-                    label fineElem = elems[i];
+                    label fineElem = elemsMap[elems[i]];
                     label coarseElem = allRestrict[fineElem];
                     if (oldToNew[coarseElem] == -1)
                     {
@@ -371,16 +373,26 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
         // the sending map
 
         labelListList tgtConstructMap(Pstream::nProcs());
-        labelList tgtCompactMap;
 
         // Local constructMap is just identity
         {
             tgtConstructMap[Pstream::myProcNo()] =
                 identity(targetCoarseSize);
-
-            tgtCompactMap = targetRestrictAddressing;
         }
-        tgtCompactMap.setSize(map.constructSize());
+
+        labelList tgtCompactMap(map.constructSize());
+
+        {
+            const labelList& elemsMap =
+                map.constructMap()[Pstream::myProcNo()];
+            forAll(elemsMap, i)
+            {
+                label fineElem = elemsMap[i];
+                label coarseElem = allRestrict[fineElem];
+                tgtCompactMap[fineElem] = coarseElem;
+            }
+        }
+
         label compactI = targetCoarseSize;
 
         // Compact data from other processors
@@ -439,7 +451,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
             }
         }
 
-
         srcAddress.setSize(sourceCoarseSize);
         srcWeights.setSize(sourceCoarseSize);
 
@@ -1095,7 +1106,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
     const TargetPatch& tgtPatch
 )
 {
-    // create a new interpolation
+    // Create a new interpolation
     autoPtr<AMIInterpolation<SourcePatch, TargetPatch> > newPtr
     (
         new AMIInterpolation<SourcePatch, TargetPatch>
@@ -1110,7 +1121,145 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
         )
     );
 
-    // combine new and current data
+    // If parallel then combine the mapDistribution and re-index
+    if (singlePatchProc_ == -1)
+    {
+        labelListList& srcSubMap = srcMapPtr_->subMap();
+        labelListList& srcConstructMap = srcMapPtr_->constructMap();
+
+        labelListList& tgtSubMap = tgtMapPtr_->subMap();
+        labelListList& tgtConstructMap = tgtMapPtr_->constructMap();
+
+        labelListList& newSrcSubMap = newPtr->srcMapPtr_->subMap();
+        labelListList& newSrcConstructMap = newPtr->srcMapPtr_->constructMap();
+
+        labelListList& newTgtSubMap = newPtr->tgtMapPtr_->subMap();
+        labelListList& newTgtConstructMap = newPtr->tgtMapPtr_->constructMap();
+
+        // Re-calculate the source indices
+        {
+            labelList mapMap(0), newMapMap(0);
+            forAll(srcSubMap, procI)
+            {
+                mapMap.append
+                (
+                    identity(srcConstructMap[procI].size())
+                  + mapMap.size() + newMapMap.size()
+                );
+                newMapMap.append
+                (
+                    identity(newSrcConstructMap[procI].size())
+                  + mapMap.size() + newMapMap.size()
+                );
+            }
+
+            forAll(srcSubMap, procI)
+            {
+                forAll(srcConstructMap[procI], srcI)
+                {
+                    srcConstructMap[procI][srcI] =
+                        mapMap[srcConstructMap[procI][srcI]];
+                }
+            }
+
+            forAll(srcSubMap, procI)
+            {
+                forAll(newSrcConstructMap[procI], srcI)
+                {
+                    newSrcConstructMap[procI][srcI] =
+                        newMapMap[newSrcConstructMap[procI][srcI]];
+                }
+            }
+
+            forAll(tgtAddress_, tgtI)
+            {
+                forAll(tgtAddress_[tgtI], tgtJ)
+                {
+                    tgtAddress_[tgtI][tgtJ] =
+                        mapMap[tgtAddress_[tgtI][tgtJ]];
+                }
+            }
+
+            forAll(newPtr->tgtAddress_, tgtI)
+            {
+                forAll(newPtr->tgtAddress_[tgtI], tgtJ)
+                {
+                    newPtr->tgtAddress_[tgtI][tgtJ] =
+                        newMapMap[newPtr->tgtAddress_[tgtI][tgtJ]];
+                }
+            }
+        }
+
+        // Re-calculate the target indices
+        {
+            labelList mapMap(0), newMapMap(0);
+            forAll(srcSubMap, procI)
+            {
+                mapMap.append
+                (
+                    identity(tgtConstructMap[procI].size())
+                  + mapMap.size() + newMapMap.size()
+                );
+                newMapMap.append
+                (
+                    identity(newTgtConstructMap[procI].size())
+                  + mapMap.size() + newMapMap.size()
+                );
+            }
+
+            forAll(srcSubMap, procI)
+            {
+                forAll(tgtConstructMap[procI], tgtI)
+                {
+                    tgtConstructMap[procI][tgtI] =
+                        mapMap[tgtConstructMap[procI][tgtI]];
+                }
+            }
+
+            forAll(srcSubMap, procI)
+            {
+                forAll(newTgtConstructMap[procI], tgtI)
+                {
+                    newTgtConstructMap[procI][tgtI] =
+                        newMapMap[newTgtConstructMap[procI][tgtI]];
+                }
+            }
+
+            forAll(srcAddress_, srcI)
+            {
+                forAll(srcAddress_[srcI], srcJ)
+                {
+                    srcAddress_[srcI][srcJ] =
+                        mapMap[srcAddress_[srcI][srcJ]];
+                }
+            }
+
+            forAll(newPtr->srcAddress_, srcI)
+            {
+                forAll(newPtr->srcAddress_[srcI], srcJ)
+                {
+                    newPtr->srcAddress_[srcI][srcJ] =
+                        newMapMap[newPtr->srcAddress_[srcI][srcJ]];
+                }
+            }
+        }
+
+        // Sum the construction sizes
+        srcMapPtr_->constructSize() += newPtr->srcMapPtr_->constructSize();
+        tgtMapPtr_->constructSize() += newPtr->tgtMapPtr_->constructSize();
+
+        // Combine the maps
+        forAll(srcSubMap, procI)
+        {
+            srcSubMap[procI].append(newSrcSubMap[procI]);
+            srcConstructMap[procI].append(newSrcConstructMap[procI]);
+
+            tgtSubMap[procI].append(newTgtSubMap[procI]);
+            tgtConstructMap[procI].append(newTgtConstructMap[procI]);
+        }
+    }
+
+    // Combine new and current source data
     forAll(srcMagSf_, srcFaceI)
     {
         srcAddress_[srcFaceI].append(newPtr->srcAddress()[srcFaceI]);
@@ -1118,6 +1267,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
         srcWeightsSum_[srcFaceI] += newPtr->srcWeightsSum()[srcFaceI];
     }
 
+    // Combine new and current target data
     forAll(tgtMagSf_, tgtFaceI)
     {
         tgtAddress_[tgtFaceI].append(newPtr->tgtAddress()[tgtFaceI]);

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

@@ -44,6 +44,151 @@ namespace Foam
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+void Foam::cyclicPeriodicAMIPolyPatch::syncTransforms() const
+{
+    if (owner())
+    {
+        // At this point we guarantee that the transformations have been
+        // updated. There is one particular case, where if the periodic patch
+        // locally has zero faces then its transformation will not be set. We
+        // need to synchronise the transforms over the zero-sized patches as
+        // well.
+        //
+        // We can't put the logic into cyclicPolyPatch as
+        // processorCyclicPolyPatch uses cyclicPolyPatch functionality.
+        // Synchronisation will fail because processor-type patches do not exist
+        // on all processors.
+        //
+        // The use in cyclicPeriodicAMI is special; we use the patch even
+        // though we have no faces in it. Ideally, in future, the transformation
+        // logic will be abstracted out, and we won't need a periodic patch
+        // here. Until then, we can just fix the behaviour of the zero-sized
+        // coupled patches here
+
+        // Get the periodic patch
+        const coupledPolyPatch& periodicPatch
+        (
+            refCast<const coupledPolyPatch>
+            (
+                boundaryMesh()[periodicPatchID()]
+            )
+        );
+
+        // If there are any zero-sized periodic patches
+        if (returnReduce((size() && !periodicPatch.size()), orOp<bool>()))
+        {
+            if (periodicPatch.separation().size() > 1)
+            {
+                FatalErrorIn
+                (
+                    "cyclicPeriodicAMIPolyPatch::resetAMI"
+                    "(const AMIPatchToPatchInterpolation::interpolationMethod&"
+                    ") const"
+                )   << "Periodic patch " << periodicPatchName_
+                    << " has non-uniform separation vector "
+                    << periodicPatch.separation()
+                    << "This is not allowed inside " << type()
+                    << " patch " << name()
+                    << exit(FatalError);
+            }
+
+            if (periodicPatch.forwardT().size() > 1)
+            {
+                FatalErrorIn
+                (
+                    "cyclicPeriodicAMIPolyPatch::resetAMI"
+                    "(const AMIPatchToPatchInterpolation::interpolationMethod&"
+                    ") const"
+                )   << "Periodic patch " << periodicPatchName_
+                    << " has non-uniform transformation tensor "
+                    << periodicPatch.forwardT()
+                    << "This is not allowed inside " << type()
+                    << " patch " << name()
+                    << exit(FatalError);
+            }
+
+            // Note that zero-sized patches will have zero-sized fields for the
+            // separation vector, forward and reverse transforms. These need
+            // replacing with the transformations from other processors.
+
+            // Parallel in this context refers to a parallel transformation,
+            // rather than a rotational transformation
+            bool isParallel = periodicPatch.parallel();
+            reduce(isParallel, orOp<bool>());
+
+            if (isParallel)
+            {
+                // Sync a list of separation vectors
+                List<vectorField> sep(Pstream::nProcs());
+                sep[Pstream::myProcNo()] = periodicPatch.separation();
+                Pstream::gatherList(sep);
+                Pstream::scatterList(sep);
+
+                List<boolList> coll(Pstream::nProcs());
+                coll[Pstream::myProcNo()] = periodicPatch.collocated();
+                Pstream::gatherList(coll);
+                Pstream::scatterList(coll);
+
+                // If locally we have zero faces pick the first one that has a
+                // separation vector
+                if (!periodicPatch.size())
+                {
+                    forAll(sep, procI)
+                    {
+                        if (sep[procI].size())
+                        {
+                            const_cast<vectorField&>
+                            (
+                                periodicPatch.separation()
+                            ) = sep[procI];
+                            const_cast<boolList&>
+                            (
+                                periodicPatch.collocated()
+                            ) = coll[procI];
+
+                            break;
+                        }
+                    }
+                }
+            }
+            else
+            {
+                // Sync a list of forward and reverse transforms
+                List<tensorField> forwardT(Pstream::nProcs());
+                forwardT[Pstream::myProcNo()] = periodicPatch.forwardT();
+                Pstream::gatherList(forwardT);
+                Pstream::scatterList(forwardT);
+
+                List<tensorField> reverseT(Pstream::nProcs());
+                reverseT[Pstream::myProcNo()] = periodicPatch.reverseT();
+                Pstream::gatherList(reverseT);
+                Pstream::scatterList(reverseT);
+
+                // If locally we have zero faces pick the first one that has a
+                // transformation vector
+                if (!periodicPatch.size())
+                {
+                    forAll(forwardT, procI)
+                    {
+                        if (forwardT[procI].size())
+                        {
+                            const_cast<tensorField&>
+                            (
+                                periodicPatch.forwardT()
+                            ) = forwardT[procI];
+                            const_cast<tensorField&>
+                            (
+                                periodicPatch.reverseT()
+                            ) = reverseT[procI];
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
 (
     const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
@@ -60,6 +205,9 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
             )
         );
 
+        // Synchronise the transforms
+        syncTransforms();
+
         // Create copies of both patches' points, transformed to the owner
         pointField thisPoints0(localPoints());
         pointField nbrPoints0(neighbPatch().localPoints());

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

@@ -72,6 +72,9 @@ private:
 
     // Private Member Functions
 
+        //- Synchronise the periodic transformations
+        void syncTransforms() const;
+
         //- Reset the AMI interpolator
         virtual void resetAMI
         (