ENH: GAMG: processor agglomeration extended for all interfaces

applications/solvers/basic/laplacianFoam/Make/files

@@ -1,3 +1,5 @@
+/* lduPrimitiveMesh.C */
+
 laplacianFoam.C
 
 EXE = $(FOAM_APPBIN)/laplacianFoam

applications/test/parallel-comm1/oneBlock_cyclic/0/T

@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    top
+    {
+        type            cyclic; //AMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    bottom
+    {
+        type            cyclic; //AMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    left
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+
+    right
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DT              4e-05;
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/blockMeshDict

@@ -0,0 +1,93 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    //- Left block
+    hex (0 1 2 3 4 5 6 7) (60 50 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    top
+    {
+        type            cyclic; //AMI;
+        neighbourPatch  bottom;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    bottom
+    {
+        type            cyclic; //AMI;
+        neighbourPatch  top;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    left
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    right
+    {
+        type patch;
+        faces
+        (
+            (1 2 5 6)
+        );
+    }
+
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/controlDict

@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs        (utilityFunctionObjects);
+
+DebugSwitches
+{
+//    mappedPatchBase 2;
+//    syncObjects     2;
+
+    disallowGenericFvPatchField 1;
+    disallowGenericFvsPatchField 1;
+    disallowGenericPointPatchField 1;
+    disallowGenericPolyPatch    1;
+    lduPrimitiveMesh            1;
+    GAMGAgglomeration           2;
+    GAMG                        2;
+}
+
+
+application     laplacianFoam;
+
+startFrom       startTime;  //latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          1;
+
+//writeControl    runTime;
+//writeInterval   0.1;
+writeControl    timeStep;
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+//functions
+//{
+//    syncObjects
+//    {
+//        type        syncObjects;
+//        libs        (utilityFunctionObjects);
+//
+//        // Where is data located relative to runTime. Given as a filename
+//        // with every '/' indicating a sub-objectRegistry w.r.t. runTime.
+//        // Local data is under <root>/send/processorXXX. After execution
+//        // data will be under the corresponding <root>/receive/processorYYY
+//        // objectRegistry
+//        //root        "level0/level1/level2";
+//    }
+//}
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/decomposeParDict

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- The total number of domains (mandatory)
+numberOfSubdomains  20;
+
+//- The decomposition method (mandatory)
+method          scotch;
+
+constraints
+{
+    patches
+    {
+        //- Keep owner and neighbour on same processor for faces in patches
+        //  (only makes sense for cyclic patches and cyclicAMI)
+        type    preservePatches;
+        patches (cyclic);
+    }
+}
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/fvSchemes

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;    //Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(T)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(DT,T) Gauss linear corrected;
+    laplacian(DTV,T) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/fvSolution

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    T
+    {
+        //solver          PCG;
+        //preconditioner  DIC;
+        solver          GAMG;
+        smoother        GaussSeidel;
+        nCellsInCoarsestLevel   1;
+        processorAgglomerator   masterCoarsest;
+
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/0/T

@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    AMI1
+    {
+        type            cyclicAMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    AMI2
+    {
+        type            cyclicAMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    bottom
+    {
+        type            zeroGradient;
+    }
+
+    left
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+
+    right
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/0/cellToRegion

@@ -0,0 +1,155 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2208                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       labelList;
+    location    "0";
+    object      cellToRegion;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+130
+(
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+)
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/0/processorID

@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2208                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0";
+    object      processorID;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    AMI1
+    {
+        type            cyclicAMI;
+        value           uniform 0;
+    }
+    AMI2
+    {
+        type            cyclicAMI;
+        value           uniform 0;
+    }
+    top
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    bottom
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    left
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    right
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DT              4e-05;
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/blockMeshDict

@@ -0,0 +1,136 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.1;
+
+vertices
+(
+    (0 0 0)
+    (0.5 0 0)
+    (0.5 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (0.5 0 0.1)
+    (0.5 1 0.1)
+    (0 1 0.1)
+
+
+    (0.5 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0.5 1 0)
+    (0.5 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0.5 1 0.1)
+
+);
+
+blocks
+(
+    //- Left block
+    hex (0 1 2 3 4 5 6 7) left (5 2 1) simpleGrading (1 1 1)
+
+    //- Right block
+    hex (8 9 10 11 12 13 14 15) right (5 2 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    AMI1
+    {
+        type            cyclicAMI;
+        neighbourPatch  AMI2;
+        transform       noOrdering;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+
+    AMI2
+    {
+        type            cyclicAMI;
+        neighbourPatch  AMI1;
+        transform       noOrdering;
+        faces
+        (
+            (8 12 15 11)
+        );
+    }
+
+    top
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+
+            (11 15 14 10)
+        );
+    }
+    bottom
+    {
+        type patch;
+        faces
+        (
+            (1 5 4 0)
+
+            (9 13 12 8)
+        );
+    }
+    left
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+
+        );
+    }
+
+    right
+    {
+        type patch;
+        faces
+        (
+            (10 14 13 9)
+        );
+    }
+
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+
+            (8 11 10 9)
+            (12 13 14 15)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/controlDict

@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs        (utilityFunctionObjects);
+
+DebugSwitches
+{
+//    mappedPatchBase 2;
+//    syncObjects     2;
+
+    disallowGenericFvPatchField 1;
+    disallowGenericFvsPatchField 1;
+    disallowGenericPointPatchField 1;
+    disallowGenericPolyPatch    1;
+    lduPrimitiveMesh            1;
+    GAMGAgglomeration           2;
+    GAMG                        2;
+}
+
+
+application     laplacianFoam;
+
+startFrom       startTime;  //latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          1;
+
+//writeControl    runTime;
+//writeInterval   0.1;
+writeControl    timeStep;
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+//functions
+//{
+//    syncObjects
+//    {
+//        type        syncObjects;
+//        libs        (utilityFunctionObjects);
+//
+//        // Where is data located relative to runTime. Given as a filename
+//        // with every '/' indicating a sub-objectRegistry w.r.t. runTime.
+//        // Local data is under <root>/send/processorXXX. After execution
+//        // data will be under the corresponding <root>/receive/processorYYY
+//        // objectRegistry
+//        //root        "level0/level1/level2";
+//    }
+//}
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/decomposeParDict

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- The total number of domains (mandatory)
+numberOfSubdomains  3;
+
+//- The decomposition method (mandatory)
+method          hierarchical;
+n               (3 1 1);
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/fvSchemes

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;    //Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(T)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(DT,T) Gauss linear corrected;
+    laplacian(DTV,T) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/fvSolution

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    T
+    {
+        //solver          PCG;
+        //preconditioner  DIC;
+        solver          GAMG;
+        smoother        GaussSeidel;
+        nCellsInCoarsestLevel   1;
+        processorAgglomerator   masterCoarsest;
+
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -409,12 +409,16 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
 
     if (debug & 2)
     {
+        const auto& coarseAddr = meshLevels_[fineLevelIndex].lduAddr();
+
         Pout<< "GAMGAgglomeration :"
             << " agglomerated level " << fineLevelIndex
             << " from nCells:" << fineMeshAddr.size()
             << " nFaces:" << upperAddr.size()
             << " to nCells:" << nCoarseCells
-            << " nFaces:" << nCoarseFaces
+            << " nFaces:" << nCoarseFaces << nl
+            << "    lower:" << flatOutput(coarseAddr.lowerAddr()) << nl
+            << "    upper:" << flatOutput(coarseAddr.upperAddr()) << nl
             << endl;
     }
 }
@@ -430,13 +434,29 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     const label levelIndex
 )
 {
-    const lduMesh& myMesh = meshLevels_[levelIndex-1];
+    // - Assemble all the procIDs in meshComm onto a single master
+    //   (procIDs[0]). This constructs a new communicator ('comm') first.
+    // - The master communicates with neighbouring masters using
+    //   allMeshComm
 
+    const lduMesh& myMesh = meshLevels_[levelIndex-1];
+    const label nOldInterfaces = myMesh.interfaces().size();
 
     procAgglomMap_.set(levelIndex, new labelList(procAgglomMap));
     agglomProcIDs_.set(levelIndex, new labelList(procIDs));
     procCommunicator_[levelIndex] = allMeshComm;
 
+    procAgglomCommunicator_.set
+    (
+        levelIndex,
+        new UPstream::communicator
+        (
+            meshComm,
+            procIDs
+        )
+    );
+    const label comm = agglomCommunicator(levelIndex);
+
     // These could only be set on the master procs but it is
     // quite convenient to also have them on the slaves
     procCellOffsets_.set(levelIndex, new labelList(0));
@@ -447,7 +467,7 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
 
     // Collect meshes
     PtrList<lduPrimitiveMesh> otherMeshes;
-    lduPrimitiveMesh::gather(meshComm, myMesh, procIDs, otherMeshes);
+    lduPrimitiveMesh::gather(comm, myMesh, otherMeshes);
 
     if (Pstream::myProcNo(meshComm) == procIDs[0])
     {
@@ -476,6 +496,41 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     }
 
 
+    // Scatter the procBoundaryMap back to the originating processor
+    // so it knows which proc boundaries are to be kept. This is used
+    // so we only send over interfaceFields on kept processors (see
+    // GAMGSolver::procAgglomerateMatrix)
+    // TBD: using sub-communicator here (instead of explicit procIDs). Should
+    //      use sub-communicators more in other places.
+    {
+        const CompactListList<label> data
+        (
+            CompactListList<label>::pack<labelList>
+            (
+                procBoundaryMap_[levelIndex]
+            )
+        );
+
+        // Make space
+        procBoundaryMap_[levelIndex].setSize(procIDs.size());
+        labelList& bMap = procBoundaryMap_[levelIndex][Pstream::myProcNo(comm)];
+        bMap.setSize(nOldInterfaces);
+
+        // Scatter relevant section to originating processor
+        UPstream::scatter
+        (
+            data.values().cdata(),
+            data.localSizes(),
+            data.offsets(),
+
+
+            bMap.data(),
+            bMap.size(),
+            comm
+        );
+    }
+
+
     // Combine restrict addressing
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -69,6 +69,7 @@ void Foam::GAMGAgglomeration::compactLevels
     {
         procAgglomMap_.setSize(nCreatedLevels);
         agglomProcIDs_.setSize(nCreatedLevels);
+        procAgglomCommunicator_.setSize(nCreatedLevels);
         procCellOffsets_.setSize(nCreatedLevels);
         procFaceMap_.setSize(nCreatedLevels);
         procBoundaryMap_.setSize(nCreatedLevels);
@@ -289,6 +290,7 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
     {
         procAgglomMap_.setSize(maxLevels_);
         agglomProcIDs_.setSize(maxLevels_);
+        procAgglomCommunicator_.setSize(maxLevels_);
         procCellOffsets_.setSize(maxLevels_);
         procFaceMap_.setSize(maxLevels_);
         procBoundaryMap_.setSize(maxLevels_);
@@ -569,6 +571,12 @@ Foam::label Foam::GAMGAgglomeration::procCommunicator(const label leveli) const
 }
 
 
+Foam::label Foam::GAMGAgglomeration::agglomCommunicator(const label leveli) const
+{
+    return procAgglomCommunicator_[leveli];
+}
+
+
 const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets
 (
     const label leveli

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H

@@ -131,6 +131,11 @@ protected:
             //- Communicator for given level
             mutable labelList procCommunicator_;
 
+            //- Communicator for collecting contributions. Note self-contained
+            //  communicator for ease of cleanup - this one gets allocated
+            //  internally, not by the GAMGProcAgglomerations
+            mutable PtrList<UPstream::communicator> procAgglomCommunicator_;
+
             //- Mapping from processor to procMeshLevel cells
             mutable PtrList<labelList> procCellOffsets_;
 
@@ -423,6 +428,17 @@ public:
                 const bool procAgglom
             ) const;
 
+            //- Prolong (interpolate by injection) cell field. Return reference
+            //  to (potentially scattered) coarse field
+            template<class Type>
+            const Field<Type>& prolongField
+            (
+                Field<Type>& ff,        // fine-level field
+                Field<Type>& work,      // work storage
+                const Field<Type>& cf,  // coarse-level field
+                const label coarseLevelIndex
+            ) const;
+
 
         // Processor agglomeration. Note that the mesh and agglomeration is
         // stored per fineLevel (even though it is the coarse level mesh that
@@ -466,6 +482,9 @@ public:
             //- Communicator for current level or -1
             label procCommunicator(const label fineLeveli) const;
 
+            //- Communicator for collecting contributions
+            label agglomCommunicator(const label fineLeveli) const;
+
             //- Mapping from processor to procMesh cells
             const labelList& cellOffsets(const label fineLeveli) const;
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerationTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -176,4 +177,56 @@ void Foam::GAMGAgglomeration::prolongField
 }
 
 
+template<class Type>
+const Foam::Field<Type>& Foam::GAMGAgglomeration::prolongField
+(
+    Field<Type>& ff,
+    Field<Type>& allCf,      // work storage
+    const Field<Type>& cf,
+    const label levelIndex
+) const
+{
+    const labelList& fineToCoarse = restrictAddressing_[levelIndex];
+
+    const label coarseLevelIndex = levelIndex+1;
+
+    if (hasProcMesh(coarseLevelIndex))
+    {
+        const label coarseComm =
+            UPstream::parent(procCommunicator_[coarseLevelIndex]);
+
+        const List<label>& procIDs = agglomProcIDs(coarseLevelIndex);
+        const labelList& offsets = cellOffsets(coarseLevelIndex);
+
+        const label localSize = nCells_[levelIndex];
+        allCf.resize_nocopy(localSize);
+
+        globalIndex::scatter
+        (
+            offsets,
+            coarseComm,
+            procIDs,
+            cf,
+            allCf,
+            UPstream::msgType(),
+            Pstream::commsTypes::nonBlocking    //Pstream::commsTypes::scheduled
+        );
+
+        forAll(fineToCoarse, i)
+        {
+            ff[i] = allCf[fineToCoarse[i]];
+        }
+        return allCf;
+    }
+    else
+    {
+        forAll(fineToCoarse, i)
+        {
+            ff[i] = cf[fineToCoarse[i]];
+        }
+        return cf;
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C

@@ -125,12 +125,12 @@ Foam::GAMGSolver::GAMGSolver
                     {
                         if (fineMeshInterfaces.set(intI))
                         {
-                            OStringStream os;
+                            OStringStream os(IOstreamOption::BINARY);
                             refCast<const GAMGInterface>
                             (
                                 fineMeshInterfaces[intI]
                             ).write(os);
-                            IStringStream is(os.str());
+                            IStringStream is(os.str(), IOstreamOption::BINARY);
 
                             dummyPrimMeshInterfaces.set
                             (
@@ -238,6 +238,8 @@ Foam::GAMGSolver::GAMGSolver
                     {
                         Pout<< "        " << i
                             << "\ttype:" << interfaces[i].type()
+                            << "\tsize:"
+                            << interfaces[i].interface().faceCells().size()
                             << endl;
                     }
                 }

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,7 +45,8 @@ Description
       - Coarse matrix scaling: performed by correction scaling, using steepest
         descent optimisation.
       - Type of cycle: V-cycle with optional pre-smoothing.
-      - Coarsest-level matrix solved using PCG or PBiCGStab.
+      - Coarsest-level matrix solved using any lduSolver (PCG, PBiCGStab,
+        smoothSolver) or direct solver on master processor
 
 SourceFiles
     GAMGSolver.C
@@ -190,9 +191,8 @@ class GAMGSolver
         //- Collect matrices from other processors
         void gatherMatrices
         (
-            const labelList& procIDs,
-            const lduMesh& dummyMesh,
-            const label meshComm,
+            const label destLevel,
+            const label comm,
 
             const lduMatrix& mat,
             const FieldField<Field, scalar>& interfaceBouCoeffs,
@@ -202,8 +202,7 @@ class GAMGSolver
             PtrList<lduMatrix>& otherMats,
             PtrList<FieldField<Field, scalar>>& otherBouCoeffs,
             PtrList<FieldField<Field, scalar>>& otherIntCoeffs,
-            List<boolList>& otherTransforms,
-            List<List<label>>& otherRanks
+            PtrList<PtrList<lduInterfaceField>>& otherInterfaces
         ) const;
 
         //- Agglomerate processor matrices

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C

@@ -283,74 +283,136 @@ void Foam::GAMGSolver::agglomerateInterfaceCoefficients
 
 void Foam::GAMGSolver::gatherMatrices
 (
-    const labelList& procIDs,
-    const lduMesh& dummyMesh,
-    const label meshComm,
+    const label destLevel,
+    const label comm,
 
+    // Local matrix
     const lduMatrix& mat,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const FieldField<Field, scalar>& interfaceIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
 
+    // Remote matrices
     PtrList<lduMatrix>& otherMats,
     PtrList<FieldField<Field, scalar>>& otherBouCoeffs,
     PtrList<FieldField<Field, scalar>>& otherIntCoeffs,
-    List<boolList>& otherTransforms,
-    List<List<label>>& otherRanks
+    PtrList<PtrList<lduInterfaceField>>& otherInterfaces
 ) const
 {
     if (debug & 2)
     {
+        const auto& procIDs = UPstream::procID(comm);
+
         Pout<< "GAMGSolver::gatherMatrices :"
             << " collecting matrices from procs:" << procIDs
-            << " using comm:" << meshComm << endl;
+            << " using comm:" << comm << endl;
     }
 
-    if (Pstream::myProcNo(meshComm) == procIDs[0])
+    const auto& boundaryMap = agglomeration_.boundaryMap(destLevel);
+
+    // Use PstreamBuffers
+    PstreamBuffers pBufs
+    (
+        UPstream::commsTypes::nonBlocking,
+        UPstream::msgType(),
+        comm
+    );
+
+    // Send to master
+    if (!UPstream::master(comm))
     {
-        // Master.
-        otherMats.setSize(procIDs.size()-1);
-        otherBouCoeffs.setSize(procIDs.size()-1);
-        otherIntCoeffs.setSize(procIDs.size()-1);
-        otherTransforms.setSize(procIDs.size()-1);
-        otherRanks.setSize(procIDs.size()-1);
+        // Mark valid interfaces
+        // -1   : not set
+        // >= 0 : coupled interface (might also be unmerged processor boundary)
+        //
+        // Note: most processor interfaces will disappear. Originally
+        // we did not know which ones were kept but this is now stored
+        // on the boundaryMap (even on the slave processors). So we can
+        // already filter here and avoid sending across typeNames etc.
 
-        for (label proci = 1; proci < procIDs.size(); proci++)
+        const label proci = UPstream::myProcNo(comm);
+
+        labelList validInterface(interfaces.size(), -1);
+        forAll(interfaces, intI)
         {
-            label otherI = proci-1;
+            const label allIntI = boundaryMap[proci][intI];
+            if (interfaces.set(intI) && allIntI != -1)
+            {
+                validInterface[intI] = intI;
+            }
+        }
 
-            IPstream fromProc
-            (
-                Pstream::commsTypes::scheduled,
-                procIDs[proci],
-                0,          // bufSize
-                Pstream::msgType(),
-                meshComm
-            );
+        UOPstream toMaster(UPstream::masterNo(), pBufs);
 
-            otherMats.set(otherI, new lduMatrix(dummyMesh, fromProc));
+        toMaster<< mat << token::SPACE << validInterface;
+
+        forAll(validInterface, intI)
+        {
+            if (validInterface[intI] != -1)
+            {
+                const auto& interface = refCast<const GAMGInterfaceField>
+                (
+                    interfaces[intI]
+                );
+
+                toMaster
+                    << interfaceBouCoeffs[intI]
+                    << interfaceIntCoeffs[intI]
+                    << interface.type();
+                interface.write(toMaster);
+            }
+        }
+    }
+
+    // Wait for finish
+    pBufs.finishedGathers();
+
+    // Consume
+    if (UPstream::master(comm))
+    {
+        const label nProcs = UPstream::nProcs(comm);
+
+        const lduMesh& destMesh = agglomeration_.meshLevel(destLevel);
+        lduInterfacePtrsList destInterfaces = destMesh.interfaces();
+
+        // Master.
+        otherMats.setSize(nProcs-1);
+        otherBouCoeffs.setSize(nProcs-1);
+        otherIntCoeffs.setSize(nProcs-1);
+        otherInterfaces.setSize(nProcs-1);
+
+        for (const int proci : UPstream::subProcs(comm))
+        {
+            const label otherI = proci-1;
+
+            UIPstream fromProc(proci, pBufs);
+
+            otherMats.set(otherI, new lduMatrix(destMesh, fromProc));
 
             // Receive number of/valid interfaces
-            boolList& procTransforms = otherTransforms[otherI];
-            List<label>& procRanks = otherRanks[otherI];
+            // >= 0 : remote interface index
+            // -1   : invalid interface
+            const labelList validInterface(fromProc);
 
-            fromProc >> procTransforms;
-            fromProc >> procRanks;
-
-            // Size coefficients
             otherBouCoeffs.set
             (
                 otherI,
-                new FieldField<Field, scalar>(procRanks.size())
+                new FieldField<Field, scalar>(validInterface.size())
             );
             otherIntCoeffs.set
             (
                 otherI,
-                new FieldField<Field, scalar>(procRanks.size())
+                new FieldField<Field, scalar>(validInterface.size())
             );
-            forAll(procRanks, intI)
+            otherInterfaces.set
+            (
+                otherI,
+                new PtrList<lduInterfaceField>(validInterface.size())
+            );
+
+            forAll(validInterface, intI)
             {
-                if (procRanks[intI] != -1)
+                if (validInterface[intI] != -1)
                 {
                     otherBouCoeffs[otherI].set
                     (
@@ -362,49 +424,25 @@ void Foam::GAMGSolver::gatherMatrices
                         intI,
                         new scalarField(fromProc)
                     );
-                }
-            }
-        }
-    }
-    else
-    {
-        // Send to master
 
-        // Count valid interfaces
-        boolList procTransforms(interfaceBouCoeffs.size(), false);
-        List<label> procRanks(interfaceBouCoeffs.size(), -1);
-        forAll(interfaces, intI)
-        {
-            if (interfaces.set(intI))
-            {
-                const processorLduInterfaceField& interface =
-                    refCast<const processorLduInterfaceField>
+                    const word coupleType(fromProc);
+
+                    const label allIntI = boundaryMap[proci][intI];
+
+                    otherInterfaces[otherI].set
                     (
-                        interfaces[intI]
+                        intI,
+                        GAMGInterfaceField::New
+                        (
+                            coupleType,
+                            refCast<const GAMGInterface>
+                            (
+                                destInterfaces[allIntI]
+                            ),
+                            fromProc
+                        ).release()
                     );
-
-                procTransforms[intI] = interface.doTransform();
-                procRanks[intI] = interface.rank();
-            }
-        }
-
-        OPstream toMaster
-        (
-            Pstream::commsTypes::scheduled,
-            procIDs[0],
-            0,
-            Pstream::msgType(),
-            meshComm
-        );
-
-        toMaster << mat << procTransforms << procRanks;
-        forAll(procRanks, intI)
-        {
-            if (procRanks[intI] != -1)
-            {
-                toMaster
-                    << interfaceBouCoeffs[intI]
-                    << interfaceIntCoeffs[intI];
+                }
             }
         }
     }
@@ -434,9 +472,11 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         interfaceLevelsBouCoeffs_[levelI];
     const FieldField<Field, scalar>& coarsestIntCoeffs =
         interfaceLevelsIntCoeffs_[levelI];
-    const lduMesh& coarsestMesh = coarsestMatrix.mesh();
 
-    label coarseComm = coarsestMesh.comm();
+    // Communicator containing all processors to combine (=agglomProcIDs).
+    // Result will be on master of communicator.
+    const label agglomComm = agglomeration_.agglomCommunicator(levelI+1);
+
 
     // Gather all matrix coefficients onto agglomProcIDs[0]
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -444,43 +484,29 @@ void Foam::GAMGSolver::procAgglomerateMatrix
     PtrList<lduMatrix> otherMats;
     PtrList<FieldField<Field, scalar>> otherBouCoeffs;
     PtrList<FieldField<Field, scalar>> otherIntCoeffs;
-    List<boolList> otherTransforms;
-    List<List<label>> otherRanks;
+    PtrList<PtrList<lduInterfaceField>> otherInterfaces;
     gatherMatrices
     (
-        agglomProcIDs,
-        coarsestMesh,
-        coarseComm,
+        levelI+1,                       // allMesh level (only on master)
+        agglomComm,
 
-        coarsestMatrix,
-        coarsestBouCoeffs,
-        coarsestIntCoeffs,
-        coarsestInterfaces,
+        coarsestMatrix,                 // master before gathering
+        coarsestBouCoeffs,              // master before gathering
+        coarsestIntCoeffs,              // master before gathering
+        coarsestInterfaces,             // master before gathering
 
-        otherMats,
-        otherBouCoeffs,
-        otherIntCoeffs,
-        otherTransforms,
-        otherRanks
+        otherMats,                      // slave info
+        otherBouCoeffs,                 // slave info
+        otherIntCoeffs,                 // slave info
+        otherInterfaces                 // slave info
     );
 
 
-    if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
+    if (UPstream::master(agglomComm))
     {
         // Agglomerate all matrix
         // ~~~~~~~~~~~~~~~~~~~~~~
 
-        //Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
-        //
-        //forAll(otherMats, i)
-        //{
-        //    Pout<< "** otherMats " << i << " "
-        //        << otherMats[i].info()
-        //        << endl;
-        //}
-        //Pout<< endl;
-
-
         const lduMesh& allMesh = agglomeration_.meshLevel(levelI+1);
         const labelList& cellOffsets = agglomeration_.cellOffsets(levelI+1);
         const labelListList& faceMap = agglomeration_.faceMap(levelI+1);
@@ -566,7 +592,8 @@ void Foam::GAMGSolver::procAgglomerateMatrix
             allInterfaceIntCoeffs.set(intI, new scalarField(size));
         }
 
-        labelList nBounFaces(allMeshInterfaces.size());
+        UPtrList<lduInterfaceField> otherFlds(0);
+
         forAll(boundaryMap, proci)
         {
             const FieldField<Field, scalar>& procBouCoeffs
@@ -582,6 +609,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
               : otherIntCoeffs[proci-1]
             );
 
+
             const labelList& bMap = boundaryMap[proci];
             forAll(bMap, procIntI)
             {
@@ -594,47 +622,76 @@ void Foam::GAMGSolver::procAgglomerateMatrix
 
                     if (!allInterfaces.set(allIntI))
                     {
-                        // Construct lduInterfaceField
+                        const GAMGInterface& intf = refCast<const GAMGInterface>
+                        (
+                            allMeshInterfaces[allIntI]
+                        );
 
-                        bool doTransform = false;
-                        int rank = -1;
                         if (proci == 0)
                         {
-                            const processorGAMGInterfaceField& procInt =
-                                refCast
-                                <
-                                    const processorGAMGInterfaceField
-                                >
+                            // Clone my local interfaceField. Since it is from
+                            // this processor it will already exist, even if it
+                            // is a processor one.
+
+                            const auto& ffld =
+                            refCast<const GAMGInterfaceField>
+                            (
+                                coarsestInterfaces[procIntI]
+                            );
+
+                            allPrimitiveInterfaces.set
+                            (
+                                allIntI,
+                                ffld.clone
                                 (
-                                    coarsestInterfaces[procIntI]
-                                );
-                            doTransform = procInt.doTransform();
-                            rank = procInt.rank();
+                                    intf,
+                                    otherFlds
+                                ).release()
+                            );
                         }
                         else
                         {
-                            doTransform =
-                                otherTransforms[proci-1][procIntI];
-                            rank = otherRanks[proci-1][procIntI];
-                        }
-
-                        allPrimitiveInterfaces.set
-                        (
-                            allIntI,
-                            GAMGInterfaceField::New
-                            (
-                                refCast<const GAMGInterface>
+                            // Recreate a remote interfaceField
+                            if (otherInterfaces[proci-1].set(procIntI))
+                            {
+                                const auto& ffld =
+                                refCast<const GAMGInterfaceField>
                                 (
-                                    allMeshInterfaces[allIntI]
-                                ),
-                                doTransform,
-                                rank
-                            ).ptr()
-                        );
+                                    otherInterfaces[proci-1][procIntI]
+                                );
+
+                                allPrimitiveInterfaces.set
+                                (
+                                    allIntI,
+                                    ffld.clone
+                                    (
+                                        intf,
+                                        otherFlds
+                                    ).release()
+                                );
+                            }
+                            else
+                            {
+                                // Recreate a default interfaceField with
+                                // sensible defaults.
+                                // Should not occur since all unmerged
+                                // processor interfaces get transferred.
+                                allPrimitiveInterfaces.set
+                                (
+                                    allIntI,
+                                    GAMGInterfaceField::New
+                                    (
+                                        intf,
+                                        false,  //doTransform,
+                                        0       //rank
+                                    ).ptr()
+                                );
+                            }
+                        }
                         allInterfaces.set
                         (
                             allIntI,
-                            &allPrimitiveInterfaces[allIntI]
+                            allPrimitiveInterfaces.get(allIntI)
                         );
                     }
 
@@ -669,7 +726,6 @@ void Foam::GAMGSolver::procAgglomerateMatrix
                     const scalarField& procBou = procBouCoeffs[procIntI];
                     const scalarField& procInt = procIntCoeffs[procIntI];
 
-
                     forAll(map, i)
                     {
                         if (map[i] >= 0)
@@ -703,6 +759,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
             }
         }
 
+
         //Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl;
         //
         //forAll(allInterfaces, intI)
@@ -712,19 +769,18 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         //        Pout<< "    patch:" << intI
         //            << " type:" << allInterfaces[intI].type()
         //            << " size:"
-        //            << allInterfaces[intI].interface().
-        //                faceCells().size()
+        //            << allInterfaces[intI].interface().faceCells().size()
         //            << endl;
         //
-        //        //const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
-        //        //const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
-        //        //forAll(bouCoeffs, facei)
-        //        //{
-        //        //    Pout<< "        " << facei
-        //        //        << "\tbou:" << bouCoeffs[facei]
-        //        //        << "\tint:" << intCoeffs[facei]
-        //        //        << endl;
-        //        //}
+        //        const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
+        //        const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
+        //        forAll(bouCoeffs, facei)
+        //        {
+        //            Pout<< "        " << facei
+        //                << "\tbou:" << bouCoeffs[facei]
+        //                << "\tint:" << intCoeffs[facei]
+        //                << endl;
+        //        }
         //    }
         //}
     }

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -8,6 +8,7 @@
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2016-2021,2023 OpenCFD Ltd.
     Copyright (C) 2023 Huawei (Yu Ankun)
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -301,6 +302,9 @@ void Foam::GAMGSolver::Vcycle
 
     solveScalarField dummyField(0);
 
+    // Work storage for prolongation
+    solveScalarField work;
+
     for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
     {
         if (coarseCorrFields.set(leveli))
@@ -321,16 +325,18 @@ void Foam::GAMGSolver::Vcycle
                 preSmoothedCoarseCorrField = coarseCorrFields[leveli];
             }
 
-            agglomeration_.prolongField
+
+            // Prolong correction to leveli
+            const auto& cf = agglomeration_.prolongField
             (
-                coarseCorrFields[leveli],
+                coarseCorrFields[leveli],   // current level
+                work,
                 (
                     coarseCorrFields.set(leveli + 1)
                   ? coarseCorrFields[leveli + 1]
                   : dummyField              // dummy value
                 ),
-                leveli + 1,
-                true
+                leveli + 1
             );
 
 
@@ -346,43 +352,21 @@ void Foam::GAMGSolver::Vcycle
                     ACf.operator const solveScalarField&()
                 );
 
-            if (interpolateCorrection_) //&& leveli < coarsestLevel - 2)
+            if (interpolateCorrection_)
             {
-                if
+                // Normal operation : have both coarse level and fine
+                // level. No processor agglomeration
+                interpolate
                 (
-                    coarseCorrFields.set(leveli+1)
-                 && (
-                        matrixLevels_[leveli].mesh().comm()
-                     == matrixLevels_[leveli+1].mesh().comm()
-                    )
-                )
-                {
-                    // Normal operation : have both coarse level and fine
-                    // level. No processor agglomeration
-                    interpolate
-                    (
-                        coarseCorrFields[leveli],
-                        ACfRef,
-                        matrixLevels_[leveli],
-                        interfaceLevelsBouCoeffs_[leveli],
-                        interfaceLevels_[leveli],
-                        agglomeration_.restrictAddressing(leveli + 1),
-                        coarseCorrFields[leveli + 1],
-                        cmpt
-                    );
-                }
-                else
-                {
-                    interpolate
-                    (
-                        coarseCorrFields[leveli],
-                        ACfRef,
-                        matrixLevels_[leveli],
-                        interfaceLevelsBouCoeffs_[leveli],
-                        interfaceLevels_[leveli],
-                        cmpt
-                    );
-                }
+                    coarseCorrFields[leveli],
+                    ACfRef,
+                    matrixLevels_[leveli],
+                    interfaceLevelsBouCoeffs_[leveli],
+                    interfaceLevels_[leveli],
+                    agglomeration_.restrictAddressing(leveli + 1),
+                    cf,
+                    cmpt
+                );
             }
 
             // Scale coarse-grid correction field

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceField.C

@@ -35,6 +35,7 @@ namespace Foam
     defineTypeNameAndDebug(GAMGInterfaceField, 0);
     defineRunTimeSelectionTable(GAMGInterfaceField, lduInterface);
     defineRunTimeSelectionTable(GAMGInterfaceField, lduInterfaceField);
+    defineRunTimeSelectionTable(GAMGInterfaceField, Istream);
 }
 
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,6 +104,18 @@ public:
             (GAMGCp, doTransform, rank)
         );
 
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            GAMGInterfaceField,
+            Istream,
+            (
+                const GAMGInterface& GAMGCp,
+                Istream& is
+            ),
+            (GAMGCp, is)
+        );
+
 
     // Selectors
 
@@ -124,6 +136,15 @@ public:
             const int rank
         );
 
+        //- Return a pointer to a new interface created on freestore given
+        //  the fine interface and stream
+        static autoPtr<GAMGInterfaceField> New
+        (
+            const word& patchFieldType,
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
 
     // Constructors
 
@@ -150,6 +171,28 @@ public:
             interface_(GAMGCp)
         {}
 
+        //- Construct from GAMG interface and fine level interface field
+        GAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        )
+        :
+            lduInterfaceField(GAMGCp),
+            interface_(GAMGCp)
+        {}
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other // other
+        ) const = 0;
+
+
+    //- Destructor
+    virtual ~GAMGInterfaceField() = default;
+
 
     // Member Functions
 
@@ -161,6 +204,11 @@ public:
                 return interface_;
             }
 
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const = 0;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceFieldNew.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2021,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -79,4 +79,27 @@ Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
 }
 
 
+Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
+(
+    const word& patchFieldType,
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+{
+    auto* ctorPtr = IstreamConstructorTable(patchFieldType);
+
+    if (!ctorPtr)
+    {
+        FatalErrorInLookup
+        (
+            "GAMGInterfaceField",
+            patchFieldType,
+            *IstreamConstructorTablePtr_
+        ) << exit(FatalError);
+    }
+
+    return autoPtr<GAMGInterfaceField>(ctorPtr(GAMGCp, is));
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,6 +47,12 @@ namespace Foam
         cyclicGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        cyclicGAMGInterfaceField,
+        Istream
+    );
 
     // Add under name cyclicSlip
     addNamedToRunTimeSelectionTable
@@ -63,6 +69,13 @@ namespace Foam
         lduInterfaceField,
         cyclicSlip
     );
+    addNamedToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        cyclicGAMGInterfaceField,
+        Istream,
+        cyclicSlip
+    );
 }
 
 
@@ -101,6 +114,38 @@ Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
 {}
 
 
+Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    cyclicInterface_(refCast<const cyclicGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
+Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const lduInterfaceField& local,
+    const UPtrList<lduInterfaceField>& other
+)
+:
+    GAMGInterfaceField(GAMGCp, local),
+    cyclicInterface_(refCast<const cyclicGAMGInterface>(GAMGCp)),
+    doTransform_(false),
+    rank_(0)
+{
+    const auto& p = refCast<const cyclicLduInterfaceField>(local);
+
+    doTransform_ = p.doTransform();
+    rank_ = p.rank();
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
@@ -133,4 +178,12 @@ void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
 }
 
 
+void Foam::cyclicGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.H

@@ -100,6 +100,39 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        cyclicGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct from GAMG interface and local and remote fields
+        cyclicGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const lduInterfaceField& local, // local
+            const UPtrList<lduInterfaceField>& other // other
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other // other
+        ) const
+        {
+            return autoPtr<GAMGInterfaceField>
+            (
+                new cyclicGAMGInterfaceField
+                (
+                    GAMGCp,
+                    *this,      // local field
+                    other       // other fields
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicGAMGInterfaceField() = default;
@@ -157,6 +190,12 @@ public:
                 const direction cmpt,
                 const Pstream::commsTypes commsType
             ) const;
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorCyclicGAMGInterfaceField/processorCyclicGAMGInterfaceField.C

@@ -46,6 +46,12 @@ namespace Foam
         processorCyclicGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        processorCyclicGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -72,6 +78,16 @@ Foam::processorCyclicGAMGInterfaceField::processorCyclicGAMGInterfaceField
 {}
 
 
+Foam::processorCyclicGAMGInterfaceField::processorCyclicGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    processorGAMGInterfaceField(GAMGCp, is)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::processorCyclicGAMGInterfaceField::~processorCyclicGAMGInterfaceField()

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorCyclicGAMGInterfaceField/processorCyclicGAMGInterfaceField.H

@@ -87,6 +87,24 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        processorCyclicGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            NotImplemented;
+            return autoPtr<GAMGInterfaceField>(nullptr);
+        }
+
 
     // Destructor
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C

@@ -47,6 +47,12 @@ namespace Foam
         processorGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        processorGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -88,6 +94,19 @@ Foam::processorGAMGInterfaceField::processorGAMGInterfaceField
 {}
 
 
+Foam::processorGAMGInterfaceField::processorGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    procInterface_(refCast<const processorGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
@@ -193,4 +212,12 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
 }
 
 
+void Foam::processorGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H

@@ -119,6 +119,31 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        processorGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& patch,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            return autoPtr<GAMGInterfaceField>
+            (
+                new processorGAMGInterfaceField
+                (
+                    patch,
+                    doTransform_,
+                    rank_
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~processorGAMGInterfaceField() = default;
@@ -201,6 +226,12 @@ public:
             {
                 return rank_;
             }
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H

@@ -197,6 +197,26 @@ public:
             Istream& is
         );
 
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,  // current to fine interface
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            NotImplemented;
+            return nullptr;
+        }
+
+
 
     // Member Functions
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -162,6 +162,38 @@ Foam::cyclicGAMGInterface::cyclicGAMGInterface
 {}
 
 
+Foam::cyclicGAMGInterface::cyclicGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    const lduInterface& fineInterface,
+    const labelList& interfaceMap,
+    const labelUList& faceCells,
+    const labelUList& faceRestrictAddresssing,
+    const labelUList& faceOffsets,
+    const lduInterfacePtrsList& allInterfaces
+)
+:
+    GAMGInterface
+    (
+        index,
+        coarseInterfaces,
+        faceCells,
+        faceRestrictAddresssing
+    ),
+    neighbPatchID_
+    (
+        interfaceMap.find
+        (
+            refCast<const cyclicLduInterface>(fineInterface).neighbPatchID()
+        )
+    ),
+    owner_(refCast<const cyclicLduInterface>(fineInterface).owner()),
+    forwardT_(refCast<const cyclicLduInterface>(fineInterface).forwardT()),
+    reverseT_(refCast<const cyclicLduInterface>(fineInterface).reverseT())
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::labelField> Foam::cyclicGAMGInterface::internalFieldTransfer

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.H

@@ -108,6 +108,51 @@ public:
             Istream& is
         );
 
+        //- Construct from fine level interface,
+        //- local and neighbour restrict addressing
+        cyclicGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const lduInterface& fineInterface,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces
+        );
+
+        //- Construct by assembling and returning a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            return autoPtr<GAMGInterface>
+            (
+                new cyclicGAMGInterface
+                (
+                    index,
+                    coarseInterfaces,
+                    *this,
+                    interfaceMap,
+                    faceCells,
+                    faceRestrictAddresssing,
+                    faceOffsets,
+                    allInterfaces
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicGAMGInterface() = default;

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -28,6 +28,7 @@ License
 
 #include "lduPrimitiveMesh.H"
 #include "processorLduInterface.H"
+#include "processorCyclicGAMGInterface.H"
 #include "edgeHashes.H"
 #include "labelPair.H"
 #include "processorGAMGInterface.H"
@@ -308,6 +309,56 @@ Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
 }
 
 
+Foam::label Foam::lduPrimitiveMesh::findConnectedInterface
+(
+    const lduMesh& myMesh,
+    const PtrList<lduPrimitiveMesh>& otherMeshes,
+    const labelPairList& procAndInterfaces,
+
+    const label nbrProci,
+    const label myRank
+) const
+{
+    // Find mesh, interfacei in procAndInterfaces
+    label nbrInti = -1;
+    for (const auto& procAndInterface : procAndInterfaces)
+    {
+        const label proci = procAndInterface[0];
+
+        if (proci == nbrProci)
+        {
+            const label interfacei = procAndInterface[1];
+            const lduInterfacePtrsList interfaces =
+                mesh
+                (
+                    myMesh,
+                    otherMeshes,
+                    proci
+                ).interfaces();
+            const processorLduInterface& pldui =
+                refCast<const processorLduInterface>
+                (
+                    interfaces[interfacei]
+                );
+
+            if (pldui.neighbProcNo() == myRank)
+            {
+                nbrInti = procAndInterface[1];
+                break;
+            }
+        }
+    }
+
+
+    if (nbrInti == -1)
+    {
+        FatalErrorInFunction
+            << "procAndInterfaces:" << procAndInterfaces << abort(FatalError);
+    }
+    return nbrInti;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
@@ -472,7 +523,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     const label nMeshes = otherMeshes.size()+1;
 
     const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
-
     if (lduPrimitiveMesh::debug)
     {
         Pout<< "I am " << UPstream::myProcNo(currentComm)
@@ -481,6 +531,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
             << endl;
     }
 
+    const lduInterfacePtrsList myInterfaces = myMesh.interfaces();
 
     forAll(procIDs, i)
     {
@@ -522,24 +573,28 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
           + procMesh.lduAddr().lowerAddr().size();
     }
 
-    // Count how faces get added. Interfaces inbetween get merged.
+    // Count how faces get added. Proc interfaces inbetween get merged.
 
-    // Merged interfaces: map from two coarse processors back to
+    // Merged proc interfaces: map from two coarse processors back to
     // - procMeshes
     // - interface in procMesh
     // (estimate size from number of patches of mesh0)
-    EdgeMap<labelPairList> mergedMap(2*myMesh.interfaces().size());
+    EdgeMap<labelPairList> mergedMap(2*myInterfaces.size());
 
-    // Unmerged interfaces: map from two coarse processors back to
+    // Unmerged proc interfaces: map from two coarse processors back to
     // - procMeshes
     // - interface in procMesh
     EdgeMap<labelPairList> unmergedMap(mergedMap.size());
 
+    // (unmerged) global interfaces. These are present on all processors
+    // in the same order (and keep the order in the merged mesh)
+    List<DynamicList<label>> procToGlobal(nMeshes);
+
+
     boundaryMap.setSize(nMeshes);
     boundaryFaceMap.setSize(nMeshes);
 
 
-    label nOtherInterfaces = 0;
     labelList nCoupledFaces(nMeshes, Zero);
 
     for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
@@ -560,6 +615,16 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
                 if (isA<processorLduInterface>(ldui))
                 {
+                    if (isA<processorCyclicGAMGInterface>(ldui))
+                    {
+                        FatalErrorInFunction
+                            << "At mesh from processor " << procIDs[procMeshI]
+                            << " have interface " << intI
+                            << " of unhandled type " << ldui.type()
+                            << ". Adapt decomposition to avoid these"
+                            << exit(FatalError);
+                    }
+
                     const processorLduInterface& pldui =
                         refCast<const processorLduInterface>(ldui);
 
@@ -582,7 +647,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                         // Merged interface
                         if (debug)
                         {
-                            Pout<< "merged interface: myProcNo:"
+                            Pout<< "merged proc interface: myProcNo:"
                                 << pldui.myProcNo()
                                 << " nbr:" << pldui.neighbProcNo()
                                 << " size:" << ldui.faceCells().size()
@@ -607,7 +672,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                     {
                         if (debug)
                         {
-                            Pout<< "external interface: myProcNo:"
+                            Pout<< "external proc interface: myProcNo:"
                                 << pldui.myProcNo()
                                 << " nbr:" << pldui.neighbProcNo()
                                 << " size:" << ldui.faceCells().size()
@@ -623,53 +688,90 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 else
                 {
                     // Still external (non proc) interface
-                    FatalErrorInFunction
-                        << "At mesh from processor " << procIDs[procMeshI]
-                        << " have interface " << intI
-                        << " of unhandled type " << interfaces[intI].type()
-                        << exit(FatalError);
+                    //FatalErrorInFunction
+                    //WarningInFunction
+                    //    << "At mesh from processor " << procIDs[procMeshI]
+                    //    << " have interface " << intI
+                    //    << " of unhandled type " << ldui.type()
+                    //    << " with size:" << ldui.faceCells().size()
+                    //    //<< exit(FatalError);
+                    //    << endl;
 
-                    ++nOtherInterfaces;
+                    procToGlobal[procMeshI].append(intI);
                 }
             }
         }
     }
 
 
-
+    // Check that all processors have any global patches in the same order
     if (debug)
     {
-        Pout<< "Remaining interfaces:" << endl;
-        forAllConstIters(unmergedMap, iter)
+        const auto& global0 = procToGlobal[0];
+        for (label procMeshI = 1; procMeshI < nMeshes; ++procMeshI)
         {
-            Pout<< "    agglom procEdge:" << iter.key() << endl;
-            const labelPairList& elems = iter.val();
-            forAll(elems, i)
+            const auto& global = procToGlobal[procMeshI];
+            if (global != global0)
             {
-                label procMeshI = elems[i][0];
-                label interfacei = elems[i][1];
-                const lduInterfacePtrsList interfaces =
-                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
-
-                const processorLduInterface& pldui =
-                    refCast<const processorLduInterface>
-                    (
-                        interfaces[interfacei]
-                    );
-
-                Pout<< "        proc:" << procIDs[procMeshI]
-                    << " interfacei:" << interfacei
-                    << " between:" << pldui.myProcNo()
-                    << " and:" << pldui.neighbProcNo()
-                    << endl;
+                FatalErrorInFunction
+                    << "At mesh from processor " << procIDs[procMeshI]
+                    << " have global interfaces " << global
+                    << " which differ from those on processor "
+                    << procIDs[procMeshI]
+                    << " : " << global0
+                    << exit(FatalError);
             }
-            Pout<< endl;
         }
     }
+
+
     if (debug)
     {
+        Pout<< "Global interfaces:" << endl;
+        const auto& global0 = procToGlobal[0];
+        for (const label intI : global0)
+        {
+            Pout<< "        interfacei:" << intI
+                << " type:" << myInterfaces[intI].type()
+                << endl;
+
+        }
+        Pout<< endl;
+
+        Pout<< "Remaining interfaces:" << endl;
+        //forAllConstIters(unmergedMap, iter)
+        for (const auto& iter : unmergedMap.csorted())
+        {
+            Pout<< "    agglom procEdge:" << iter.key() << endl;
+            const labelPairList& elems = iter.val();
+            forAll(elems, i)
+            {
+                label procMeshI = elems[i][0];
+                label interfacei = elems[i][1];
+                const lduInterfacePtrsList interfaces =
+                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
+
+                const processorLduInterface& pldui =
+                    refCast<const processorLduInterface>
+                    (
+                        interfaces[interfacei]
+                    );
+
+                Pout<< "        proc:" << procIDs[procMeshI]
+                    << " interfacei:" << interfacei
+                    << " between:" << pldui.myProcNo()
+                    << " and:" << pldui.neighbProcNo()
+                    << " localsize:"
+                    << interfaces[interfacei].faceCells().size()
+                    << endl;
+            }
+            Pout<< endl;
+        }
+        Pout<< endl;
+
         Pout<< "Merged interfaces:" << endl;
-        forAllConstIters(mergedMap, iter)
+        //forAllConstIters(mergedMap, iter)
+        for (const auto& iter : mergedMap.csorted())
         {
             Pout<< "    agglom procEdge:" << iter.key() << endl;
             const labelPairList& elems = iter.val();
@@ -690,10 +792,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                     << " interfacei:" << interfacei
                     << " between:" << pldui.myProcNo()
                     << " and:" << pldui.neighbProcNo()
+                    << " localsize:"
+                    << interfaces[interfacei].faceCells().size()
                     << endl;
             }
             Pout<< endl;
         }
+        Pout<< endl;
     }
 
 
@@ -780,43 +885,14 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                         {
                             const labelPairList& elems = fnd();
 
-                            // Find nbrP in elems
-                            label nbrIntI = -1;
-                            forAll(elems, i)
-                            {
-                                label proci = elems[i][0];
-                                label interfacei = elems[i][1];
-                                const lduInterfacePtrsList interfaces =
-                                    mesh
-                                    (
-                                        myMesh,
-                                        otherMeshes,
-                                        proci
-                                    ).interfaces();
-                                const processorLduInterface& pldui =
-                                    refCast<const processorLduInterface>
-                                    (
-                                        interfaces[interfacei]
-                                    );
-
-                                if
-                                (
-                                    elems[i][0] == nbrProcMeshI
-                                 && pldui.neighbProcNo() == procIDs[procMeshI]
-                                )
-                                {
-                                    nbrIntI = elems[i][1];
-                                    break;
-                                }
-                            }
-
-
-                            if (nbrIntI == -1)
-                            {
-                                FatalErrorInFunction
-                                    << "elems:" << elems << abort(FatalError);
-                            }
-
+                            const label nbrIntI = findConnectedInterface
+                            (
+                                myMesh,
+                                otherMeshes,
+                                elems,
+                                nbrProcMeshI,
+                                procIDs[procMeshI]
+                            );
 
                             const lduInterfacePtrsList nbrInterfaces = mesh
                             (
@@ -930,12 +1006,118 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // Kept interfaces
     // ~~~~~~~~~~~~~~~
 
-    interfaces_.setSize(unmergedMap.size() + nOtherInterfaces);
+    interfaces_.setSize(unmergedMap.size() + procToGlobal[0].size());
     primitiveInterfaces_.setSize(interfaces_.size());
 
     label allInterfacei = 0;
 
-    forAllConstIters(unmergedMap, iter)
+
+    // Global interfaces
+    // ~~~~~~~~~~~~~~~~~
+    // (e.g. cyclicAMI)
+
+    {
+        const auto& global0 = procToGlobal[0];
+
+        for (const label interfacei : global0)
+        {
+            //Pout<< "        interfacei:" << interfacei
+            //    << " type:" << interfaces[interfacei].type()
+            //    << endl;
+
+            // Just add all individual face-cells in processor order
+
+            // Count
+            label n = 0;
+            for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
+            {
+                const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+                n += procMesh.interfaces()[interfacei].faceCells().size();
+            }
+
+            // Size
+            labelField allFaceCells(n);
+            labelField allFaceRestrictAddressing(n);
+            labelList faceOffsets(nMeshes+1, 0);
+            lduInterfacePtrsList allProcInterfaces(nMeshes);
+            n = 0;
+
+            // Fill
+            for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
+            {
+                faceOffsets[procMeshI] = n;
+
+                const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+                const lduInterfacePtrsList interfaces = procMesh.interfaces();
+
+                allProcInterfaces.set(procMeshI, &interfaces[interfacei]);
+                boundaryMap[procMeshI][interfacei] = allInterfacei;
+                labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];
+
+                const labelUList& l = interfaces[interfacei].faceCells();
+                bfMap.setSize(l.size());
+
+                forAll(l, facei)
+                {
+                    allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
+                    allFaceRestrictAddressing[n] = n;
+                    bfMap[facei] = n;
+                    n++;
+                }
+            }
+
+            // For convenience populate last element with size
+            faceOffsets.last() = n;
+
+
+            const auto& myFineInterface =
+                refCast<const GAMGInterface>(myInterfaces[interfacei]);
+
+            autoPtr<GAMGInterface> ppPtr
+            (
+                myFineInterface.clone
+                (
+                    allInterfacei,
+                    interfaces_,
+                    global0,
+                    allFaceCells,
+                    allFaceRestrictAddressing,
+                    faceOffsets,
+                    allProcInterfaces,
+                    comm_,
+                    myAgglom,
+                    procAgglomMap
+                )
+            );
+
+            primitiveInterfaces_.set
+            (
+                allInterfacei,
+                ppPtr.ptr()
+            );
+
+            interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);
+
+            if (debug)
+            {
+                Pout<< "Created " << interfaces_[allInterfacei].type()
+                    << " interface at " << allInterfacei
+                    << " comm:" << comm_
+                    << " myProcNo:" << myAgglom
+                    << " nFaces:" << allFaceCells.size()
+                    << endl;
+            }
+
+            allInterfacei++;
+        }
+    }
+
+
+    // Unmerged processor interfaces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //forAllConstIters(unmergedMap, iter)
+    for (const auto& iter : unmergedMap.csorted())
     {
         const labelPairList& elems = iter.val();
 
@@ -1085,6 +1267,58 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     }
 
 
+    if (allInterfacei != interfaces_.size())
+    {
+        FatalErrorInFunction << "allInterfacei:" << allInterfacei
+            << " interfaces_:" << interfaces_.size() << exit(FatalError);
+    }
+
+
+    if (debug)
+    {
+        Pout<< endl
+            << "Created new lduPrimitiveMesh:" << nl
+            << "    cells:" << this->lduAddr().size() << nl
+            << "    internal face lower:"
+            << this->lduAddr().lowerAddr().size() << nl
+            << "    internal faces upper:"
+            << this->lduAddr().upperAddr().size()
+            << endl;
+        forAll(interfaces_, i)
+        {
+            if (interfaces_.set(i))
+            {
+                Pout<< "    interface:" << i << " type:" << interfaces_[i].type()
+                    << nl
+                    << "    faceCells:" << flatOutput(interfaces_[i].faceCells())
+                    << endl;
+            }
+        }
+
+
+        Pout<< "Original input meshes:" << endl;
+        forAll(boundaryMap, procMeshI)
+        {
+            const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+            const lduInterfacePtrsList interfaces = procMesh.interfaces();
+
+            Pout<< "    proc:" << procMeshI
+                << " interfaces:" << interfaces.size() << endl;
+            forAll(interfaces, inti)
+            {
+                if (interfaces.set(inti))
+                {
+                    Pout<< "        int:" << inti
+                        << " type:" << interfaces[inti].type()
+                        << " size:" << interfaces[inti].faceCells().size()
+                        << " maps to:" << boundaryMap[procMeshI][inti]
+                        << endl;
+                }
+            }
+        }
+    }
+
+
     patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
 
     if (debug)
@@ -1111,89 +1345,25 @@ void Foam::lduPrimitiveMesh::gather
 (
     const label comm,
     const lduMesh& mesh,
-    const labelList& procIDs,
     PtrList<lduPrimitiveMesh>& otherMeshes
 )
 {
     // Force calculation of schedule (since does parallel comms)
     (void)mesh.lduAddr().patchSchedule();
 
-    if (Pstream::myProcNo(comm) == procIDs[0])
+    // Use PstreamBuffers
+    PstreamBuffers pBufs
+    (
+        Pstream::commsTypes::nonBlocking,
+        UPstream::msgType(),
+        comm
+    );
+
+    // Send to master
+    if (!Pstream::master(comm))
     {
-        // Master.
-        otherMeshes.setSize(procIDs.size()-1);
-
-        for (label i = 1; i < procIDs.size(); ++i)
-        {
-            //Pout<< "on master :"
-            //    << " receiving from proc " << procIDs[i] << endl;
-
-            IPstream fromProc
-            (
-                Pstream::commsTypes::scheduled,
-                procIDs[i],
-                0,          // bufSize
-                Pstream::msgType(),
-                comm
-            );
-
-            label nCells = readLabel(fromProc);
-            labelList lowerAddr(fromProc);
-            labelList upperAddr(fromProc);
-            boolList validInterface(fromProc);
-
-
-            // Construct mesh without interfaces
-            otherMeshes.set
-            (
-                i-1,
-                new lduPrimitiveMesh
-                (
-                    nCells,
-                    lowerAddr,
-                    upperAddr,
-                    comm,
-                    true    // reuse
-                )
-            );
-
-            // Construct GAMGInterfaces
-            lduInterfacePtrsList newInterfaces(validInterface.size());
-            forAll(validInterface, intI)
-            {
-                if (validInterface[intI])
-                {
-                    word coupleType(fromProc);
-
-                    newInterfaces.set
-                    (
-                        intI,
-                        GAMGInterface::New
-                        (
-                            coupleType,
-                            intI,
-                            otherMeshes[i-1].rawInterfaces(),
-                            fromProc
-                        ).ptr()
-                    );
-                }
-            }
-
-            otherMeshes[i-1].addInterfaces
-            (
-                newInterfaces,
-                nonBlockingSchedule<processorGAMGInterface>
-                (
-                    newInterfaces
-                )
-            );
-       }
-    }
-    else if (procIDs.found(Pstream::myProcNo(comm)))
-    {
-        // Send to master
-
         const lduAddressing& addressing = mesh.lduAddr();
+
         lduInterfacePtrsList interfaces(mesh.interfaces());
         boolList validInterface(interfaces.size());
         forAll(interfaces, intI)
@@ -1201,14 +1371,7 @@ void Foam::lduPrimitiveMesh::gather
             validInterface[intI] = interfaces.set(intI);
         }
 
-        OPstream toMaster
-        (
-            Pstream::commsTypes::scheduled,
-            procIDs[0],
-            0,
-            Pstream::msgType(),
-            comm
-        );
+        UOPstream toMaster(UPstream::masterNo(), pBufs);
 
         toMaster
             << addressing.size()
@@ -1230,6 +1393,74 @@ void Foam::lduPrimitiveMesh::gather
             }
         }
     }
+
+    // Wait for finish
+    pBufs.finishedGathers();
+
+    // Consume
+    if (Pstream::master(comm))
+    {
+        const label nProcs = UPstream::nProcs(comm);
+
+        // Master.
+        otherMeshes.setSize(nProcs-1);
+
+        for (const int proci : UPstream::subProcs(comm))
+        {
+            UIPstream fromProc(proci, pBufs);
+
+            const label nCells = readLabel(fromProc);
+            labelList lowerAddr(fromProc);
+            labelList upperAddr(fromProc);
+            const boolList validInterface(fromProc);
+
+
+            // Construct mesh without interfaces
+            otherMeshes.set
+            (
+                proci-1,
+                new lduPrimitiveMesh
+                (
+                    nCells,
+                    lowerAddr,
+                    upperAddr,
+                    mesh.comm(),
+                    true    // reuse
+                )
+            );
+
+            // Construct GAMGInterfaces
+            lduInterfacePtrsList newInterfaces(validInterface.size());
+            forAll(validInterface, intI)
+            {
+                if (validInterface[intI])
+                {
+                    word coupleType(fromProc);
+
+                    newInterfaces.set
+                    (
+                        intI,
+                        GAMGInterface::New
+                        (
+                            coupleType,
+                            intI,
+                            otherMeshes[proci-1].rawInterfaces(),
+                            fromProc
+                        ).ptr()
+                    );
+                }
+            }
+
+            otherMeshes[proci-1].addInterfaces
+            (
+                newInterfaces,
+                nonBlockingSchedule<processorGAMGInterface>
+                (
+                    newInterfaces
+                )
+            );
+        }
+    }
 }
 
 

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -82,6 +82,17 @@ class lduPrimitiveMesh
         //- Get size of all meshes
         static label totalSize(const PtrList<lduPrimitiveMesh>&);
 
+        label findConnectedInterface
+        (
+            const lduMesh& myMesh,
+            const PtrList<lduPrimitiveMesh>& otherMeshes,
+            const labelPairList& procAndInterfaces,
+
+            const label nbrProci,
+            const label myRank
+        ) const;
+
+
         //- No copy construct
         lduPrimitiveMesh(const lduPrimitiveMesh&) = delete;
 
@@ -277,13 +288,12 @@ public:
                 const label meshI
             );
 
-            //- Gather meshes from other processors onto procIDs[0].
-            //  Received meshes get GAMGInterface and communicator comm
+            //- Gather meshes from other processors using agglomComm.
+            //  Received meshes get GAMGInterface.
             static void gather
             (
-                const label comm,
+                const label agglomComm,
                 const lduMesh& mesh,
-                const labelList& procIDs,
                 PtrList<lduPrimitiveMesh>& otherMeshes
             );
 

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMeshTemplates.C

@@ -84,6 +84,8 @@ Foam::lduSchedule Foam::lduPrimitiveMesh::nonBlockingSchedule
     // (where interfaces are defined) but must retain the full list length
     // for later (external) bookkeeping
 
+    schedule.setSize(patchEvali);
+
     return schedule;
 }
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -87,7 +87,12 @@ void Foam::globalMeshData::initProcAddr()
 
     if (Pstream::parRun())
     {
-        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+        PstreamBuffers pBufs
+        (
+            Pstream::commsTypes::nonBlocking,
+            UPstream::msgType(),
+            mesh_.comm()
+        );
 
         // Send indices of my processor patches to my neighbours
         for (const label patchi : processorPatches_)
@@ -1738,7 +1743,7 @@ Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
         processorTopology::New<processorPolyPatch>
         (
             mesh.boundaryMesh(),
-            UPstream::worldComm
+            mesh_.comm()
         )
     ),
     processorPatches_(),

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -146,6 +146,21 @@ void Foam::mapDistribute::printLayout(Ostream& os) const
 }
 
 
+Foam::UPtrList<const Foam::mapDistributeBase> Foam::mapDistribute::extractBase
+(
+    const UPtrList<const mapDistribute>& maps
+)
+{
+    UPtrList<const mapDistributeBase> baseMaps(maps.size());
+    forAll(maps, i)
+    {
+        const mapDistributeBase& map = maps[i];
+        baseMaps.set(i, &map);
+    }
+    return baseMaps;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::mapDistribute::mapDistribute() noexcept
@@ -437,6 +452,32 @@ Foam::mapDistribute::mapDistribute
 }
 
 
+Foam::mapDistribute::mapDistribute
+(
+    const UPtrList<const mapDistribute>& maps,
+    const labelList& localRanks,
+    const label newComm,
+    const labelListList& newToOldRanks, // from rank in newComm to
+                                        // ranks in (old)comm
+    labelList& startOfLocal,
+    List<Map<label>>& compactMaps
+)
+:
+    mapDistributeBase
+    (
+        extractBase(maps),
+        localRanks,
+        newComm,
+        newToOldRanks,
+        startOfLocal,
+        compactMaps
+    )
+{
+    // TBD. -have mapDistributeBase::add or something
+    //      -set transforms from individual maps
+}
+
+
 Foam::autoPtr<Foam::mapDistribute> Foam::mapDistribute::clone() const
 {
     return autoPtr<mapDistribute>::New(*this);

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H

@@ -198,6 +198,12 @@ class mapDistribute
             const TransformOp& top
         ) const;
 
+        //- Helper: convert mapDistribute to mapDistributeBase
+        static UPtrList<const mapDistributeBase> extractBase
+        (
+            const UPtrList<const mapDistribute>& maps
+        );
+
 
 public:
 
@@ -393,6 +399,27 @@ public:
             const label comm = UPstream::worldComm
         );
 
+        //- Construct from multiple maps and processor collation
+        //  Assumes all local data first. Sorts contributions of maps
+        //  in processor order i.e.  constructed map has all local data first.
+        //  Returns
+        //  - startOfLocal : per input map the start of the local data. Extends
+        //                   one beyond number of maps so overall local size
+        //                   is startOfLocal.last()
+        //  - compactMaps  : per input map from slot position in the input map
+        //                   to new slot position. (note there is no information
+        //                   returned about which processor it is from)
+        mapDistribute
+        (
+            const UPtrList<const mapDistribute>& maps,
+            const labelList& localRanks,
+            const label newComm,
+            const labelListList& newToOldRanks, // from rank in newComm to
+                                                // ranks in (old)comm
+            labelList& startOfLocal,            // per map start of local data
+            List<Map<label>>& compactMaps       // per map old slot to new slot
+        );
+
         //- Construct from Istream
         explicit mapDistribute(Istream& is);
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBase.C

@@ -399,9 +399,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
 
     for (const label globalIdx : elements)
     {
-        if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+        if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
         {
-            label proci = globalNumbering.whichProcID(globalIdx);
+            label proci = globalNumbering.whichProcID(myRank, globalIdx);
             nNonLocal[proci]++;
         }
     }
@@ -421,9 +421,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
     // Collect all (non-local) elements needed.
     for (const label globalIdx : elements)
     {
-        if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+        if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
         {
-            label proci = globalNumbering.whichProcID(globalIdx);
+            label proci = globalNumbering.whichProcID(myRank, globalIdx);
             label index = globalNumbering.toLocal(proci, globalIdx);
             label nCompact = compactMap[proci].size();
             compactMap[proci].insert(index, nCompact);
@@ -449,9 +449,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
     {
         for (const label globalIdx : cCells)
         {
-            if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+            if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
             {
-                label proci = globalNumbering.whichProcID(globalIdx);
+                label proci = globalNumbering.whichProcID(myRank, globalIdx);
                 nNonLocal[proci]++;
             }
         }
@@ -474,9 +474,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
     {
         for (const label globalIdx : cCells)
         {
-            if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+            if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
             {
-                label proci = globalNumbering.whichProcID(globalIdx);
+                label proci = globalNumbering.whichProcID(myRank, globalIdx);
                 label index = globalNumbering.toLocal(proci, globalIdx);
                 label nCompact = compactMap[proci].size();
                 compactMap[proci].insert(index, nCompact);
@@ -504,7 +504,7 @@ void Foam::mapDistributeBase::exchangeAddressing
 
     compactStart.setSize(nProcs);
     compactStart[myRank] = 0;
-    constructSize_ = globalNumbering.localSize();
+    constructSize_ = globalNumbering.localSize(myRank);
     forAll(compactStart, proci)
     {
         if (proci != myRank)
@@ -526,7 +526,7 @@ void Foam::mapDistributeBase::exchangeAddressing
         if (proci == myRank)
         {
             // All my own elements are used
-            label nLocal = globalNumbering.localSize();
+            label nLocal = globalNumbering.localSize(myRank);
             wantedRemoteElements[proci] = identity(nLocal);
             constructMap_[proci] = identity(nLocal);
         }
@@ -561,7 +561,7 @@ void Foam::mapDistributeBase::exchangeAddressing
     // Renumber elements
     for (label& elem : elements)
     {
-        elem = renumber(globalNumbering, compactMap, elem);
+        elem = renumber(globalNumbering, comm_, compactMap, elem);
     }
 }
 
@@ -584,7 +584,7 @@ void Foam::mapDistributeBase::exchangeAddressing
 
     compactStart.setSize(nProcs);
     compactStart[myRank] = 0;
-    constructSize_ = globalNumbering.localSize();
+    constructSize_ = globalNumbering.localSize(myRank);
     forAll(compactStart, proci)
     {
         if (proci != myRank)
@@ -606,7 +606,7 @@ void Foam::mapDistributeBase::exchangeAddressing
         if (proci == myRank)
         {
             // All my own elements are used
-            label nLocal = globalNumbering.localSize();
+            label nLocal = globalNumbering.localSize(myRank);
             wantedRemoteElements[proci] = identity(nLocal);
             constructMap_[proci] = identity(nLocal);
         }
@@ -643,7 +643,7 @@ void Foam::mapDistributeBase::exchangeAddressing
     {
         for (label& celli : cCells)
         {
-            celli = renumber(globalNumbering, compactMap, celli);
+            celli = renumber(globalNumbering, comm_, compactMap, celli);
         }
     }
 }
@@ -1005,6 +1005,820 @@ Foam::mapDistributeBase::mapDistributeBase
 {}
 
 
+Foam::mapDistributeBase::mapDistributeBase
+(
+    const UPtrList<const mapDistributeBase>& maps,
+    const labelList& localRanks,
+    const label newComm,
+    const labelListList& newToOldRanks,// from newComm to comm_
+    labelList& startOfLocal,
+    List<Map<label>>& compactMaps
+)
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    comm_(-1),
+    schedulePtr_(nullptr)
+{
+    if (maps.empty())
+    {
+        return;
+    }
+
+    comm_ = newComm;
+    subHasFlip_ = maps[0].subHasFlip();
+    constructHasFlip_ = maps[0].constructHasFlip();
+
+    const label nNewRanks = newToOldRanks.size();
+    const label myNewRank = UPstream::myProcNo(newComm);
+    if (nNewRanks != UPstream::nProcs(newComm))
+    {
+        FatalErrorInFunction<< "nNewRanks:" << nNewRanks
+            << " nProcs:" << UPstream::nProcs(newComm)
+            << exit(FatalError);
+    }
+
+    if (localRanks.size() != maps.size())
+    {
+        FatalErrorInFunction
+            << "Number of maps:" << maps.size()
+            << " number of localRanks:" << localRanks.size()
+            << exit(FatalError);
+    }
+
+    // Sanity checks
+    const auto& map0 = maps[0];
+    forAll(maps, mapi)
+    {
+        const auto& map = maps[mapi];
+
+        if
+        (
+            (map.comm() != map0.comm())
+         || (map.subHasFlip() != map0.subHasFlip())
+         || (map.constructHasFlip() != map0.constructHasFlip())
+        )
+        {
+            FatalErrorInFunction
+                << "Maps should all be the same form"
+                << " Map " << mapi
+                << " has comm:" << map.comm()
+                << " subHasFlip:" << map.subHasFlip()
+                << " constructHasFlip:" << map.constructHasFlip()
+                << " which is different from map 0"
+                << exit(FatalError);
+        }
+
+        const label localRank = localRanks[mapi];
+        const auto& constructOwn = maps[mapi].constructMap()[localRank];
+        forAll(constructOwn, i)
+        {
+            if (constructOwn[i] != i)
+            {
+                FatalErrorInFunction
+                    << "Maps constructMap not identity."
+                    << " Map " << mapi
+                    << " constructMap:" << flatOutput(constructOwn)
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    constructMap_.resize_nocopy(nNewRanks);
+    subMap_.resize_nocopy(nNewRanks);
+
+
+    // Store starts
+    startOfLocal.setSize(maps.size()+1);
+    compactMaps.resize_nocopy(maps.size());
+
+    label constructi = 0;
+    forAll(maps, mapi)
+    {
+        startOfLocal[mapi] = constructi;
+        const label localRank = localRanks[mapi];
+        const auto& map = maps[mapi].constructMap()[localRank];
+
+        // Presize compaction array
+        const label nRemote = maps[mapi].constructSize()-map.size();
+        compactMaps[mapi].resize(2*nRemote);
+
+        constructi += map.size();
+    }
+    startOfLocal.last() = constructi;
+
+
+    // Determine start of constructed remote data. This is used to get the
+    // local offset which can then be used to get the relative subMap location.
+    labelListList startOfRemote(maps.size());
+    forAll(maps, mapi)
+    {
+        const label nOldProcs = maps[mapi].constructMap().size();
+        labelList& starts = startOfRemote[mapi];
+
+        starts.setSize(nOldProcs, labelMax);
+        forAll(maps[mapi].constructMap(), oldProci)
+        {
+            const labelList& map = maps[mapi].constructMap()[oldProci];
+            forAll(map, i)
+            {
+                const label index
+                (
+                    constructHasFlip_
+                  ? mag(map[i])-1
+                  : map[i]
+                );
+                starts[oldProci] = min(starts[oldProci], index);
+            }
+        }
+    }
+
+
+    // Construct map
+    // ~~~~~~~~~~~~~
+    // - all localRanks:
+    //      - data gets appended in map order
+    //      - map is just an offset (startOfLocal)
+    // - all previously remote ranks:
+    //      - data is already present according to startOfLocal
+    //      - map is old-to-new index
+    // - all still remote ranks:
+    //      - data gets appended in map order after the startOfLocal
+    //      - map is old-to-new index
+
+
+    // Append local (= myRank) data. TBD: assumes subMap and constructMap
+    // are identity maps.
+    {
+        labelList& myConstruct = constructMap_[myNewRank];
+        myConstruct.resize_nocopy(constructi);
+        constructi = 0;
+        forAll(maps, mapi)
+        {
+            const label localRank = localRanks[mapi];
+            const auto& map = maps[mapi].constructMap()[localRank];
+            const label offset = startOfLocal[mapi];
+
+            forAll(map, i)
+            {
+                if (constructHasFlip_)
+                {
+                    forAll(map, i)
+                    {
+                        if (map[i] < 0)
+                        {
+                            myConstruct[constructi++] = map[i]-offset;
+                        }
+                        else
+                        {
+                            myConstruct[constructi++] = map[i]+offset;
+                        }
+                    }
+                }
+                else
+                {
+                    myConstruct[constructi++] = map[i]+offset;
+                }
+            }
+        }
+    }
+
+    // Filter remote construct data
+    {
+        // Remote ranks that are now local
+        //  - store new index for mapping stencils
+        //  - no need to construct since already
+        const auto& oldProcs = newToOldRanks[myNewRank];
+
+        forAll(maps, mapi)
+        {
+            for (const label oldProci : oldProcs)
+            {
+                if (oldProci != localRanks[mapi])
+                {
+                    const auto& map = maps[mapi].constructMap()[oldProci];
+
+                    if (!map.size())
+                    {
+                        continue;
+                    }
+
+
+                    // The slots come from a local map so we can look up the
+                    // new location
+                    const label sourceMapi = localRanks.find(oldProci);
+                    const auto& subMap =
+                        maps[sourceMapi].subMap()[localRanks[mapi]];
+
+                    //Pout<< "From oldRank:" << oldProci
+                    //    << " sending to masterRank:" << localRanks[mapi]
+                    //    << " elements:" << flatOutput(subMap)
+                    //    << nl
+                    //    << "   received as elements:" << flatOutput(map)
+                    //    << endl;
+
+                    if (map.size() != subMap.size())
+                    {
+                        FatalErrorInFunction << "Problem:"
+                            << "oldProci:" << oldProci
+                            << " mapi:" << mapi
+                            << " constructMap:" << map.size()
+                            << " sourceMapi:" << sourceMapi
+                            << " subMap:" << subMap.size()
+                            << exit(FatalError);
+                    }
+
+                    const label offset = startOfLocal[sourceMapi];
+                    // Construct map starts after the local data
+                    const label nMapLocal = startOfRemote[mapi][oldProci];
+
+                    auto& cptMap = compactMaps[mapi];
+                    forAll(map, i)
+                    {
+                        // old slot position to new slot position
+                        const label index
+                        (
+                            constructHasFlip_
+                          ? mag(map[i])-1
+                          : map[i]
+                        );
+                        const label newIndex = subMap[index-nMapLocal]+offset;
+
+                        // Note: should always warn for duplicates? Or only if
+                        // different?
+                        if
+                        (
+                           !cptMap.insert(index, newIndex)
+                         && cptMap[index] != newIndex
+                        )
+                        {
+                            FatalErrorInFunction<< "Duplicate insertion"
+                                << "From oldProc:" << oldProci
+                                << " on map:" << mapi
+                                << " at index:" << i
+                                << " have construct slot:" << index
+                                << " new index:" << newIndex
+                                << " but already have entry:" << cptMap[index]
+                                << " on for that slot"
+                                << exit(FatalError);
+                        }
+                    }
+                }
+            }
+        }
+
+
+        // Remote ranks that are still remote
+        //  - store new index for mapping stencils
+        //  - append to construction
+
+        // Either loop over all old ranks and filter out ones already handled
+        // or loop over all new ranks and avoid myNewRank
+
+        forAll(newToOldRanks, newProci)
+        {
+            if (newProci != myNewRank)
+            {
+                const auto& oldProcs = newToOldRanks[newProci];
+
+                label allSize = 0;
+                forAll(maps, mapi)
+                {
+                    for (const label oldProci : oldProcs)
+                    {
+                        allSize += maps[mapi].constructMap()[oldProci].size();
+                    }
+                }
+
+                labelList& myConstruct = constructMap_[newProci];
+                myConstruct.resize_nocopy(allSize);
+
+                allSize = 0;
+                forAll(maps, mapi)
+                {
+                    for (const label oldProci : oldProcs)
+                    {
+                        const auto& map = maps[mapi].constructMap()[oldProci];
+                        // Construct map starts after the local data
+                        const label nMapLocal = startOfRemote[mapi][oldProci];
+                        SubList<label> slice(myConstruct, map.size(), allSize);
+
+                        if (constructHasFlip_)
+                        {
+                            forAll(map, i)
+                            {
+                                if (map[i] < 0)
+                                {
+                                    slice[i] = map[i]+nMapLocal-constructi;
+                                }
+                                else
+                                {
+                                    slice[i] = map[i]-nMapLocal+constructi;
+                                }
+                            }
+
+                            auto& cptMap = compactMaps[mapi];
+                            forAll(map, i)
+                            {
+                                cptMap.insert(mag(map[i])-1,mag(slice[i])-1);
+                            }
+                        }
+                        else
+                        {
+                            forAll(map, i)
+                            {
+                                slice[i] = map[i]-nMapLocal+constructi;
+                                compactMaps[mapi].insert(map[i], slice[i]);
+                            }
+                        }
+                        allSize += map.size();
+                        constructi += map.size();
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Sub (=send) map
+    // ~~~~~~~~~~~~~~~
+    // - all localRanks:
+    //      - get appended in map order
+    // - all previously remote ranks:
+    //      - not needed. Stay empty
+    // - all still remote ranks:
+    //      - convert to new local index
+
+    // Append local (= myRank) data
+    {
+        label allSize = 0;
+        forAll(maps, mapi)
+        {
+            const label localRank = localRanks[mapi];
+            allSize += maps[mapi].subMap()[localRank].size();
+        }
+
+        labelList& mySub = subMap_[myNewRank];
+        mySub.resize_nocopy(allSize);
+        allSize = 0;
+        forAll(maps, mapi)
+        {
+            const label localRank = localRanks[mapi];
+            const auto& map = maps[mapi].subMap()[localRank];
+            SubList<label> slice(mySub, map.size(), allSize);
+
+            if (subHasFlip_)
+            {
+                forAll(slice, i)
+                {
+                    if (map[i] < 0)
+                    {
+                        slice[i] = map[i]-startOfLocal[mapi];
+                    }
+                    else
+                    {
+                        slice[i] = map[i]+startOfLocal[mapi];
+                    }
+                }
+            }
+            else
+            {
+                forAll(slice, i)
+                {
+                    slice[i] = map[i]+startOfLocal[mapi];
+                }
+            }
+            allSize += map.size();
+        }
+    }
+    // Filter remote sub data
+    forAll(newToOldRanks, newProci)
+    {
+        if (newProci != myNewRank)
+        {
+            const auto& oldProcs = newToOldRanks[newProci];
+
+            label allSize = 0;
+            forAll(maps, mapi)
+            {
+                for (const label oldProci : oldProcs)
+                {
+                    allSize += maps[mapi].subMap()[oldProci].size();
+                }
+            }
+
+            labelList& mySub = subMap_[newProci];
+            mySub.resize_nocopy(allSize);
+
+            allSize = 0;
+            for (const label oldProci : oldProcs)
+            {
+                forAll(maps, mapi)
+                {
+                    const auto& map = maps[mapi].subMap()[oldProci];
+                    SubList<label> slice(mySub, map.size(), allSize);
+                    if (subHasFlip_)
+                    {
+                        forAll(map, i)
+                        {
+                            if (map[i] < 0)
+                            {
+                                slice[i] = map[i]-startOfLocal[mapi];
+                            }
+                            else
+                            {
+                                slice[i] = map[i]+startOfLocal[mapi];
+                            }
+                        }
+                    }
+                    else
+                    {
+                        forAll(map, i)
+                        {
+                            slice[i] = map[i]+startOfLocal[mapi];
+                        }
+                    }
+                    allSize += map.size();
+                }
+            }
+        }
+    }
+
+
+    constructSize_ = constructi;
+}
+
+//XXXXXXXX
+/*
+Foam::mapDistributeBase::mapDistributeBase
+(
+    const UPtrList<const mapDistributeBase>& maps,
+    const labelList& localRanks,
+    const label newComm,
+    const labelListList& newToOldRanks,// from newComm to comm_
+    labelList& startOfLocal,
+    List<Map<label>>& compactMaps
+)
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    comm_(-1),
+    schedulePtr_(nullptr)
+{
+    if (maps.empty())
+    {
+        return;
+    }
+
+    comm_ = newComm;
+    subHasFlip_ = maps[0].subHasFlip();
+    constructHasFlip_ = maps[0].constructHasFlip();
+
+    const label nNewRanks = newToOldRanks.size();
+    const label myNewRank = UPstream::myProcNo(newComm);
+    if (nNewRanks != UPstream::nProcs(newComm))
+    {
+        FatalErrorInFunction<< "nNewRanks:" << nNewRanks
+            << " nProcs:" << UPstream::nProcs(newComm)
+            << exit(FatalError);
+    }
+
+    if (localRanks.size() != maps.size())
+    {
+        FatalErrorInFunction
+            << "Number of maps:" << maps.size()
+            << " number of localRanks:" << localRanks.size()
+            << exit(FatalError);
+    }
+
+    // Sanity checks
+    const auto& map0 = maps[0];
+    forAll(maps, mapi)
+    {
+        const auto& map = maps[mapi];
+
+        if
+        (
+            (map.comm() != map0.comm())
+         || (map.subHasFlip() != map0.subHasFlip())
+         || (map.constructHasFlip() != map0.constructHasFlip())
+        )
+        {
+            FatalErrorInFunction
+                << "Maps should all be the same form"
+                << " Map " << mapi
+                << " has comm:" << map.comm()
+                << " subHasFlip:" << map.subHasFlip()
+                << " constructHasFlip:" << map.constructHasFlip()
+                << " which is different from map 0"
+                << exit(FatalError);
+        }
+
+        const label localRank = localRanks[mapi];
+        const auto& constructOwn = maps[mapi].constructMap()[localRank];
+        forAll(constructOwn, i)
+        {
+            if (constructOwn[i] != i)
+            {
+                FatalErrorInFunction
+                    << "Maps constructMap not identity."
+                    << " Map " << mapi
+                    << " constructMap:" << flatOutput(constructOwn)
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    constructMap_.resize_nocopy(nNewRanks);
+    subMap_.resize_nocopy(nNewRanks);
+
+
+    // Store starts. All local data gets appended.
+    startOfLocal.setSize(maps.size()+1);
+    compactMaps.resize_nocopy(maps.size());
+
+    label constructi = 0;
+    forAll(maps, mapi)
+    {
+        startOfLocal[mapi] = constructi;
+
+        const label localRank = localRanks[mapi];
+        const auto& map = maps[mapi].constructMap()[localRank];
+
+        //Pout<< "map:" << mapi << " has localRank:" << localRank
+        //    << " and local size:" << map.size()
+        //    << " and starts at:" << constructi << endl;
+
+        // Presize compaction array
+        const label nRemote = maps[mapi].constructSize()-map.size();
+        compactMaps[mapi].resize(2*nRemote);
+
+        constructi += map.size();
+    }
+    startOfLocal.last() = constructi;
+
+
+    // Determine start of constructed remote data. This is used to get the
+    // local offset which can then be used to get the relative subMap location.
+    labelListList startOfRemote(maps.size());
+    forAll(maps, mapi)
+    {
+        const label nOldProcs = maps[mapi].constructMap().size();
+        labelList& starts = startOfRemote[mapi];
+
+        starts.setSize(nOldProcs, labelMax);
+        forAll(maps[mapi].constructMap(), oldProci)
+        {
+            const labelList& map = maps[mapi].constructMap()[oldProci];
+            forAll(map, i)
+            {
+                const label index
+                (
+                    constructHasFlip_
+                  ? mag(map[i])-1
+                  : map[i]
+                );
+                starts[oldProci] = min(starts[oldProci], index);
+            }
+        }
+    }
+
+
+    // Adjust the submaps. Only the ones going to now remote processors
+    // need to be kept.
+
+    subMap_[myNewRank] = identity(startOfLocal.last());
+    constructMap_[myNewRank] = subMap_[myNewRank];
+
+    // Do (still) remote sub data. Get sent in processor order
+    forAll(newToOldRanks, newProci)
+    {
+        if (newProci != myNewRank)
+        {
+            const auto& oldProcs = newToOldRanks[newProci];
+
+            label allSize = 0;
+            for (const label oldProci : oldProcs)
+            {
+                forAll(maps, mapi)
+                {
+                    const auto& map = maps[mapi].subMap()[oldProci];
+                    allSize += map.size();
+                }
+            }
+
+            labelList& mySub = subMap_[newProci];
+            mySub.resize_nocopy(allSize);
+
+            allSize = 0;
+            forAll(maps, mapi)
+            {
+                for (const label oldProci : oldProcs)
+                {
+                    const auto& map = maps[mapi].subMap()[oldProci];
+
+                    if (!map.size())
+                    {
+                        continue;
+                    }
+
+                    SubList<label> slice(mySub, map.size(), allSize);
+                    forAll(map, i)
+                    {
+                        // old slot (map[i]) becomes new slot by adding the
+                        // offset
+                        if (subHasFlip_ && map[i] < 0)
+                        {
+                            slice[i] = map[i]-startOfLocal[mapi];
+                        }
+                        else
+                        {
+                            slice[i] = map[i]+startOfLocal[mapi];
+                        }
+                    }
+
+                    allSize += map.size();
+                }
+
+                //Pout<< "for map:" << mapi
+                //    << " done all oldProcs:" << oldProcs
+                //    << " and have combined sub:"
+                //    << flatOutput(SubList<label>(mySub, allSize)) << endl;
+            }
+        }
+    }
+
+
+    // Adjust mapping for now local construct data
+    {
+        const auto& oldProcs = newToOldRanks[myNewRank];
+
+        for (const label oldProci : oldProcs)
+        {
+            forAll(maps, mapi)
+            {
+                if (oldProci != localRanks[mapi])
+                {
+                    // Old-local never in compactMap since already in
+                    // startOfLocal
+
+                    const auto& map = maps[mapi].constructMap()[oldProci];
+                    const label sourceMapi = localRanks.find(oldProci);
+                    const auto& subMap = maps[sourceMapi].subMap()[oldProci];
+
+
+                    // So:
+                    // - sourceMapi is sending to oldProci the subMap[oldProci]
+                    //   elements
+                    // - mapi is receiving (in same order) and inserting them
+                    //   into constructMap[oldProci]
+
+                    const label offset = startOfLocal[sourceMapi];
+                    // Construct map starts after the local data
+                    const label nMapLocal = startOfRemote[mapi][oldProci];
+
+                    auto& cptMap = compactMaps[mapi];
+                    forAll(map, i)
+                    {
+                        // old slot position
+                        const label index
+                        (
+                            constructHasFlip_
+                          ? mag(map[i])-1
+                          : map[i]
+                        );
+                        // this was say the nth element in the data received
+                        // from sourceMapi so work out the relative index in
+                        // the constructMap and work out what element was
+                        // received and put it in the correct new location.
+                        const label newIndex = subMap[index-nMapLocal]+offset;
+
+                        // Note: should always warn for duplicates? Or only if
+                        // different?
+                        if
+                        (
+                           !cptMap.insert(index, newIndex)
+                         && cptMap[index] != newIndex
+                        )
+                        {
+                            FatalErrorInFunction<< "Duplicate insertion"
+                                << "From oldProc:" << oldProci
+                                << " on map:" << mapi
+                                << " at index:" << i
+                                << " have construct slot:" << index
+                                << " new index:" << newIndex
+                                << " but already have entry:" << cptMap[index]
+                                << " on for that slot"
+                                << exit(FatalError);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Do (still) remote construct data. (construct-maps are where the received
+    // data gets stored in the flat list). Store mapping from old slot to new
+    // slot.
+    forAll(newToOldRanks, newProci)
+    {
+        if (newProci != myNewRank)
+        {
+            const auto& oldProcs = newToOldRanks[newProci];
+
+            label allSize = 0;
+            for (const label oldProci : oldProcs)
+            {
+                forAll(maps, mapi)
+                {
+                    const auto& map = maps[mapi].constructMap()[oldProci];
+                    allSize += map.size();
+                }
+            }
+
+            labelList& myConstruct = constructMap_[newProci];
+            myConstruct.resize_nocopy(allSize);
+
+            allSize = 0;
+            for (const label oldProci : oldProcs)
+            {
+                forAll(maps, mapi)
+                {
+                    const auto& map = maps[mapi].constructMap()[oldProci];
+
+                    if (!map.size())
+                    {
+                        continue;
+                    }
+
+                    const label nMapLocal = startOfRemote[mapi][oldProci];
+                    SubList<label> slice(myConstruct, map.size(), allSize);
+                    auto& cptMap = compactMaps[mapi];
+                    forAll(map, i)
+                    {
+                        // Construct maps are slots into the old data
+                        const label index
+                        (
+                            constructHasFlip_
+                          ? mag(map[i])-1
+                          : map[i]
+                        );
+                        // Equivalent index into the new data
+                        const label newIndex = (map[i]-nMapLocal) + constructi;
+
+                        slice[i] = newIndex;
+
+                        if
+                        (
+                           !cptMap.insert(index, newIndex)
+                         && cptMap[index] != newIndex
+                        )
+                        {
+                            FatalErrorInFunction<< "Duplicate insertion"
+                                << "From oldProc:" << oldProci
+                                << " on map:" << mapi
+                                << " at index:" << i
+                                << " have construct slot:" << index
+                                << " new index:" << newIndex
+                                << " but already have entry:" << cptMap[index]
+                                << " on for that slot"
+                                << exit(FatalError);
+                        }
+                    }
+
+                    constructi += map.size();
+                    allSize += map.size();
+                }
+            }
+        }
+    }
+    constructSize_ = constructi;
+
+//Pout<< "constructSize:" << constructSize_ << endl;
+//Pout<< "startOfLocal:" << flatOutput(startOfLocal) << endl;
+//Pout<< "compactMaps:" << flatOutput(compactMaps) << endl;
+//Pout<< "subMap:" << endl;
+//forAll(subMap_, newRanki)
+//{
+//    Pout<< "    Sending to " << newRanki << " "
+//        << flatOutput(subMap_[newRanki]) << endl;
+//}
+//Pout<< "constructMap:" << endl;
+//forAll(constructMap_, newRanki)
+//{
+//    Pout<< "    Recieving from " << newRanki << " "
+//        << flatOutput(constructMap_[newRanki]) << endl;
+//}
+}
+*/
+//XXXXXXX
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::labelList Foam::mapDistributeBase::subMapSizes() const
@@ -1088,21 +1902,24 @@ void Foam::mapDistributeBase::transfer(mapDistributeBase& rhs)
 Foam::label Foam::mapDistributeBase::renumber
 (
     const globalIndex& globalNumbering,
+    const label comm,
     const List<Map<label>>& compactMap,
     const label globalI
 )
 {
+    const label myRank = Pstream::myProcNo(comm);
+
     if (globalI == -1)
     {
         return globalI;
     }
-    if (globalNumbering.isLocal(globalI))
+    if (globalNumbering.isLocal(myRank, globalI))
     {
-        return globalNumbering.toLocal(globalI);
+        return globalNumbering.toLocal(myRank, globalI);
     }
     else
     {
-        label proci = globalNumbering.whichProcID(globalI);
+        label proci = globalNumbering.whichProcID(myRank, globalI);
         label index = globalNumbering.toLocal(proci, globalI);
         return compactMap[proci][index];
     }

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBase.H

@@ -484,6 +484,27 @@ public:
             const label comm = UPstream::worldComm
         );
 
+        //- Construct from multiple maps and processor collation
+        //  Assumes all local data first. Sorts contributions of maps
+        //  in processor order i.e.  constructed map has all local data first.
+        //  Returns
+        //  - startOfLocal : per input map the start of the local data. Extends
+        //                   one beyond number of maps so overall local size
+        //                   is startOfLocal.last()
+        //  - compactMaps  : per input map from slot position in the input map
+        //                   to new slot position. (note there is no information
+        //                   returned about which processor it is from)
+        mapDistributeBase
+        (
+            const UPtrList<const mapDistributeBase>& maps,
+            const labelList& localRanks,
+            const label newComm,
+            const labelListList& newToOldRanks, // from rank in newComm to
+                                                // ranks in (old)comm
+            labelList& startOfLocal,            // per map start of local data
+            List<Map<label>>& compactMaps       // per map old slot to new slot
+        );
+
         //- Construct from Istream
         explicit mapDistributeBase(Istream& is);
 
@@ -645,6 +666,7 @@ public:
         static label renumber
         (
             const globalIndex&,
+            const label comm,
             const List<Map<label>>& compactMap,
             const label globalElement
         );
@@ -665,6 +687,36 @@ public:
             const bool hasFlip
         );
 
+        //- Helper for renumbering the map elements. Assumes local elements
+        //- are first, followed by any remote. Local elements get offset,
+        //- remote elements are mapped.
+        //
+        //  \param[in,out] map  The map to be renumbered
+        //  \param localSize    elements < localSize are offset
+        //  \param offset       offset
+        //  \param cMap         map for non-local elements
+        //  \param hasFlip      True if map has flip addressing
+        //
+        //  \return max-size needed for new addressing (eg, constructSize)
+        static label renumberMap
+        (
+            labelList& map,
+            const label localSize,
+            const label offset,
+            const Map<label>& cMap,
+            const bool hasFlip
+        );
+
+        //- Helper for a list of maps. Calls above renumberMap for all elements
+        static label renumberMap
+        (
+            labelListList& mapElements,
+            const label localSize,
+            const label offset,
+            const Map<label>& cMap,
+            const bool hasFlip
+        );
+
 
     // Compaction
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseSubset.C

@@ -329,6 +329,104 @@ Foam::label Foam::mapDistributeBase::renumberMap
 }
 
 
+Foam::label Foam::mapDistributeBase::renumberMap
+(
+    labelList& map,
+    const label localSize,
+    const label offset,
+    const Map<label>& cMap,
+    const bool hasFlip
+)
+{
+    label maxIndex = -1;
+
+    // Transcribe the map
+    if (hasFlip)
+    {
+        for (label& val : map)
+        {
+            // Unflip indexed value
+            const label index = mag(val)-1;
+            if (index < localSize)
+            {
+                // Local element
+                if (val < 0)
+                {
+                    val -= offset;
+                }
+                else
+                {
+                    val += offset;
+                }
+            }
+            else
+            {
+                // Remote element
+                if (val < 0)
+                {
+                    val = -cMap[index]-1;
+                }
+                else
+                {
+                    val = cMap[index]+1;
+                }
+            }
+            maxIndex = max(maxIndex, mag(val)-1);
+        }
+    }
+    else
+    {
+        for (label& val : map)
+        {
+            // Get indexed value (no flipping)
+            if (val < localSize)
+            {
+                val += offset;
+            }
+            else
+            {
+                val = cMap[val];
+            }
+            maxIndex = max(maxIndex, val);
+        }
+    }
+
+    return (maxIndex+1);
+}
+
+
+Foam::label Foam::mapDistributeBase::renumberMap
+(
+    labelListList& mapElements,
+    const label localSize,
+    const label offset,
+    const Map<label>& cMap,
+    const bool hasFlip
+)
+{
+    label maxIndex = -1;
+
+    // Transcribe the map
+    for (labelList& map : mapElements)
+    {
+        maxIndex = max
+        (
+            maxIndex,
+            renumberMap
+            (
+                map,
+                localSize,
+                offset,
+                cMap,
+                hasFlip
+            )
+        );
+    }
+
+    return (maxIndex+1);
+}
+
+
 void Foam::mapDistributeBase::renumberVisitOrder
 (
     const labelUList& origElements,

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -89,7 +89,8 @@ Foam::label Foam::AMIInterpolation::calcDistribution
         (
             UPstream::listGatherValues<bool>
             (
-                srcPatch.size() > 0 || tgtPatch.size() > 0
+                (srcPatch.size() > 0 || tgtPatch.size() > 0),
+                comm_
             )
         );
 
@@ -108,7 +109,7 @@ Foam::label Foam::AMIInterpolation::calcDistribution
                 << "AMI split across multiple processors" << endl;
         }
 
-        Pstream::broadcast(proci);
+        Pstream::broadcast(proci, comm_);
     }
 
     return proci;
@@ -164,7 +165,8 @@ void Foam::AMIInterpolation::normaliseWeights
     scalarField& wghtSum,
     const bool conformal,
     const bool output,
-    const scalar lowWeightTol
+    const scalar lowWeightTol,
+    const label comm
 )
 {
     addProfiling(ami, "AMIInterpolation::normaliseWeights");
@@ -207,6 +209,11 @@ void Foam::AMIInterpolation::normaliseWeights
 
     if (output)
     {
+        // Note: change global communicator since gMin,gAverage etc don't
+        // support user communicator
+        const label oldWorldComm(UPstream::worldComm);
+        UPstream::worldComm = comm;
+
         if (returnReduceOr(wght.size()))
         {
             Info<< indent
@@ -227,6 +234,8 @@ void Foam::AMIInterpolation::normaliseWeights
                     << endl;
             }
         }
+
+        UPstream::worldComm = oldWorldComm;
     }
 }
 
@@ -245,7 +254,8 @@ void Foam::AMIInterpolation::agglomerate
     labelListList& srcAddress,
     scalarListList& srcWeights,
     scalarField& srcWeightsSum,
-    autoPtr<mapDistribute>& tgtMap
+    autoPtr<mapDistribute>& tgtMap,
+    const label comm
 )
 {
     addProfiling(ami, "AMIInterpolation::agglomerate");
@@ -266,7 +276,9 @@ void Foam::AMIInterpolation::agglomerate
 
     // Agglomerate face areas
     {
-        srcMagSf.setSize(sourceRestrictAddressing.size(), 0.0);
+        //srcMagSf.setSize(sourceRestrictAddressing.size(), 0.0);
+        srcMagSf.setSize(sourceCoarseSize, 0.0);
+
         forAll(sourceRestrictAddressing, facei)
         {
             label coarseFacei = sourceRestrictAddressing[facei];
@@ -309,16 +321,16 @@ void Foam::AMIInterpolation::agglomerate
         // - a subMap : these are face indices
         // - a constructMap : these are from 'transferred-data' to slots
 
-        labelListList tgtSubMap(Pstream::nProcs());
+        labelListList tgtSubMap(Pstream::nProcs(comm));
 
         // Local subMap is just identity
         {
-            tgtSubMap[Pstream::myProcNo()] = identity(targetCoarseSize);
+            tgtSubMap[Pstream::myProcNo(comm)] = identity(targetCoarseSize);
         }
 
         forAll(map.subMap(), proci)
         {
-            if (proci != Pstream::myProcNo())
+            if (proci != Pstream::myProcNo(comm))
             {
                 // Combine entries that point to the same coarse element.
                 // The important bit is to loop over the data (and hand out
@@ -330,7 +342,7 @@ void Foam::AMIInterpolation::agglomerate
 
                 const labelList& elems = map.subMap()[proci];
                 const labelList& elemsMap =
-                    map.constructMap()[Pstream::myProcNo()];
+                    map.constructMap()[Pstream::myProcNo(comm)];
                 labelList& newSubMap = tgtSubMap[proci];
                 newSubMap.resize_nocopy(elems.size());
 
@@ -356,11 +368,12 @@ void Foam::AMIInterpolation::agglomerate
         // of handing out indices should be the same as loop above to compact
         // the sending map
 
-        labelListList tgtConstructMap(Pstream::nProcs());
+        labelListList tgtConstructMap(Pstream::nProcs(comm));
 
         // Local constructMap is just identity
         {
-            tgtConstructMap[Pstream::myProcNo()] = identity(targetCoarseSize);
+            tgtConstructMap[Pstream::myProcNo(comm)] =
+                identity(targetCoarseSize);
         }
 
         labelList tgtCompactMap(map.constructSize());
@@ -372,7 +385,8 @@ void Foam::AMIInterpolation::agglomerate
             // 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()];
+            const labelList& elemsMap =
+                map.constructMap()[Pstream::myProcNo(comm)];
             for (const label fineElem : elemsMap)
             {
                 label coarseElem = allRestrict[fineElem];
@@ -385,7 +399,7 @@ void Foam::AMIInterpolation::agglomerate
         // Compact data from other processors
         forAll(map.constructMap(), proci)
         {
-            if (proci != Pstream::myProcNo())
+            if (proci != Pstream::myProcNo(comm))
             {
                 // Combine entries that point to the same coarse element. All
                 // elements now are remote data so we cannot use any local
@@ -477,7 +491,10 @@ void Foam::AMIInterpolation::agglomerate
             (
                 compacti,
                 std::move(tgtSubMap),
-                std::move(tgtConstructMap)
+                std::move(tgtConstructMap),
+                false,      //subHasFlip
+                false,      //constructHasFlip
+                comm
             )
         );
     }
@@ -528,7 +545,8 @@ void Foam::AMIInterpolation::agglomerate
         srcWeightsSum,
         true,
         false,
-        -1
+        -1,
+        comm
     );
 }
 
@@ -545,6 +563,7 @@ Foam::AMIInterpolation::AMIInterpolation
     reverseTarget_(dict.getOrDefault("reverseTarget", reverseTarget)),
     lowWeightCorrection_(dict.getOrDefault<scalar>("lowWeightCorrection", -1)),
     singlePatchProc_(-999),
+    comm_(UPstream::worldComm),
     srcMagSf_(),
     srcAddress_(),
     srcWeights_(),
@@ -572,6 +591,7 @@ Foam::AMIInterpolation::AMIInterpolation
     reverseTarget_(reverseTarget),
     lowWeightCorrection_(lowWeightCorrection),
     singlePatchProc_(-999),
+    comm_(UPstream::worldComm),
     srcMagSf_(),
     srcAddress_(),
     srcWeights_(),
@@ -601,6 +621,7 @@ Foam::AMIInterpolation::AMIInterpolation
     reverseTarget_(fineAMI.reverseTarget_),
     lowWeightCorrection_(-1.0),
     singlePatchProc_(fineAMI.singlePatchProc_),
+    comm_(fineAMI.comm_),
     srcMagSf_(),
     srcAddress_(),
     srcWeights_(),
@@ -634,6 +655,7 @@ Foam::AMIInterpolation::AMIInterpolation
         Pout<< "AMI: Creating addressing and weights as agglomeration of AMI :"
             << " source:" << fineAMI.srcAddress().size()
             << " target:" << fineAMI.tgtAddress().size()
+            << " fineComm:" << fineAMI.comm()
             << " coarse source size:" << sourceCoarseSize
             << " neighbour source size:" << neighbourCoarseSize
             << endl;
@@ -673,7 +695,8 @@ Foam::AMIInterpolation::AMIInterpolation
         srcAddress_,
         srcWeights_,
         srcWeightsSum_,
-        tgtMapPtr_
+        tgtMapPtr_,
+        comm_
     );
 
     agglomerate
@@ -690,7 +713,8 @@ Foam::AMIInterpolation::AMIInterpolation
         tgtAddress_,
         tgtWeights_,
         tgtWeightsSum_,
-        srcMapPtr_
+        srcMapPtr_,
+        comm_
     );
 }
 
@@ -701,6 +725,7 @@ Foam::AMIInterpolation::AMIInterpolation(const AMIInterpolation& ami)
     reverseTarget_(ami.reverseTarget_),
     lowWeightCorrection_(ami.lowWeightCorrection_),
     singlePatchProc_(ami.singlePatchProc_),
+    comm_(ami.comm_),
     srcMagSf_(ami.srcMagSf_),
     srcAddress_(ami.srcAddress_),
     srcWeights_(ami.srcWeights_),
@@ -717,6 +742,40 @@ Foam::AMIInterpolation::AMIInterpolation(const AMIInterpolation& ami)
 {}
 
 
+Foam::AMIInterpolation::AMIInterpolation(Istream& is)
+:
+    requireMatch_(readBool(is)),
+    reverseTarget_(readBool(is)),
+    lowWeightCorrection_(readScalar(is)),
+    singlePatchProc_(readLabel(is)),
+    comm_(readLabel(is)),
+
+    srcMagSf_(is),
+    srcAddress_(is),
+    srcWeights_(is),
+    srcWeightsSum_(is),
+    srcCentroids_(is),
+    //srcPatchPts_(is),
+    srcMapPtr_(nullptr),
+
+    tgtMagSf_(is),
+    tgtAddress_(is),
+    tgtWeights_(is),
+    tgtWeightsSum_(is),
+    tgtCentroids_(is),
+    //tgtPatchPts_(is),
+    tgtMapPtr_(nullptr),
+
+    upToDate_(readBool(is))
+{
+    if (singlePatchProc_ == -1)
+    {
+        srcMapPtr_.reset(new mapDistribute(is));
+        tgtMapPtr_.reset(new mapDistribute(is));
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::AMIInterpolation::calculate
@@ -757,8 +816,20 @@ bool Foam::AMIInterpolation::calculate
         ttgtPatch0_.cref(tgtPatch);
     }
 
-    label srcTotalSize = returnReduce(srcPatch.size(), sumOp<label>());
-    label tgtTotalSize = returnReduce(tgtPatch.size(), sumOp<label>());
+    label srcTotalSize = returnReduce
+    (
+        srcPatch.size(),
+        sumOp<label>(),
+        UPstream::msgType(),
+        comm_
+    );
+    label tgtTotalSize = returnReduce
+    (
+        tgtPatch.size(),
+        sumOp<label>(),
+        UPstream::msgType(),
+        comm_
+    );
 
     if (srcTotalSize == 0)
     {
@@ -794,7 +865,8 @@ void Foam::AMIInterpolation::reset
     labelListList&& srcAddress,
     scalarListList&& srcWeights,
     labelListList&& tgtAddress,
-    scalarListList&& tgtWeights
+    scalarListList&& tgtWeights,
+    const label singlePatchProc
 )
 {
     DebugInFunction<< endl;
@@ -820,6 +892,8 @@ void Foam::AMIInterpolation::reset
     srcMapPtr_ = std::move(srcToTgtMap);
     tgtMapPtr_ = std::move(tgtToSrcMap);
 
+    singlePatchProc_ = singlePatchProc;
+
     upToDate_ = true;
 }
 
@@ -1030,7 +1104,8 @@ void Foam::AMIInterpolation::normaliseWeights
         srcWeightsSum_,
         conformal,
         output,
-        lowWeightCorrection_
+        lowWeightCorrection_,
+        comm_
     );
 
     normaliseWeights
@@ -1042,7 +1117,8 @@ void Foam::AMIInterpolation::normaliseWeights
         tgtWeightsSum_,
         conformal,
         output,
-        lowWeightCorrection_
+        lowWeightCorrection_,
+        comm_
     );
 }
 
@@ -1261,4 +1337,36 @@ void Foam::AMIInterpolation::write(Ostream& os) const
 }
 
 
+bool Foam::AMIInterpolation::writeData(Ostream& os) const
+{
+    os  << requireMatch()
+        << token::SPACE<< reverseTarget()
+        << token::SPACE<< lowWeightCorrection()
+        << token::SPACE<< singlePatchProc()
+        << token::SPACE<< comm()
+
+        << token::SPACE<< srcMagSf()
+        << token::SPACE<< srcAddress()
+        << token::SPACE<< srcWeights()
+        << token::SPACE<< srcWeightsSum()
+        << token::SPACE<< srcCentroids()
+
+        << token::SPACE<< tgtMagSf()
+        << token::SPACE<< tgtAddress()
+        << token::SPACE<< tgtWeights()
+        << token::SPACE<< tgtWeightsSum()
+        << token::SPACE<< tgtCentroids_
+
+        << token::SPACE<< upToDate();
+
+    if (distributed())
+    {
+        os  << token::SPACE<< srcMap()
+            << token::SPACE<< tgtMap();
+    }
+
+    return os.good();
+}
+
+
 // ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -52,20 +52,14 @@ SourceFiles
 #ifndef Foam_AMIInterpolation_H
 #define Foam_AMIInterpolation_H
 
-#include "className.H"
 #include "searchableSurface.H"
-#include "treeBoundBoxList.H"
 #include "boolList.H"
 #include "primitivePatch.H"
-#include "faceAreaIntersect.H"
-#include "globalIndex.H"
 #include "ops.H"
 #include "refPtr.H"
-#include "Enum.H"
 #include "pointList.H"
 #include "indexedOctree.H"
 #include "treeDataPrimitivePatch.H"
-
 #include "runTimeSelectionTables.H"
 
 
@@ -109,6 +103,9 @@ protected:
         //- The value is -1 for distributed cases
         label singlePatchProc_;
 
+        //- Communicator to use for parallel operations.
+        label comm_;
+
 
         // Source patch
 
@@ -224,7 +221,8 @@ protected:
                 scalarField& wghtSum,
                 const bool conformal,
                 const bool output,
-                const scalar lowWeightTol
+                const scalar lowWeightTol,
+                const label comm
             );
 
 
@@ -244,7 +242,8 @@ protected:
                 labelListList& srcAddress,
                 scalarListList& srcWeights,
                 scalarField& srcWeightsSum,
-                autoPtr<mapDistribute>& tgtMap
+                autoPtr<mapDistribute>& tgtMap,
+                const label comm
             );
 
 
@@ -342,6 +341,10 @@ public:
             return autoPtr<AMIInterpolation>::New(*this);
         }
 
+        //- Construct from Istream
+        //  Note: does not transfer src(tgt)PatchPts, tsrc(tgt)Patch0_
+        AMIInterpolation(Istream& is);
+
 
     //- Destructor
     virtual ~AMIInterpolation() = default;
@@ -371,7 +374,7 @@ public:
             //- Set the require match flag, return the \em new value
             inline bool setRequireMatch(bool flag) noexcept;
 
-            //- Return true if requireMatch and lowWeightCorrectionin active
+            //- Return true if requireMatch and but not lowWeightCorrection
             inline bool mustMatchFaces() const;
 
             //- Access to the reverseTarget flag
@@ -387,6 +390,13 @@ public:
             //- or -1 if distributed
             inline label singlePatchProc() const noexcept;
 
+            //- Communicator to use for parallel operations
+            inline label comm() const;
+
+            //- Set communicator to use for parallel operations. Return
+            //  old value
+            inline label comm(const label newComm);
+
 
             // Source patch
 
@@ -488,7 +498,8 @@ public:
                 labelListList&& srcAddress,
                 scalarListList&& srcWeights,
                 labelListList&& tgtAddress,
-                scalarListList&& tgtWeights
+                scalarListList&& tgtWeights,
+                const label singlePatchProc
             );
 
             //- Append additional addressing and weights
@@ -664,9 +675,12 @@ public:
 
         // I-O
 
-            //- Write
+            //- Write AMI as a dictionary
             virtual void write(Ostream& os) const;
 
+            //- Write AMI raw
+            virtual bool writeData(Ostream& os) const;
+
 
     // Housekeeping
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H

@@ -102,6 +102,20 @@ inline Foam::label Foam::AMIInterpolation::singlePatchProc() const noexcept
 }
 
 
+inline Foam::label Foam::AMIInterpolation::comm() const
+{
+    return comm_;
+}
+
+
+inline Foam::label Foam::AMIInterpolation::comm(const label newComm)
+{
+    const label oldComm(comm_);
+    comm_ = newComm;
+    return oldComm;
+}
+
+
 inline const Foam::List<Foam::scalar>& Foam::AMIInterpolation::srcMagSf() const
 {
     return srcMagSf_;

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C

@@ -68,8 +68,36 @@ void Foam::advancingFrontAMI::checkPatches() const
         const scalar maxBoundsError = 0.05;
 
         // Check bounds of source and target
-        boundBox bbSrc(src.points(), src.meshPoints(), true);
-        boundBox bbTgt(tgt.points(), tgt.meshPoints(), true);
+        boundBox bbSrc(src.points(), src.meshPoints(), false);
+        Foam::reduce
+        (
+            bbSrc.min(),
+            minOp<point>(),
+            UPstream::msgType(),
+            comm_
+        );
+        Foam::reduce
+        (
+            bbSrc.max(),
+            maxOp<point>(),
+            UPstream::msgType(),
+            comm_
+        );
+        boundBox bbTgt(tgt.points(), tgt.meshPoints(), false);
+        Foam::reduce
+        (
+            bbTgt.min(),
+            minOp<point>(),
+            UPstream::msgType(),
+            comm_
+        );
+        Foam::reduce
+        (
+            bbTgt.max(),
+            maxOp<point>(),
+            UPstream::msgType(),
+            comm_
+        );
 
         boundBox bbTgtInf(bbTgt);
         bbTgtInf.inflate(maxBoundsError);
@@ -150,7 +178,7 @@ void Foam::advancingFrontAMI::createExtendedTgtPatch()
 
     // Original faces from tgtPatch
     // Note: in globalIndexing since might be remote
-    globalIndex globalTgtFaces(tgtPatch0().size());
+    globalIndex globalTgtFaces(tgtPatch0().size(), comm_);
     distributeAndMergePatches
     (
         map,

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C

@@ -74,14 +74,18 @@ void Foam::advancingFrontAMI::distributePatches
     List<labelList>& faceIDs
 ) const
 {
-    faces.resize_nocopy(Pstream::nProcs());
-    points.resize_nocopy(Pstream::nProcs());
-    faceIDs.resize_nocopy(Pstream::nProcs());
+    faces.setSize(Pstream::nProcs(comm_));
+    points.setSize(Pstream::nProcs(comm_));
+    faceIDs.setSize(Pstream::nProcs(comm_));
 
+    PstreamBuffers pBufs
+    (
+        Pstream::commsTypes::nonBlocking,
+        UPstream::msgType(),
+        comm_
+    );
 
-    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
-
-    for (const int domain : Pstream::allProcs())
+    for (const int domain : Pstream::allProcs(comm_))
     {
         const labelList& sendElems = map.subMap()[domain];
 
@@ -104,12 +108,13 @@ void Foam::advancingFrontAMI::distributePatches
             }
 
 
-            if (domain == Pstream::myProcNo())
+            if (domain == Pstream::myProcNo(comm_))
             {
                 // Do send/receive for myself
                 faces[domain] = subPatch.localFaces();
                 points[domain] = subPatch.localPoints();
-                faceIDs[domain] = gi.toGlobal(sendElems);
+                faceIDs[domain] =
+                    gi.toGlobal(Pstream::myProcNo(comm_), sendElems);
             }
             else
             {
@@ -118,19 +123,20 @@ void Foam::advancingFrontAMI::distributePatches
                 str
                     << subPatch.localFaces()
                     << subPatch.localPoints()
-                    << gi.toGlobal(sendElems);
+                    << gi.toGlobal(Pstream::myProcNo(comm_), sendElems);
             }
         }
     }
 
     pBufs.finishedSends();
 
+
     // Consume
-    for (const int domain : Pstream::allProcs())
+    for (const int domain : Pstream::allProcs(comm_))
     {
         const labelList& recvElems = map.constructMap()[domain];
 
-        if (domain != Pstream::myProcNo() && recvElems.size())
+        if (domain != Pstream::myProcNo(comm_) && recvElems.size())
         {
             UIPstream is(domain, pBufs);
 
@@ -176,10 +182,10 @@ void Foam::advancingFrontAMI::distributeAndMergePatches
 
     // My own data first
     {
-        const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
+        const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo(comm_)];
         SubList<label>(tgtFaceIDs, faceIDs.size()) = faceIDs;
 
-        const faceList& fcs = allFaces[Pstream::myProcNo()];
+        const faceList& fcs = allFaces[Pstream::myProcNo(comm_)];
         for (const face& f : fcs)
         {
             face& newF = tgtFaces[nFaces++];
@@ -190,7 +196,7 @@ void Foam::advancingFrontAMI::distributeAndMergePatches
             }
         }
 
-        const pointField& pts = allPoints[Pstream::myProcNo()];
+        const pointField& pts = allPoints[Pstream::myProcNo(comm_)];
         for (const point& pt: pts)
         {
             tgtPoints[nPoints++] = pt;
@@ -201,7 +207,7 @@ void Foam::advancingFrontAMI::distributeAndMergePatches
     // Other proc data follows
     forAll(allFaces, proci)
     {
-        if (proci != Pstream::myProcNo())
+        if (proci != Pstream::myProcNo(comm_))
         {
             const labelList& faceIDs = allTgtFaceIDs[proci];
             SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces) = faceIDs;
@@ -258,11 +264,11 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
 ) const
 {
     // Get decomposition of patch
-    List<treeBoundBoxList> procBb(Pstream::nProcs());
+    List<treeBoundBoxList> procBb(Pstream::nProcs(comm_));
 
     if (srcPatch.size())
     {
-        procBb[Pstream::myProcNo()] =
+        procBb[Pstream::myProcNo(comm_)] =
             AABBTree<face>
             (
                 srcPatch.localFaces(),
@@ -272,10 +278,10 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
     }
     else
     {
-        procBb[Pstream::myProcNo()] = treeBoundBoxList();
+        procBb[Pstream::myProcNo(comm_)] = treeBoundBoxList();
     }
 
-    Pstream::allGatherList(procBb);
+    Pstream::allGatherList(procBb, UPstream::msgType(), comm_);
 
     if (debug)
     {
@@ -295,10 +301,10 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
 
     {
         // Per processor indices into all segments to send
-        List<DynamicList<label>> dynSendMap(Pstream::nProcs());
+        List<DynamicList<label>> dynSendMap(Pstream::nProcs(comm_));
 
         // Work array - whether processor bb overlaps the face bounds
-        boolList procBbOverlaps(Pstream::nProcs());
+        boolList procBbOverlaps(Pstream::nProcs(comm_));
 
         forAll(faces, facei)
         {
@@ -320,7 +326,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
         }
 
         // Convert dynamicList to labelList
-        sendMap.setSize(Pstream::nProcs());
+        sendMap.setSize(Pstream::nProcs(comm_));
         forAll(sendMap, proci)
         {
             sendMap[proci].transfer(dynSendMap[proci]);
@@ -338,7 +344,13 @@ Foam::autoPtr<Foam::mapDistribute> Foam::advancingFrontAMI::calcProcMap
         }
     }
 
-    return autoPtr<mapDistribute>::New(std::move(sendMap));
+    return autoPtr<mapDistribute>::New
+    (
+        std::move(sendMap),
+        false,      //subHasFlip
+        false,      //constructHasFlip
+        comm_
+    );
 }
 
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C

@@ -713,12 +713,14 @@ bool Foam::faceAreaWeightAMI::calculate
 
     if (distributed())
     {
+        const label myRank = UPstream::myProcNo(comm_);
+
         const primitivePatch& srcPatch0 = this->srcPatch0();
         const primitivePatch& tgtPatch0 = this->tgtPatch0();
 
         // Create global indexing for each original patch
-        globalIndex globalSrcFaces(srcPatch0.size());
-        globalIndex globalTgtFaces(tgtPatch0.size());
+        globalIndex globalSrcFaces(srcPatch0.size(), comm_);
+        globalIndex globalTgtFaces(tgtPatch0.size(), comm_);
 
         for (labelList& addressing : srcAddress_)
         {
@@ -730,7 +732,7 @@ bool Foam::faceAreaWeightAMI::calculate
 
         for (labelList& addressing : tgtAddress_)
         {
-            globalSrcFaces.inplaceToGlobal(addressing);
+            globalSrcFaces.inplaceToGlobal(myRank, addressing);
         }
 
         // Send data back to originating procs. Note that contributions
@@ -748,7 +750,9 @@ bool Foam::faceAreaWeightAMI::calculate
             tgtAddress_,
             labelList(),
             ListOps::appendEqOp<label>(),
-            flipOp()                    // flip operation
+            flipOp(),                   // flip operation
+            UPstream::msgType(),
+            comm_
         );
 
         mapDistributeBase::distribute
@@ -763,7 +767,9 @@ bool Foam::faceAreaWeightAMI::calculate
             tgtWeights_,
             scalarList(),
             ListOps::appendEqOp<scalar>(),
-            flipOp()
+            flipOp(),
+            UPstream::msgType(),
+            comm_
         );
 
         // Note: using patch face areas calculated by the AMI method
@@ -773,13 +779,27 @@ bool Foam::faceAreaWeightAMI::calculate
         List<Map<label>> cMapSrc;
         srcMapPtr_.reset
         (
-            new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
+            new mapDistribute
+            (
+                globalSrcFaces,
+                tgtAddress_,
+                cMapSrc,
+                UPstream::msgType(),
+                comm_
+            )
         );
 
         List<Map<label>> cMapTgt;
         tgtMapPtr_.reset
         (
-            new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
+            new mapDistribute
+            (
+                globalTgtFaces,
+                srcAddress_,
+                cMapTgt,
+                UPstream::msgType(),
+                comm_
+            )
         );
     }
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C

@@ -427,12 +427,14 @@ bool Foam::faceAreaWeightAMI2D::calculate
 
     if (distributed())
     {
+        const label myRank = UPstream::myProcNo(comm_);
+
         const primitivePatch& srcPatch0 = this->srcPatch0();
         const primitivePatch& tgtPatch0 = this->tgtPatch0();
 
         // Create global indexing for each original patch
-        globalIndex globalSrcFaces(srcPatch0.size());
-        globalIndex globalTgtFaces(tgtPatch0.size());
+        const globalIndex globalSrcFaces(srcPatch0.size(), comm_);
+        const globalIndex globalTgtFaces(tgtPatch0.size(), comm_);
 
         for (labelList& addressing : srcAddress_)
         {
@@ -444,7 +446,7 @@ bool Foam::faceAreaWeightAMI2D::calculate
 
         for (labelList& addressing : tgtAddress_)
         {
-            globalSrcFaces.inplaceToGlobal(addressing);
+            globalSrcFaces.inplaceToGlobal(myRank, addressing);
         }
 
         // Send data back to originating procs. Note that contributions
@@ -462,7 +464,9 @@ bool Foam::faceAreaWeightAMI2D::calculate
             tgtAddress_,
             labelList(),
             ListOps::appendEqOp<label>(),
-            flipOp()
+            flipOp(),                   // flip operation
+            UPstream::msgType(),
+            comm_
         );
 
         mapDistributeBase::distribute
@@ -477,7 +481,9 @@ bool Foam::faceAreaWeightAMI2D::calculate
             tgtWeights_,
             scalarList(),
             ListOps::appendEqOp<scalar>(),
-            flipOp()
+            flipOp(),                   // flip operation
+            UPstream::msgType(),
+            comm_
         );
 
         // Note: using patch face areas calculated by the AMI method
@@ -487,13 +493,27 @@ bool Foam::faceAreaWeightAMI2D::calculate
         List<Map<label>> cMapSrc;
         srcMapPtr_.reset
         (
-            new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
+            new mapDistribute
+            (
+                globalSrcFaces,
+                tgtAddress_,
+                cMapSrc,
+                UPstream::msgType(),
+                comm_
+            )
         );
 
         List<Map<label>> cMapTgt;
         tgtMapPtr_.reset
         (
-            new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
+            new mapDistribute
+            (
+                globalTgtFaces,
+                srcAddress_,
+                cMapTgt,
+                UPstream::msgType(),
+                comm_
+            )
         );
     }
 

src/meshTools/AMIInterpolation/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C

@@ -111,6 +111,9 @@ Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
     globalIndex globalTgtCells(tgt.size());
 
 
+    const label myRank = UPstream::myProcNo(comm_);
+
+
     // First pass
     // ==========
     // For each srcPatch face, determine local match on tgtPatch
@@ -134,7 +137,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
             {
                 // With a search radius2 of GREAT all cells should receive a hit
                 localInfo[srcCelli].second() = test.point().distSqr(srcCc);
-                test.setIndex(globalTgtCells.toGlobal(test.index()));
+                test.setIndex(globalTgtCells.toGlobal(myRank, test.index()));
             }
         }
     }
@@ -176,7 +179,7 @@ Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
             test = tgtTree.findNearest(srcCcs[i], remoteInfo[i].second());
             if (test.hit())
             {
-                test.setIndex(globalTgtCells.toGlobal(test.index()));
+                test.setIndex(globalTgtCells.toGlobal(myRank, test.index()));
                 testInfo.second() = test.point().distSqr(srcCcs[i]);
                 nearestEqOp()(remoteInfo[i], testInfo);
             }
@@ -200,7 +203,9 @@ Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
         remoteInfo,
         nearestZero,
         nearestEqOp(),
-        identityOp()  // No flipping
+        identityOp(),       // No flipping
+        UPstream::msgType(),
+        comm_
     );
 
 
@@ -223,7 +228,14 @@ Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
     }
 
     List<Map<label>> cMap;
-    return autoPtr<mapDistribute>::New(globalTgtCells, srcToTgtAddr, cMap);
+    return autoPtr<mapDistribute>::New
+    (
+        globalTgtCells,
+        srcToTgtAddr,
+        cMap,
+        UPstream::msgType(),
+        comm_
+    );
 }
 
 

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.C

@@ -47,6 +47,12 @@ namespace Foam
         cyclicACMIGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        cyclicACMIGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -85,6 +91,38 @@ Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
 {}
 
 
+Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
+Foam::cyclicACMIGAMGInterfaceField::cyclicACMIGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const lduInterfaceField& local,
+    const UPtrList<lduInterfaceField>& other
+)
+:
+    GAMGInterfaceField(GAMGCp, local),
+    cyclicACMIInterface_(refCast<const cyclicACMIGAMGInterface>(GAMGCp)),
+    doTransform_(false),
+    rank_(0)
+{
+    const auto& p = refCast<const cyclicACMILduInterfaceField>(local);
+
+    doTransform_ = p.doTransform();
+    rank_ = p.rank();
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cyclicACMIGAMGInterfaceField::~cyclicACMIGAMGInterfaceField()
@@ -132,4 +170,12 @@ void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
 }
 
 
+void Foam::cyclicACMIGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
+}
+
+
 // ************************************************************************* //

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.H

@@ -104,6 +104,39 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        cyclicACMIGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct from GAMG interface and local and remote fields
+        cyclicACMIGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const lduInterfaceField& local,
+            const UPtrList<lduInterfaceField>& other
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            return autoPtr<GAMGInterfaceField>
+            (
+                new cyclicACMIGAMGInterfaceField
+                (
+                    GAMGCp,
+                    *this,      // local field
+                    other       // other fields
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicACMIGAMGInterfaceField();
@@ -161,6 +194,12 @@ public:
             {
                 return rank_;
             }
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,6 +47,12 @@ namespace Foam
         cyclicAMIGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        cyclicAMIGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -85,6 +91,38 @@ Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
 {}
 
 
+Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
+Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const lduInterfaceField& local,
+    const UPtrList<lduInterfaceField>& other
+)
+:
+    GAMGInterfaceField(GAMGCp, local),
+    cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
+    doTransform_(false),
+    rank_(0)
+{
+    const auto& p = refCast<const cyclicAMILduInterfaceField>(local);
+
+    doTransform_ = p.doTransform();
+    rank_ = p.rank();
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cyclicAMIGAMGInterfaceField::~cyclicAMIGAMGInterfaceField()
@@ -107,6 +145,8 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
 {
     // Get neighbouring field
 
+    const label oldWarnComm = UPstream::warnComm;
+
     const labelList& nbrFaceCells =
         lduAddr.patchAddr
         (
@@ -120,16 +160,38 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
 
     if (cyclicAMIInterface_.owner())
     {
-        pnf = cyclicAMIInterface_.AMI().interpolateToSource(pnf);
+        const auto& AMI = cyclicAMIInterface_.AMI();
+
+        // Switch on warning if using wrong communicator. Can be removed if
+        // sure all is correct
+        UPstream::warnComm = AMI.comm();
+
+        pnf = AMI.interpolateToSource(pnf);
     }
     else
     {
-        pnf = cyclicAMIInterface_.neighbPatch().AMI().interpolateToTarget(pnf);
+        const auto& AMI = cyclicAMIInterface_.neighbPatch().AMI();
+
+        // Switch on warning if using wrong communicator. Can be removed if
+        // sure all is correct
+        UPstream::warnComm = AMI.comm();
+
+        pnf = AMI.interpolateToTarget(pnf);
     }
 
     const labelUList& faceCells = lduAddr.patchAddr(patchId);
 
     this->addToInternalField(result, !add, faceCells, coeffs, pnf);
+
+    UPstream::warnComm = oldWarnComm;
+}
+
+
+void Foam::cyclicAMIGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
 }
 
 

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.H

@@ -103,6 +103,39 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        cyclicAMIGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct from GAMG interface and local and remote fields
+        cyclicAMIGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const lduInterfaceField& local,
+            const UPtrList<lduInterfaceField>& other
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            return autoPtr<GAMGInterfaceField>
+            (
+                new cyclicAMIGAMGInterfaceField
+                (
+                    GAMGCp,
+                    *this,      // local field
+                    other       // other fields
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicAMIGAMGInterfaceField();
@@ -160,6 +193,12 @@ public:
             {
                 return rank_;
             }
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,12 @@ namespace Foam
         cyclicACMIGAMGInterface,
         lduInterface
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterface,
+        cyclicACMIGAMGInterface,
+        Istream
+    );
 }
 
 
@@ -63,11 +69,26 @@ Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
         index,
         coarseInterfaces
     ),
-    fineCyclicACMIInterface_
+    neighbPatchID_
     (
-        refCast<const cyclicACMILduInterface>(fineInterface)
+        refCast<const cyclicACMILduInterface>(fineInterface).neighbPatchID()
+    ),
+    owner_
+    (
+        refCast<const cyclicACMILduInterface>(fineInterface).owner()
+    ),
+    forwardT_
+    (
+        refCast<const cyclicACMILduInterface>(fineInterface).forwardT()
+    ),
+    reverseT_
+    (
+        refCast<const cyclicACMILduInterface>(fineInterface).reverseT()
     )
 {
+    const auto& fineCyclicACMIInterface =
+        refCast<const cyclicACMILduInterface>(fineInterface);
+
     // Construct face agglomeration from cell agglomeration
     {
         // From coarse face to cell
@@ -107,7 +128,7 @@ Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
 
     // On the owner side construct the AMI
 
-    if (fineCyclicACMIInterface_.owner())
+    if (fineCyclicACMIInterface.owner())
     {
         // Construct the neighbour side agglomeration (as the neighbour would
         // do it so it the exact loop above using neighbourRestrictAddressing
@@ -144,11 +165,12 @@ Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
             nbrFaceRestrictAddressing.transfer(dynNbrFaceRestrictAddressing);
         }
 
+
         amiPtr_.reset
         (
             new AMIPatchToPatchInterpolation
             (
-                fineCyclicACMIInterface_.AMI(),
+                fineCyclicACMIInterface.AMI(),
                 faceRestrictAddressing_,
                 nbrFaceRestrictAddressing
             )
@@ -157,10 +179,520 @@ Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    Istream& is
+)
+:
+    GAMGInterface(index, coarseInterfaces, is),
+    neighbPatchID_(readLabel(is)),
+    owner_(readBool(is)),
+    forwardT_(is),
+    reverseT_(is)
+{
+    const bool hasAMI(readBool(is));
 
-Foam::cyclicACMIGAMGInterface::~cyclicACMIGAMGInterface()
-{}
+    if (hasAMI)
+    {
+        amiPtr_.reset(new AMIPatchToPatchInterpolation(is));
+    }
+}
+
+
+Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    const lduInterface& fineInterface,
+    const labelList& interfaceMap,
+    const labelUList& faceCells,
+    const labelUList& faceRestrictAddresssing,
+    const labelUList& faceOffsets,
+    const lduInterfacePtrsList& allInterfaces,
+    const label coarseComm,
+    const label myProcNo,
+    const labelList& procAgglomMap
+)
+:
+    GAMGInterface
+    (
+        index,
+        coarseInterfaces,
+        faceCells,
+        faceRestrictAddresssing
+    ),
+    neighbPatchID_
+    (
+        interfaceMap.find
+        (
+            refCast
+            <
+                const cyclicACMILduInterface
+            >(fineInterface).neighbPatchID()
+        )
+    ),
+    owner_
+    (
+        refCast<const cyclicACMILduInterface>(fineInterface).owner()
+    ),
+    forwardT_
+    (
+        refCast<const cyclicACMILduInterface>(fineInterface).forwardT()
+    ),
+    reverseT_
+    (
+       refCast<const cyclicACMILduInterface>(fineInterface).reverseT()
+    )
+{
+    const auto& fineCyclicACMIInterface =
+        refCast<const cyclicACMIGAMGInterface>(fineInterface);
+
+    if (fineCyclicACMIInterface.amiPtr_)
+    {
+        const auto& AMI = const_cast<AMIPatchToPatchInterpolation&>
+        (
+            fineCyclicACMIInterface.AMI()
+        );
+
+        label singlePatchProc = AMI.singlePatchProc();
+
+
+        // Get some sizes
+        label nSrc = 0;
+        label nTgt = 0;
+        bool hasSrcMagSf = false;
+        bool hasSrcCentroids = false;
+        bool hasTgtMagSf = false;
+
+        forAll(allInterfaces, inti)
+        {
+            if (allInterfaces.set(inti))
+            {
+                const auto& intf =
+                    refCast<const cyclicACMIGAMGInterface>(allInterfaces[inti]);
+                const auto& AMI = intf.AMI();
+                nSrc += AMI.srcAddress().size();
+                nTgt += AMI.tgtAddress().size();
+
+                if (AMI.srcMagSf().size())
+                {
+                    hasSrcMagSf = true;
+                    if (AMI.srcMagSf().size() != AMI.srcAddress().size())
+                    {
+                        FatalErrorInFunction
+                            << "srcMagSf size:" << AMI.srcMagSf().size()
+                            << "srcAddress size:" << AMI.srcAddress().size()
+                            << exit(FatalError);
+                    }
+                }
+                if (AMI.srcCentroids().size())
+                {
+                    hasSrcCentroids = true;
+                    if (AMI.srcCentroids().size() != AMI.srcAddress().size())
+                    {
+                        FatalErrorInFunction
+                            << "srcCentroids size:" << AMI.srcCentroids().size()
+                            << "srcAddress size:" << AMI.srcAddress().size()
+                            << exit(FatalError);
+                    }
+                }
+                if (AMI.tgtMagSf().size())
+                {
+                    hasTgtMagSf = true;
+                    if (AMI.tgtMagSf().size() != AMI.tgtAddress().size())
+                    {
+                        FatalErrorInFunction
+                            << "tgtMagSf size:" << AMI.tgtMagSf().size()
+                            << "tgtAddress size:" << AMI.tgtAddress().size()
+                            << exit(FatalError);
+                    }
+                }
+            }
+        }
+
+
+        labelListList srcAddress;
+        scalarListList srcWeights;
+        scalarList srcMagSf;
+        // Needed?
+        pointListList srcCentroids;
+
+        labelListList tgtAddress;
+        scalarListList tgtWeights;
+        scalarList tgtMagSf;
+
+
+        // Map to send src side data to tgt side
+        autoPtr<mapDistribute> srcToTgtMap;
+
+        // Map to send tgt side data to src side
+        autoPtr<mapDistribute> tgtToSrcMap;
+
+        if (AMI.distributed())
+        {
+            // Create combined maps
+            UPtrList<const mapDistribute> srcMaps(allInterfaces.size());
+            UPtrList<const mapDistribute> tgtMaps(allInterfaces.size());
+            forAll(allInterfaces, inti)
+            {
+                if (allInterfaces.set(inti))
+                {
+                    const auto& intf = refCast<const cyclicACMIGAMGInterface>
+                    (
+                        allInterfaces[inti]
+                    );
+                    const auto& AMI = intf.AMI();
+                    srcMaps.set(inti, &AMI.srcMap());
+                    tgtMaps.set(inti, &AMI.tgtMap());
+                }
+            }
+
+
+            // Find ranks that agglomerate together
+            const label myAgglom =
+                procAgglomMap[UPstream::myProcNo(AMI.comm())];
+
+            // Invert procAgglomMap
+            const labelListList newToOldRanks
+            (
+                invertOneToMany
+                (
+                    UPstream::nProcs(coarseComm),
+                    procAgglomMap
+                )
+            );
+            const labelList& localRanks = newToOldRanks[myAgglom];
+
+
+            // Offsets for slots into results of srcToTgtMap
+            labelList srcStartOfLocal;
+            List<Map<label>> srcCompactMaps;
+
+            srcToTgtMap.reset
+            (
+                new mapDistribute
+                (
+                    srcMaps,
+                    localRanks,     // per src map which rank it is from
+                    coarseComm,
+                    newToOldRanks,  // destination rank to source ranks
+                    srcStartOfLocal,
+                    srcCompactMaps
+                )
+            );
+
+            // Assemble tgtAddress
+            tgtAddress.setSize(nTgt);
+            if (tgtAddress.size())
+            {
+                label alli = 0;
+                forAll(allInterfaces, inti)
+                {
+                    if (allInterfaces.set(inti))
+                    {
+                        const auto& intf =
+                            refCast<const cyclicACMIGAMGInterface>
+                            (
+                                allInterfaces[inti]
+                            );
+                        const auto& AMI = intf.AMI();
+                        const auto& tgtSlots = AMI.tgtAddress();
+                        const label localSize =
+                            srcStartOfLocal[inti+1]
+                          - srcStartOfLocal[inti];
+
+                        forAll(tgtSlots, tgti)
+                        {
+                            // Append old slots: copy old values and adapt
+                            auto& newSlots = tgtAddress[alli++];
+                            newSlots = tgtSlots[tgti];
+
+                            // Renumber to new indices
+                            mapDistributeBase::renumberMap
+                            (
+                                newSlots,
+                                localSize,
+                                srcStartOfLocal[inti],
+                                srcCompactMaps[inti],
+                                AMI.srcMap().constructHasFlip() //hasFlip
+                            );
+
+                            for (const label slot : newSlots)
+                            {
+                                if
+                                (
+                                    slot < 0
+                                 || slot >= srcToTgtMap().constructSize()
+                                )
+                                {
+                                    FatalErrorInFunction << " newSlots:"
+                                        << newSlots << exit(FatalError);
+                                }
+                            }
+                        }
+                    }
+                }
+
+                if (nTgt != alli)
+                {
+                    FatalErrorInFunction << "nTgt:" << nTgt
+                        << " alli:" << alli << exit(FatalError);
+                }
+            }
+
+            // Offsets for slots into results of tgtToSrcMap
+            labelList tgtStartOfLocal;
+            List<Map<label>> tgtCompactMaps;
+
+            tgtToSrcMap.reset
+            (
+                new mapDistribute
+                (
+                    tgtMaps,
+                    localRanks,
+                    coarseComm,
+                    newToOldRanks,
+                    tgtStartOfLocal,
+                    tgtCompactMaps
+                )
+            );
+
+            // Assemble srcAddress
+            srcAddress.setSize(nSrc);
+            if (srcAddress.size())
+            {
+                label alli = 0;
+                forAll(allInterfaces, inti)
+                {
+                    if (allInterfaces.set(inti))
+                    {
+                        const auto& intf =
+                            refCast<const cyclicACMIGAMGInterface>
+                            (
+                                allInterfaces[inti]
+                            );
+                        const auto& AMI = intf.AMI();
+                        const auto& srcSlots = AMI.srcAddress();
+                        const label localSize =
+                            tgtStartOfLocal[inti+1]
+                          - tgtStartOfLocal[inti];
+
+                        forAll(srcSlots, srci)
+                        {
+                            // Append old slots: copy old values and adapt
+                            auto& newSlots = srcAddress[alli++];
+                            newSlots = srcSlots[srci];
+                            // Renumber to new indices
+                            mapDistributeBase::renumberMap
+                            (
+                                newSlots,
+                                localSize,
+                                tgtStartOfLocal[inti],
+                                tgtCompactMaps[inti],
+                                AMI.tgtMap().constructHasFlip() //hasFlip
+                            );
+
+                            for (const label slot : newSlots)
+                            {
+                                if
+                                (
+                                    slot < 0
+                                 || slot >= tgtToSrcMap().constructSize()
+                                )
+                                {
+                                    FatalErrorInFunction << " newSlots:"
+                                        << newSlots << exit(FatalError);
+                                }
+                            }
+                        }
+                    }
+                }
+
+                if (nSrc != alli)
+                {
+                    FatalErrorInFunction << "nSrc:" << nSrc
+                        << " alli:" << alli << exit(FatalError);
+                }
+            }
+
+
+            // Clean up: if no remote elements sent/received mark as
+            // non-distributed. We could do this at the start but this
+            // needs to take all the internal transport into account. Easier
+            // (but less efficient) to do afterwards now all is compacted.
+            {
+                const auto& map = srcToTgtMap().subMap();
+
+                bool usesRemote = false;
+                forAll(map, proci)
+                {
+                    if (proci != myAgglom)
+                    {
+                        const auto& ss = srcToTgtMap().subMap()[proci];
+                        const auto& sc = srcToTgtMap().constructMap()[proci];
+                        const auto& ts = tgtToSrcMap().subMap()[proci];
+                        const auto& tc = tgtToSrcMap().constructMap()[proci];
+
+                        if (ss.size() || sc.size() || ts.size() || tc.size())
+                        {
+                            usesRemote = true;
+                            break;
+                        }
+                    }
+                }
+
+                if (!usesRemote)
+                {
+                    //Pout<< "** making fully local on new rank "
+                    //    << myAgglom << " in comm:" << coarseComm << endl;
+                    singlePatchProc = myAgglom;
+                    srcToTgtMap.clear();
+                    tgtToSrcMap.clear();
+                }
+            }
+        }
+        else
+        {
+            // src/tgt address are straight indices
+
+            srcAddress.setSize(nSrc);
+            tgtAddress.setSize(nTgt);
+
+            nSrc = 0;
+            nTgt = 0;
+            forAll(allInterfaces, inti)
+            {
+                if (allInterfaces.set(inti))
+                {
+                    const auto& intf = refCast<const cyclicACMIGAMGInterface>
+                    (
+                        allInterfaces[inti]
+                    );
+                    const auto& AMI = intf.AMI();
+
+                    const auto& srcA = AMI.srcAddress();
+                    if (srcAddress.size())
+                    {
+                        label srci = nSrc;
+                        forAll(srcA, i)
+                        {
+                            srcAddress[srci++] = srcA[i]+nTgt;
+                        }
+                    }
+
+                    const auto& tgtA = AMI.tgtAddress();
+                    if (tgtAddress.size())
+                    {
+                        label tgti = nTgt;
+                        forAll(tgtA, i)
+                        {
+                            tgtAddress[tgti++] = tgtA[i]+nSrc;
+                        }
+                    }
+
+                    nSrc += srcA.size();
+                    nTgt += tgtA.size();
+                }
+            }
+        }
+
+        srcWeights.setSize(nSrc);
+        if (hasSrcMagSf)
+        {
+            srcMagSf.setSize(nSrc);
+        }
+        if (hasSrcCentroids)
+        {
+            srcCentroids.setSize(nSrc);
+        }
+        tgtWeights.setSize(nTgt);
+        if (hasTgtMagSf)
+        {
+            tgtMagSf.setSize(nTgt);
+        }
+
+
+        // Append individual data
+        nSrc = 0;
+        nTgt = 0;
+        forAll(allInterfaces, inti)
+        {
+            if (allInterfaces.set(inti))
+            {
+                const auto& intf = refCast<const cyclicACMIGAMGInterface>
+                (
+                    allInterfaces[inti]
+                );
+                const auto& AMI = intf.AMI();
+
+                const auto& srcA = AMI.srcAddress();
+                {
+                    // weights
+                    SubList<scalarList>(srcWeights, srcA.size(), nSrc) =
+                        AMI.srcWeights();
+
+                    // magSf
+                    if (hasSrcMagSf)
+                    {
+                        SubList<scalar>(srcMagSf, srcA.size(), nSrc) =
+                            AMI.srcMagSf();
+                    }
+
+                    // centroids
+                    if (hasSrcCentroids)
+                    {
+                        SubList<pointList>(srcCentroids, srcA.size(), nSrc) =
+                            AMI.srcCentroids();
+                    }
+                }
+
+                const auto& tgtA = AMI.tgtAddress();
+                {
+                    // weights
+                    SubList<scalarList>(tgtWeights, tgtA.size(), nTgt) =
+                        AMI.tgtWeights();
+
+                    if (hasTgtMagSf)
+                    {
+                        SubList<scalar>(tgtMagSf, tgtA.size(), nTgt) =
+                            AMI.tgtMagSf();
+                    }
+                }
+
+                nSrc += srcA.size();
+                nTgt += tgtA.size();
+            }
+        }
+
+
+        // Construct with same arguments as original
+        amiPtr_.reset
+        (
+            new AMIPatchToPatchInterpolation
+            (
+                AMI.requireMatch(),
+                AMI.reverseTarget(),
+                AMI.lowWeightCorrection()
+            )
+        );
+        amiPtr_().comm(coarseComm),
+        amiPtr_().reset
+        (
+            std::move(srcToTgtMap),
+            std::move(tgtToSrcMap),
+            std::move(srcAddress),
+            std::move(srcWeights),
+            std::move(tgtAddress),
+            std::move(tgtWeights),
+            singlePatchProc
+        );
+        amiPtr_().srcMagSf() = std::move(srcMagSf);
+        amiPtr_().srcCentroids() = std::move(srcCentroids);
+        amiPtr_().tgtMagSf() = std::move(tgtMagSf);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -188,4 +720,24 @@ Foam::cyclicACMIGAMGInterface::internalFieldTransfer
 }
 
 
+void Foam::cyclicACMIGAMGInterface::write(Ostream& os) const
+{
+    GAMGInterface::write(os);
+
+    const bool hasAMI = amiPtr_.valid();
+
+    os  << token::SPACE << neighbPatchID_
+        << token::SPACE << owner_
+        << token::SPACE << forwardT_
+        << token::SPACE << reverseT_
+        << token::SPACE << hasAMI;
+
+    if (hasAMI)
+    {
+        os  << token::SPACE;
+        AMI().writeData(os);
+    }
+}
+
+
 // ************************************************************************* //

src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,9 +57,13 @@ class cyclicACMIGAMGInterface
 {
     // Private data
 
-        //- Reference for the cyclicLduInterface from which this is
-        //  agglomerated
-        const cyclicACMILduInterface& fineCyclicACMIInterface_;
+        label neighbPatchID_;
+
+        bool owner_;
+
+        tensorField forwardT_;
+
+        tensorField reverseT_;
 
         //- AMI interface
         autoPtr<AMIPatchToPatchInterpolation> amiPtr_;
@@ -94,9 +99,68 @@ public:
             const label coarseComm
         );
 
+        //- Construct from Istream
+        cyclicACMIGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            Istream& is
+        );
+
+        //- Construct from fine level interface,
+        //- local and neighbour restrict addressing
+        cyclicACMIGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const lduInterface& fineInterface,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        );
+
+        //- Construct by assembling and returning a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            return autoPtr<GAMGInterface>
+            (
+                new cyclicACMIGAMGInterface
+                (
+                    index,
+                    coarseInterfaces,
+                    *this,
+                    interfaceMap,
+                    faceCells,
+                    faceRestrictAddresssing,
+                    faceOffsets,
+                    allInterfaces,
+                    coarseComm,
+                    myProcNo,
+                    procAgglomMap
+                )
+            );
+        }
+
 
     //- Destructor
-    virtual ~cyclicACMIGAMGInterface();
+    virtual ~cyclicACMIGAMGInterface() = default;
 
 
     // Member Functions
@@ -116,12 +180,12 @@ public:
             //- Return neighbour processor number
             virtual label neighbPatchID() const
             {
-                return fineCyclicACMIInterface_.neighbPatchID();
+                return neighbPatchID_;
             }
 
             virtual bool owner() const
             {
-                return fineCyclicACMIInterface_.owner();
+                return owner_;
             }
 
             virtual const cyclicACMIGAMGInterface& neighbPatch() const
@@ -140,25 +204,20 @@ public:
             //- Return face transformation tensor
             virtual const tensorField& forwardT() const
             {
-                return fineCyclicACMIInterface_.forwardT();
+                return forwardT_;
             }
 
             //- Return neighbour-cell transformation tensor
             virtual const tensorField& reverseT() const
             {
-                return fineCyclicACMIInterface_.reverseT();
+                return reverseT_;
             }
 
 
         // I/O
 
             //- Write to stream
-            virtual void write(Ostream&) const
-            {
-                //TBD. How to serialise the AMI such that we can stream
-                // cyclicACMIGAMGInterface.
-                NotImplemented;
-            }
+            virtual void write(Ostream&) const;
 };
 
 

src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,12 @@ namespace Foam
         cyclicAMIGAMGInterface,
         lduInterface
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterface,
+        cyclicAMIGAMGInterface,
+        Istream
+    );
 }
 
 
@@ -63,11 +69,26 @@ Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
         index,
         coarseInterfaces
     ),
-    fineCyclicAMIInterface_
+    neighbPatchID_
     (
-        refCast<const cyclicAMILduInterface>(fineInterface)
+        refCast<const cyclicAMILduInterface>(fineInterface).neighbPatchID()
+    ),
+    owner_
+    (
+        refCast<const cyclicAMILduInterface>(fineInterface).owner()
+    ),
+    forwardT_
+    (
+        refCast<const cyclicAMILduInterface>(fineInterface).forwardT()
+    ),
+    reverseT_
+    (
+        refCast<const cyclicAMILduInterface>(fineInterface).reverseT()
     )
 {
+    const auto& fineCyclicAMIInterface =
+        refCast<const cyclicAMILduInterface>(fineInterface);
+
     // Construct face agglomeration from cell agglomeration
     {
         // From coarse face to cell
@@ -107,7 +128,7 @@ Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
 
     // On the owner side construct the AMI
 
-    if (fineCyclicAMIInterface_.owner())
+    if (fineCyclicAMIInterface.owner())
     {
         // Construct the neighbour side agglomeration (as the neighbour would
         // do it so it the exact loop above using neighbourRestrictAddressing
@@ -144,23 +165,775 @@ Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
             nbrFaceRestrictAddressing.transfer(dynNbrFaceRestrictAddressing);
         }
 
+
         amiPtr_.reset
         (
             new AMIPatchToPatchInterpolation
             (
-                fineCyclicAMIInterface_.AMI(),
+                fineCyclicAMIInterface.AMI(),
                 faceRestrictAddressing_,
                 nbrFaceRestrictAddressing
             )
         );
+
+
+        const auto& AMI = amiPtr_();
+        if (debug & 2)
+        {
+            const auto oldWarnComm = UPstream::warnComm;
+            UPstream::warnComm = AMI.comm();
+
+            const label myRank = UPstream::myProcNo(AMI.comm());
+            Pout<< "At level:" << fineLevelIndex
+                << " agglomerating from ownsize:"
+                << fineInterface.faceCells().size()
+                << " nbrSize:" << neighbourRestrictAddressing.size()
+                << " down to ownsize:" << AMI.srcAddress().size()
+                << " nbrsize:" << AMI.tgtAddress().size()
+                << " Patch:" << index << " comm:" << AMI.comm()
+                << " nProcs:" << UPstream::nProcs(AMI.comm())
+                << " myRank:" << myRank << " agglomerated AMI:"
+                << endl;
+
+            const label nbrSize = AMI.tgtAddress().size();
+            // From from nbr to owner side
+            {
+                Pout<< "From nbr:" << nbrSize << " to owner:" << this->size()
+                    << endl;
+
+                const auto& addresses = AMI.srcAddress();
+                const auto& weights = AMI.srcWeights();
+
+                labelList globalIDs(identity(nbrSize));
+                if (AMI.distributed())
+                {
+                    const auto& map = AMI.tgtMap();
+                    forAll(map.subMap(), proci)
+                    {
+                        Pout<< "    TGTMap: sending to rank:" << proci
+                            << " elements:" << flatOutput(map.subMap()[proci])
+                            << endl;
+                    }
+                    forAll(map.constructMap(), proci)
+                    {
+                        Pout<< "    TGTMap: receiving from rank:" << proci
+                            << " elements:"
+                            << flatOutput(map.constructMap()[proci])
+                            << endl;
+                    }
+                    // Fetch remote global IDs
+                    const globalIndex globalFaces(nbrSize, AMI.comm());
+                    Pout<< "    localNbrSize:" << nbrSize
+                        << " globalSize:" << globalFaces.totalSize() << endl;
+
+                    //const label myOffset = globalFaces.offsets()[myRank];
+                    for (label& id : globalIDs)
+                    {
+                        id = globalFaces.toGlobal(myRank, id);
+                    }
+                    map.distribute(globalIDs);
+                }
+
+                // Renumber my slots so they are now global face numbers
+                forAll(addresses, facei)
+                {
+                    Pout<< "    source face:" << facei
+                        << " have weights:"
+                        << flatOutput(weights[facei])
+                        << " from slots:" << flatOutput(addresses[facei])
+                        << " from global tgt faces:"
+                        << UIndirectList<label>(globalIDs, addresses[facei])
+                        << endl;
+                }
+            }
+            // From from owner to nbr side
+            {
+                Pout<< "From owner:" << this->size() << " to nbr:" << nbrSize
+                    << endl;
+
+                const auto& addresses = AMI.tgtAddress();
+                const auto& weights = AMI.tgtWeights();
+
+                labelList globalIDs(identity(this->size()));
+                if (AMI.distributed())
+                {
+                    const auto& map = AMI.srcMap();
+                    forAll(map.subMap(), proci)
+                    {
+                        Pout<< "    SRCMap: sending to rank:" << proci
+                            << " elements:" << flatOutput(map.subMap()[proci])
+                            << endl;
+                    }
+                    forAll(map.constructMap(), proci)
+                    {
+                        Pout<< "    SRCMap: receiving from rank:" << proci
+                            << " elements:"
+                            << flatOutput(map.constructMap()[proci])
+                            << endl;
+                    }
+                    // Fetch remote global IDs
+                    const globalIndex globalFaces(this->size(), AMI.comm());
+                    Pout<< "    localSize:" << this->size()
+                        << " globalSize:" << globalFaces.totalSize() << endl;
+
+                    for (label& id : globalIDs)
+                    {
+                        id = globalFaces.toGlobal(myRank, id);
+                    }
+                    map.distribute(globalIDs);
+                }
+
+                // Renumber my slots so they are now global face numbers
+                forAll(addresses, facei)
+                {
+                    Pout<< "    target face:" << facei
+                        << " have weights:"
+                        << flatOutput(weights[facei])
+                        << " from slots:" << flatOutput(addresses[facei])
+                        << " from global src faces:"
+                        << UIndirectList<label>(globalIDs, addresses[facei])
+                        << endl;
+                }
+            }
+            Pout<< "DONE agglomerating at level:" << fineLevelIndex << endl;
+
+            UPstream::warnComm = oldWarnComm;
+        }
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    Istream& is
+)
+:
+    GAMGInterface(index, coarseInterfaces, is),
+    neighbPatchID_(readLabel(is)),
+    owner_(readBool(is)),
+    forwardT_(is),
+    reverseT_(is)
+{
+    const bool hasAMI(readBool(is));
 
-Foam::cyclicAMIGAMGInterface::~cyclicAMIGAMGInterface()
-{}
+    if (hasAMI)
+    {
+        amiPtr_.reset(new AMIPatchToPatchInterpolation(is));
+    }
+}
+
+
+Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    const lduInterface& fineInterface,
+    const labelList& interfaceMap,
+    const labelUList& faceCells,
+    const labelUList& faceRestrictAddresssing,
+    const labelUList& faceOffsets,
+    const lduInterfacePtrsList& allInterfaces,
+    const label coarseComm,
+    const label myProcNo,
+    const labelList& procAgglomMap
+)
+:
+    GAMGInterface
+    (
+        index,
+        coarseInterfaces,
+        faceCells,
+        faceRestrictAddresssing
+    ),
+    neighbPatchID_
+    (
+        interfaceMap.find
+        (
+            refCast
+            <
+                const cyclicAMILduInterface
+            >(fineInterface).neighbPatchID()
+        )
+    ),
+    owner_
+    (
+        refCast<const cyclicAMILduInterface>(fineInterface).owner()
+    ),
+    forwardT_
+    (
+        refCast<const cyclicAMILduInterface>(fineInterface).forwardT()
+    ),
+    reverseT_
+    (
+       refCast<const cyclicAMILduInterface>(fineInterface).reverseT()
+    )
+{
+    const auto& fineCyclicAMIInterface =
+        refCast<const cyclicAMIGAMGInterface>(fineInterface);
+
+    if (fineCyclicAMIInterface.amiPtr_)
+    {
+        const auto& fineAMI = const_cast<AMIPatchToPatchInterpolation&>
+        (
+            fineCyclicAMIInterface.AMI()
+        );
+
+        label singlePatchProc = fineAMI.singlePatchProc();
+
+
+
+        // Get some sizes
+        label nSrc = 0;
+        label nTgt = 0;
+        bool hasSrcMagSf = false;
+        bool hasSrcCentroids = false;
+        bool hasTgtMagSf = false;
+
+        forAll(allInterfaces, inti)
+        {
+            if (allInterfaces.set(inti))
+            {
+                const auto& intf =
+                    refCast<const cyclicAMIGAMGInterface>(allInterfaces[inti]);
+                const auto& AMI = intf.AMI();
+                nSrc += AMI.srcAddress().size();
+                nTgt += AMI.tgtAddress().size();
+
+                if (AMI.srcMagSf().size())
+                {
+                    hasSrcMagSf = true;
+                    if (AMI.srcMagSf().size() != AMI.srcAddress().size())
+                    {
+                        FatalErrorInFunction
+                            << "srcMagSf size:" << AMI.srcMagSf().size()
+                            << "srcAddress size:" << AMI.srcAddress().size()
+                            << exit(FatalError);
+                    }
+                }
+                if (AMI.srcCentroids().size())
+                {
+                    hasSrcCentroids = true;
+                    if (AMI.srcCentroids().size() != AMI.srcAddress().size())
+                    {
+                        FatalErrorInFunction
+                            << "srcCentroids size:" << AMI.srcCentroids().size()
+                            << "srcAddress size:" << AMI.srcAddress().size()
+                            << exit(FatalError);
+                    }
+                }
+                if (AMI.tgtMagSf().size())
+                {
+                    hasTgtMagSf = true;
+                    if (AMI.tgtMagSf().size() != AMI.tgtAddress().size())
+                    {
+                        FatalErrorInFunction
+                            << "tgtMagSf size:" << AMI.tgtMagSf().size()
+                            << "tgtAddress size:" << AMI.tgtAddress().size()
+                            << exit(FatalError);
+                    }
+                }
+            }
+        }
+
+
+        labelListList srcAddress;
+        scalarListList srcWeights;
+        scalarList srcMagSf;
+        // Needed?
+        pointListList srcCentroids;
+
+        labelListList tgtAddress;
+        scalarListList tgtWeights;
+        scalarList tgtMagSf;
+
+
+        // Map to send src side data to tgt side
+        autoPtr<mapDistribute> srcToTgtMap;
+
+        // Map to send tgt side data to src side
+        autoPtr<mapDistribute> tgtToSrcMap;
+
+        if (fineAMI.distributed())
+        {
+            // Create combined maps
+            UPtrList<const mapDistribute> srcMaps(allInterfaces.size());
+            UPtrList<const mapDistribute> tgtMaps(allInterfaces.size());
+            forAll(allInterfaces, inti)
+            {
+                if (allInterfaces.set(inti))
+                {
+                    const auto& intf = refCast<const cyclicAMIGAMGInterface>
+                    (
+                        allInterfaces[inti]
+                    );
+                    const auto& AMI = intf.AMI();
+                    srcMaps.set(inti, &AMI.srcMap());
+                    tgtMaps.set(inti, &AMI.tgtMap());
+                }
+            }
+
+
+            // Find ranks that agglomerate together
+            const label myAgglom =
+                procAgglomMap[UPstream::myProcNo(fineAMI.comm())];
+
+            // Invert procAgglomMap
+            const labelListList newToOldRanks
+            (
+                invertOneToMany
+                (
+                    UPstream::nProcs(coarseComm),
+                    procAgglomMap
+                )
+            );
+            const labelList& localRanks = newToOldRanks[myAgglom];
+
+
+            // Offsets for slots into results of srcToTgtMap
+            labelList srcStartOfLocal;
+            List<Map<label>> srcCompactMaps;
+
+            srcToTgtMap.reset
+            (
+                new mapDistribute
+                (
+                    srcMaps,
+                    localRanks,     // per src map which rank it is from
+                    coarseComm,
+                    newToOldRanks,  // destination rank to source ranks
+                    srcStartOfLocal,
+                    srcCompactMaps
+                )
+            );
+
+            // Assemble tgtAddress
+            tgtAddress.setSize(nTgt);
+            if (tgtAddress.size())
+            {
+                label alli = 0;
+                forAll(allInterfaces, inti)
+                {
+                    if (allInterfaces.set(inti))
+                    {
+                        const auto& intf = refCast<const cyclicAMIGAMGInterface>
+                        (
+                            allInterfaces[inti]
+                        );
+                        const auto& AMI = intf.AMI();
+                        const auto& tgtSlots = AMI.tgtAddress();
+                        const label localSize =
+                            srcStartOfLocal[inti+1]
+                          - srcStartOfLocal[inti];
+
+                        forAll(tgtSlots, tgti)
+                        {
+                            // Append old slots: copy old values and adapt
+                            auto& newSlots = tgtAddress[alli++];
+                            newSlots = tgtSlots[tgti];
+
+                            // Renumber to new indices
+                            mapDistributeBase::renumberMap
+                            (
+                                newSlots,
+                                localSize,
+                                srcStartOfLocal[inti],
+                                srcCompactMaps[inti],
+                                AMI.srcMap().constructHasFlip() //hasFlip
+                            );
+
+                            for (const label slot : newSlots)
+                            {
+                                if
+                                (
+                                    slot < 0
+                                 || slot >= srcToTgtMap().constructSize()
+                                )
+                                {
+                                    FatalErrorInFunction << " newSlots:"
+                                        << newSlots << exit(FatalError);
+                                }
+                            }
+                        }
+                    }
+                }
+
+                if (nTgt != alli)
+                {
+                    FatalErrorInFunction << "nTgt:" << nTgt
+                        << " alli:" << alli << exit(FatalError);
+                }
+            }
+
+            // Offsets for slots into results of tgtToSrcMap
+            labelList tgtStartOfLocal;
+            List<Map<label>> tgtCompactMaps;
+
+            tgtToSrcMap.reset
+            (
+                new mapDistribute
+                (
+                    tgtMaps,
+                    localRanks,
+                    coarseComm,
+                    newToOldRanks,
+                    tgtStartOfLocal,
+                    tgtCompactMaps
+                )
+            );
+
+            // Assemble srcAddress
+            srcAddress.setSize(nSrc);
+            if (srcAddress.size())
+            {
+                label alli = 0;
+                forAll(allInterfaces, inti)
+                {
+                    if (allInterfaces.set(inti))
+                    {
+                        const auto& intf = refCast<const cyclicAMIGAMGInterface>
+                        (
+                            allInterfaces[inti]
+                        );
+                        const auto& AMI = intf.AMI();
+                        const auto& srcSlots = AMI.srcAddress();
+                        const label localSize =
+                            tgtStartOfLocal[inti+1]
+                          - tgtStartOfLocal[inti];
+
+                        forAll(srcSlots, srci)
+                        {
+                            // Append old slots: copy old values and adapt
+                            auto& newSlots = srcAddress[alli++];
+                            newSlots = srcSlots[srci];
+                            // Renumber to new indices
+                            mapDistributeBase::renumberMap
+                            (
+                                newSlots,
+                                localSize,
+                                tgtStartOfLocal[inti],
+                                tgtCompactMaps[inti],
+                                AMI.tgtMap().constructHasFlip() //hasFlip
+                            );
+
+                            for (const label slot : newSlots)
+                            {
+                                if
+                                (
+                                    slot < 0
+                                 || slot >= tgtToSrcMap().constructSize()
+                                )
+                                {
+                                    FatalErrorInFunction << " newSlots:"
+                                        << newSlots << exit(FatalError);
+                                }
+                            }
+                        }
+                    }
+                }
+
+                if (nSrc != alli)
+                {
+                    FatalErrorInFunction << "nSrc:" << nSrc
+                        << " alli:" << alli << exit(FatalError);
+                }
+            }
+
+
+            // Clean up: if no remote elements sent/received mark as
+            // non-distributed. We could do this at the start but this
+            // needs to take all the internal transport into account. Easier
+            // (but less efficient) to do afterwards now all is compacted.
+            {
+                const auto& map = srcToTgtMap().subMap();
+
+                bool usesRemote = false;
+                forAll(map, proci)
+                {
+                    if (proci != myAgglom)
+                    {
+                        const auto& ss = srcToTgtMap().subMap()[proci];
+                        const auto& sc = srcToTgtMap().constructMap()[proci];
+                        const auto& ts = tgtToSrcMap().subMap()[proci];
+                        const auto& tc = tgtToSrcMap().constructMap()[proci];
+
+                        if (ss.size() || sc.size() || ts.size() || tc.size())
+                        {
+                            usesRemote = true;
+                            break;
+                        }
+                    }
+                }
+
+                if (!usesRemote)
+                {
+                    //Pout<< "** making fully local on new rank "
+                    //    << myAgglom << " in comm:" << coarseComm << endl;
+                    singlePatchProc = myAgglom;
+                    srcToTgtMap.clear();
+                    tgtToSrcMap.clear();
+                }
+            }
+        }
+        else
+        {
+            // src/tgt address are straight indices
+
+            srcAddress.setSize(nSrc);
+            tgtAddress.setSize(nTgt);
+
+            nSrc = 0;
+            nTgt = 0;
+            forAll(allInterfaces, inti)
+            {
+                if (allInterfaces.set(inti))
+                {
+                    const auto& intf = refCast<const cyclicAMIGAMGInterface>
+                    (
+                        allInterfaces[inti]
+                    );
+                    const auto& AMI = intf.AMI();
+
+                    const auto& srcA = AMI.srcAddress();
+                    if (srcAddress.size())
+                    {
+                        label srci = nSrc;
+                        forAll(srcA, i)
+                        {
+                            srcAddress[srci++] = srcA[i]+nTgt;
+                        }
+                    }
+
+                    const auto& tgtA = AMI.tgtAddress();
+                    if (tgtAddress.size())
+                    {
+                        label tgti = nTgt;
+                        forAll(tgtA, i)
+                        {
+                            tgtAddress[tgti++] = tgtA[i]+nSrc;
+                        }
+                    }
+
+                    nSrc += srcA.size();
+                    nTgt += tgtA.size();
+                }
+            }
+        }
+
+        srcWeights.setSize(nSrc);
+        if (hasSrcMagSf)
+        {
+            srcMagSf.setSize(nSrc);
+        }
+        if (hasSrcCentroids)
+        {
+            srcCentroids.setSize(nSrc);
+        }
+        tgtWeights.setSize(nTgt);
+        if (hasTgtMagSf)
+        {
+            tgtMagSf.setSize(nTgt);
+        }
+
+
+        // Append individual data
+        nSrc = 0;
+        nTgt = 0;
+        forAll(allInterfaces, inti)
+        {
+            if (allInterfaces.set(inti))
+            {
+                const auto& intf = refCast<const cyclicAMIGAMGInterface>
+                (
+                    allInterfaces[inti]
+                );
+                const auto& AMI = intf.AMI();
+
+                const auto& srcA = AMI.srcAddress();
+                {
+                    // weights
+                    SubList<scalarList>(srcWeights, srcA.size(), nSrc) =
+                        AMI.srcWeights();
+
+                    // magSf
+                    if (hasSrcMagSf)
+                    {
+                        SubList<scalar>(srcMagSf, srcA.size(), nSrc) =
+                            AMI.srcMagSf();
+                    }
+
+                    // centroids
+                    if (hasSrcCentroids)
+                    {
+                        SubList<pointList>(srcCentroids, srcA.size(), nSrc) =
+                            AMI.srcCentroids();
+                    }
+                }
+
+                const auto& tgtA = AMI.tgtAddress();
+                {
+                    // weights
+                    SubList<scalarList>(tgtWeights, tgtA.size(), nTgt) =
+                        AMI.tgtWeights();
+
+                    if (hasTgtMagSf)
+                    {
+                        SubList<scalar>(tgtMagSf, tgtA.size(), nTgt) =
+                            AMI.tgtMagSf();
+                    }
+                }
+
+                nSrc += srcA.size();
+                nTgt += tgtA.size();
+            }
+        }
+
+
+        // Construct with same arguments as original
+        amiPtr_.reset
+        (
+            new AMIPatchToPatchInterpolation
+            (
+                fineAMI.requireMatch(),
+                fineAMI.reverseTarget(),
+                fineAMI.lowWeightCorrection()
+            )
+        );
+        amiPtr_().comm(coarseComm),
+        amiPtr_().reset
+        (
+            std::move(srcToTgtMap),
+            std::move(tgtToSrcMap),
+            std::move(srcAddress),
+            std::move(srcWeights),
+            std::move(tgtAddress),
+            std::move(tgtWeights),
+            singlePatchProc
+        );
+        amiPtr_().srcMagSf() = std::move(srcMagSf);
+        amiPtr_().srcCentroids() = std::move(srcCentroids);
+        amiPtr_().tgtMagSf() = std::move(tgtMagSf);
+
+
+        if (debug & 2)
+        {
+            const auto& AMI = amiPtr_();
+
+            const auto oldWarnComm = UPstream::warnComm;
+            UPstream::warnComm = AMI.comm();
+
+            const label myRank = UPstream::myProcNo(AMI.comm());
+            Pout<< "PROCAGGLOMERATED :"
+                << " Patch:" << index << " comm:" << AMI.comm()
+                << " nProcs:" << UPstream::nProcs(AMI.comm())
+                << " myRank:" << myRank << " agglomerated AMI:"
+                << endl;
+
+            const label nbrSize = AMI.tgtAddress().size();
+            // From from nbr to owner side
+            {
+                Pout<< "From nbr:" << nbrSize << " to owner:" << this->size()
+                    << endl;
+
+                const auto& addresses = AMI.srcAddress();
+                const auto& weights = AMI.srcWeights();
+
+                labelList globalIDs(identity(nbrSize));
+                if (AMI.distributed())
+                {
+                    const auto& map = AMI.tgtMap();
+                    forAll(map.subMap(), proci)
+                    {
+                        Pout<< "    TGTMap: sending to rank:" << proci
+                            << " elements:" << flatOutput(map.subMap()[proci])
+                            << endl;
+                    }
+                    forAll(map.constructMap(), proci)
+                    {
+                        Pout<< "    TGTMap: receiving from rank:" << proci
+                            << " elements:"
+                            << flatOutput(map.constructMap()[proci])
+                            << endl;
+                    }
+
+                    // Fetch remote global IDs
+                    const globalIndex globalFaces(nbrSize, AMI.comm());
+                    Pout<< "    localNbrSize:" << nbrSize
+                        << " globalSize:" << globalFaces.totalSize() << endl;
+                    for (label& id : globalIDs)
+                    {
+                        id = globalFaces.toGlobal(myRank, id);
+                    }
+                    map.distribute(globalIDs);
+                }
+
+                // Renumber my slots so they are now global face numbers
+                forAll(addresses, facei)
+                {
+                    Pout<< "    source face:" << facei
+                        << " have weights:"
+                        << flatOutput(weights[facei])
+                        << " from slots:" << flatOutput(addresses[facei])
+                        << " from global tgt faces:"
+                        << UIndirectList<label>(globalIDs, addresses[facei])
+                        << endl;
+                }
+            }
+            // From from owner to nbr side
+            {
+                Pout<< "From owner:" << this->size() << " to nbr:" << nbrSize
+                    << endl;
+
+                const auto& addresses = AMI.tgtAddress();
+                const auto& weights = AMI.tgtWeights();
+
+                labelList globalIDs(identity(this->size()));
+                if (AMI.distributed())
+                {
+                    const auto& map = AMI.srcMap();
+                    forAll(map.subMap(), proci)
+                    {
+                        Pout<< "    SRCMap: sending to rank:" << proci
+                            << " elements:" << flatOutput(map.subMap()[proci])
+                            << endl;
+                    }
+                    forAll(map.constructMap(), proci)
+                    {
+                        Pout<< "    SRCMap: receiving from rank:" << proci
+                            << " elements:"
+                            << flatOutput(map.constructMap()[proci])
+                            << endl;
+                    }
+
+                    // Fetch remote global IDs
+                    const globalIndex globalFaces(this->size(), AMI.comm());
+                    Pout<< "    localSize:" << this->size()
+                        << " globalSize:" << globalFaces.totalSize() << endl;
+                    for (label& id : globalIDs)
+                    {
+                        id = globalFaces.toGlobal(myRank, id);
+                    }
+                    map.distribute(globalIDs);
+                }
+
+                // Renumber my slots so they are now global face numbers
+                forAll(addresses, facei)
+                {
+                    Pout<< "    target face:" << facei
+                        << " have weights:"
+                        << flatOutput(weights[facei])
+                        << " from slots:" << flatOutput(addresses[facei])
+                        << " from global src faces:"
+                        << UIndirectList<label>(globalIDs, addresses[facei])
+                        << endl;
+                }
+            }
+            Pout<< "DONE PROCAGGLOMERATED" << endl;
+            UPstream::warnComm = oldWarnComm;
+        }
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -175,7 +948,7 @@ Foam::tmp<Foam::labelField> Foam::cyclicAMIGAMGInterface::internalFieldTransfer
         dynamic_cast<const cyclicAMIGAMGInterface&>(neighbPatch());
     const labelUList& nbrFaceCells = nbr.faceCells();
 
-    tmp<labelField> tpnf(new labelField(nbrFaceCells.size()));
+    auto tpnf = tmp<labelField>::New(nbrFaceCells.size());
     labelField& pnf = tpnf.ref();
 
     forAll(pnf, facei)
@@ -186,4 +959,25 @@ Foam::tmp<Foam::labelField> Foam::cyclicAMIGAMGInterface::internalFieldTransfer
     return tpnf;
 }
 
+
+void Foam::cyclicAMIGAMGInterface::write(Ostream& os) const
+{
+    GAMGInterface::write(os);
+
+    const bool hasAMI = amiPtr_.valid();
+
+    os  << token::SPACE << neighbPatchID_
+        << token::SPACE << owner_
+        << token::SPACE << forwardT_
+        << token::SPACE << reverseT_
+        << token::SPACE << hasAMI;
+
+    if (hasAMI)
+    {
+        os  << token::SPACE;
+        AMI().writeData(os);
+    }
+}
+
+
 // ************************************************************************* //

src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,9 +57,13 @@ class cyclicAMIGAMGInterface
 {
     // Private data
 
-        //- Reference for the cyclicLduInterface from which this is
-        //  agglomerated
-        const cyclicAMILduInterface& fineCyclicAMIInterface_;
+        label neighbPatchID_;
+
+        bool owner_;
+
+        tensorField forwardT_;
+
+        tensorField reverseT_;
 
         //- AMI interface
         autoPtr<AMIPatchToPatchInterpolation> amiPtr_;
@@ -94,9 +99,68 @@ public:
             const label coarseComm
         );
 
+        //- Construct from Istream
+        cyclicAMIGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            Istream& is
+        );
+
+        //- Construct from fine level interface,
+        //- local and neighbour restrict addressing
+        cyclicAMIGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const lduInterface& fineInterface,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        );
+
+        //- Construct by assembling and returning a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            return autoPtr<GAMGInterface>
+            (
+                new cyclicAMIGAMGInterface
+                (
+                    index,
+                    coarseInterfaces,
+                    *this,
+                    interfaceMap,
+                    faceCells,
+                    faceRestrictAddresssing,
+                    faceOffsets,
+                    allInterfaces,
+                    coarseComm,
+                    myProcNo,
+                    procAgglomMap
+                )
+            );
+        }
+
 
     //- Destructor
-    virtual ~cyclicAMIGAMGInterface();
+    virtual ~cyclicAMIGAMGInterface() = default;
 
 
     // Member Functions
@@ -116,12 +180,12 @@ public:
             //- Return neighbour processor number
             virtual label neighbPatchID() const
             {
-                return fineCyclicAMIInterface_.neighbPatchID();
+                return neighbPatchID_;
             }
 
             virtual bool owner() const
             {
-                return fineCyclicAMIInterface_.owner();
+                return owner_;
             }
 
             virtual const cyclicAMIGAMGInterface& neighbPatch() const
@@ -140,25 +204,20 @@ public:
             //- Return face transformation tensor
             virtual const tensorField& forwardT() const
             {
-                return fineCyclicAMIInterface_.forwardT();
+                return forwardT_;
             }
 
             //- Return neighbour-cell transformation tensor
             virtual const tensorField& reverseT() const
             {
-                return fineCyclicAMIInterface_.reverseT();
+                return reverseT_;
             }
 
 
         // I/O
 
             //- Write to stream
-            virtual void write(Ostream&) const
-            {
-                //TBD. How to serialise the AMI such that we can stream
-                // cyclicAMIGAMGInterface.
-                NotImplemented;
-            }
+            virtual void write(Ostream&) const;
 };
 
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C

@@ -302,6 +302,7 @@ void Foam::cyclicAMIPolyPatch::setAMIFaces()
     const labelListList& tgtToSrcAddr0 = AMIPtr_->tgtAddress();
     const pointListList& srcCtr0 = AMIPtr_->srcCentroids();
     const scalarListList& srcToTgtWght0 = AMIPtr_->srcWeights();
+    const label singlePatchProc = AMIPtr_->singlePatchProc();
 
     // New addressing on new mesh (extended by polyTopoChange)
     labelListList srcToTgtAddr1(size(), labelList());
@@ -316,12 +317,12 @@ void Foam::cyclicAMIPolyPatch::setAMIFaces()
        // Parallel running
 
         // Global index based on old patch sizes (when AMI was computed)
-        globalIndex globalSrcFaces0(srcToTgtAddr0.size());
-        globalIndex globalTgtFaces0(tgtToSrcAddr0.size());
+        globalIndex globalSrcFaces0(srcToTgtAddr0.size(), AMIPtr_().comm());
+        globalIndex globalTgtFaces0(tgtToSrcAddr0.size(), AMIPtr_().comm());
 
         // Global index based on new patch sizes
-        globalIndex globalSrcFaces1(size());
-        globalIndex globalTgtFaces1(nbr.size());
+        globalIndex globalSrcFaces1(size(), AMIPtr_().comm());
+        globalIndex globalTgtFaces1(nbr.size(), AMIPtr_().comm());
 
 
         // Gather source side info
@@ -623,7 +624,8 @@ void Foam::cyclicAMIPolyPatch::setAMIFaces()
         std::move(srcToTgtAddr1),
         std::move(newSrcToTgtWeights),
         std::move(tgtToSrcAddr1),
-        std::move(newTgtToSrcWeights)
+        std::move(newTgtToSrcWeights),
+        singlePatchProc
     );
 
     // Need to set areas, e.g. for agglomeration to (re-)normalisation weights

src/overset/lduPrimitiveProcessorInterface/GAMG/calculatedProcessorGAMGInterfaceField.C

@@ -46,6 +46,12 @@ namespace Foam
         calculatedProcessorGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        calculatedProcessorGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -89,6 +95,20 @@ calculatedProcessorGAMGInterfaceField
 {}
 
 
+Foam::calculatedProcessorGAMGInterfaceField::
+calculatedProcessorGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    procInterface_(refCast<const calculatedProcessorGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::calculatedProcessorGAMGInterfaceField::initInterfaceMatrixUpdate
@@ -194,4 +214,12 @@ void Foam::calculatedProcessorGAMGInterfaceField::updateInterfaceMatrix
 }
 
 
+void Foam::calculatedProcessorGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
+}
+
+
 // ************************************************************************* //

src/overset/lduPrimitiveProcessorInterface/GAMG/calculatedProcessorGAMGInterfaceField.H

@@ -118,6 +118,24 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        calculatedProcessorGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            NotImplemented;
+            return autoPtr<GAMGInterfaceField>(nullptr);
+        }
+
 
     //- Destructor
     virtual ~calculatedProcessorGAMGInterfaceField() = default;
@@ -200,6 +218,12 @@ public:
             {
                 return rank_;
             }
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/overset/oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,12 @@ namespace Foam
         oversetGAMGInterface,
         lduInterface
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterface,
+        oversetGAMGInterface,
+        Istream
+    );
 }
 
 
@@ -110,9 +116,36 @@ Foam::oversetGAMGInterface::oversetGAMGInterface
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+Foam::oversetGAMGInterface::oversetGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    Istream& is
+)
+:
+    GAMGInterface(index, coarseInterfaces, is)
+{}
 
-Foam::oversetGAMGInterface::~oversetGAMGInterface()
+
+Foam::oversetGAMGInterface::oversetGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    const lduInterface& fineInterface,
+    const labelList& interfaceMap,
+    const labelUList& faceCells,
+    const labelUList& faceRestrictAddresssing,
+    const labelUList& faceOffsets,
+    const lduInterfacePtrsList& allInterfaces
+)
+:
+    GAMGInterface
+    (
+        index,
+        coarseInterfaces,
+        faceCells,
+        faceRestrictAddresssing
+    )
 {}
 
 
@@ -128,4 +161,10 @@ Foam::tmp<Foam::labelField> Foam::oversetGAMGInterface::internalFieldTransfer
 }
 
 
+void Foam::oversetGAMGInterface::write(Ostream& os) const
+{
+    GAMGInterface::write(os);
+}
+
+
 // ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenFOAM Foundation
+    Copyright (C) 2021,2023 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,9 +84,62 @@ public:
             const label coarseComm
         );
 
+        //- Construct from Istream
+        oversetGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            Istream& is
+        );
+
+        //- Construct from fine level interface,
+        //- local and neighbour restrict addressing
+        oversetGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const lduInterface& fineInterface,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces
+        );
+
+        //- Construct by assembling and returning a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            return autoPtr<GAMGInterface>
+            (
+                new oversetGAMGInterface
+                (
+                    index,
+                    coarseInterfaces,
+                    *this,
+                    interfaceMap,
+                    faceCells,
+                    faceRestrictAddresssing,
+                    faceOffsets,
+                    allInterfaces
+                )
+            );
+        }
+
 
     //- Destructor
-    virtual ~oversetGAMGInterface();
+    virtual ~oversetGAMGInterface() = default;
 
 
     // Member Functions
@@ -104,10 +157,7 @@ public:
         // I/O
 
             //- Write to stream
-            virtual void write(Ostream&) const
-            {
-                NotImplemented;
-            }
+            virtual void write(Ostream&) const;
 };
 
 

src/overset/oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.C

@@ -46,6 +46,12 @@ namespace Foam
         oversetGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        oversetGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -72,6 +78,16 @@ Foam::oversetGAMGInterfaceField::oversetGAMGInterfaceField
 {}
 
 
+Foam::oversetGAMGInterfaceField::oversetGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::oversetGAMGInterfaceField::~oversetGAMGInterfaceField()
@@ -94,4 +110,8 @@ void Foam::oversetGAMGInterfaceField::updateInterfaceMatrix
 {}
 
 
+void Foam::oversetGAMGInterfaceField::write(Ostream& os) const
+{}
+
+
 // ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.H

@@ -90,6 +90,24 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        oversetGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            NotImplemented;
+            return autoPtr<GAMGInterfaceField>(nullptr);
+        }
+
 
     //- Destructor
     virtual ~oversetGAMGInterfaceField();
@@ -111,6 +129,12 @@ public:
                 const direction cmpt,
                 const Pstream::commsTypes commsType
             ) const;
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };