ENH: AMIInterpolation: make matching more robust

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2412
-patch=250528
+patch=0

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -53,6 +53,7 @@ Description
 #include "pointSet.H"
 #include "ReadFields.H"
 #include "processorMeshes.H"
+#include "IOmapDistributePolyMesh.H"
 
 using namespace Foam;
 
@@ -369,6 +370,11 @@ int main(int argc, char *argv[])
         "Subset with cellZone(s) instead of cellSet."
         " The command argument may be a list of words or regexs"
     );
+    argList::addBoolOption
+    (
+        "no-map",
+        "Suppress writing of map."
+    );
     argList::addOption
     (
         "resultTime",
@@ -393,6 +399,7 @@ int main(int argc, char *argv[])
 
     const bool useCellZone = args.found("zone");
     const bool overwrite = args.found("overwrite");
+    const bool noMap = args.found("no-map");
     const bool specifiedInstance = args.readIfPresent
     (
         "resultTime",
@@ -663,8 +670,6 @@ int main(int argc, char *argv[])
     Info<< "Writing subsetted mesh and fields to time " << runTime.timeName()
         << endl;
     subsetter.subMesh().write();
-    processorMeshes::removeFiles(subsetter.subMesh());
-
     auto* subPointMeshPtr =
         subsetter.subMesh().thisDb().findObject<pointMesh>
         (
@@ -676,6 +681,115 @@ int main(int argc, char *argv[])
         subPointMesh.setInstance(subsetter.subMesh().facesInstance());
         subPointMesh.write();
     }
+    processorMeshes::removeFiles(subsetter.subMesh());
+
+
+    if (!noMap)
+    {
+        const auto& subMesh = subsetter.subMesh();
+        const auto& pbm = mesh.boundaryMesh();
+        labelList patchStarts(pbm.size());
+        labelList patchNMeshPoints(pbm.size());
+        for (const auto& pp : pbm)
+        {
+            patchStarts[pp.index()] = pp.start();
+            patchNMeshPoints[pp.index()] = pp.nPoints();
+        }
+
+        const label myProcNo = UPstream::myProcNo(mesh.comm());
+
+        // cellMap
+        labelListList cellSubMap(UPstream::nProcs(mesh.comm()));
+        cellSubMap[myProcNo] = subsetter.cellMap();
+        labelListList cellConstructMap(UPstream::nProcs(mesh.comm()));
+        cellConstructMap[myProcNo] = identity(subMesh.nCells());
+        mapDistribute cellMap
+        (
+            subMesh.nCells(),
+            std::move(cellSubMap),
+            std::move(cellConstructMap),
+            false,
+            false,
+            mesh.comm()
+        );
+
+        // faceMap
+        labelListList faceSubMap(UPstream::nProcs(mesh.comm()));
+        faceSubMap[myProcNo] = subsetter.faceMap();
+        labelListList faceConstructMap(UPstream::nProcs(mesh.comm()));
+        faceConstructMap[myProcNo] = identity(subMesh.nFaces());
+        mapDistribute faceMap
+        (
+            subMesh.nFaces(),
+            std::move(faceSubMap),
+            std::move(faceConstructMap),
+            false,
+            false,
+            mesh.comm()
+        );
+
+        // pointMap
+        labelListList pointSubMap(UPstream::nProcs(mesh.comm()));
+        pointSubMap[myProcNo] = subsetter.pointMap();
+        labelListList pointConstructMap(UPstream::nProcs(mesh.comm()));
+        pointConstructMap[myProcNo] = identity(subMesh.nPoints());
+        mapDistribute pointMap
+        (
+            subMesh.nPoints(),
+            std::move(pointSubMap),
+            std::move(pointConstructMap),
+            false,
+            false,
+            mesh.comm()
+        );
+
+        // patchMap
+        labelListList patchSubMap(UPstream::nProcs(mesh.comm()));
+        patchSubMap[myProcNo] = identity(pbm.size());
+        labelListList patchConstructMap(UPstream::nProcs(mesh.comm()));
+        patchConstructMap[myProcNo] = identity(pbm.size()); // or subMesh?
+        mapDistribute patchMap
+        (
+            subMesh.nPoints(),
+            std::move(patchSubMap),
+            std::move(patchConstructMap),
+            false,
+            false,
+            mesh.comm()
+        );
+
+        mapDistributePolyMesh map
+        (
+            mesh.nPoints(), // old points
+            mesh.nFaces(),  // old faces
+            mesh.nCells(),  // old cells
+
+            std::move(patchStarts),
+            std::move(patchNMeshPoints),
+
+            std::move(pointMap),
+            std::move(faceMap),
+            std::move(cellMap),
+            std::move(patchMap)
+        );
+
+
+        const IOobject io
+        (
+            "parentMeshAddressing",
+            subMesh.facesInstance(),
+            fvMesh::meshSubDir,
+            subMesh.thisDb(),
+            IOobjectOption::NO_READ,
+            IOobjectOption::NO_WRITE,
+            IOobjectOption::NO_REGISTER
+        );
+
+        Info<< "Writing map from subsetted to original mesh to "
+            << io.objectRelPath() << endl;
+
+        IOmapDistributePolyMeshRef(io, map).write();
+    }
 
 
     // Volume fields

src/finiteVolume/fvMesh/fvMesh.C

@@ -51,6 +51,26 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::refPtr<Foam::fvMesh>
+Foam::fvMesh::parentMesh(const objectRegistry& obr) const
+{
+    const fvMesh* meshPtr = isA<fvMesh>(obr);
+
+    if (meshPtr)
+    {
+        return refPtr<fvMesh>(*meshPtr);
+    }
+    else if (obr.isTimeDb())
+    {
+        return refPtr<fvMesh>();
+    }
+    else
+    {
+        return parentMesh(obr.parent());
+    }
+}
+
+
 void Foam::fvMesh::clearGeomNotOldVol()
 {
     meshObject::clearUpto

src/finiteVolume/fvMesh/fvMesh.H

@@ -164,6 +164,9 @@ protected:
             void makeCf() const;
 
 
+        //- Helper: search for parent mesh
+        refPtr<fvMesh> parentMesh(const objectRegistry& obr) const;
+
         //- No copy construct
         fvMesh(const fvMesh&) = delete;
 

src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C

@@ -183,8 +183,7 @@ void Foam::zoneDistribute::setUpCommforZone
 Foam::List<Foam::label> Foam::zoneDistribute::getCyclicPatches
 (
     const label celli,
-    const label globalIdx,
-    const vector globalIdxCellCentre
+    const label globalIdx
 ) const
 {
     // Initialise cyclic patch label list
@@ -198,26 +197,26 @@ Foam::List<Foam::label> Foam::zoneDistribute::getCyclicPatches
 
     const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
 
-    // Making list of cyclic patches to which celli belongs
-    List<label> celliCyclicPatches;
-    forAll(bMesh, patchi)
-    {
-        if (isA<cyclicPolyPatch>(bMesh[patchi]))
-        {
-            // Note: Probably not efficient due to use of found(celli) but
-            // typically only used for very few cells (interface cells and their
-            // point neighbours on cyclic boundaries).
-            if (bMesh[patchi].faceCells().found(celli))
-            {
-                celliCyclicPatches.append(patchi);
-            }
-        }
-    }
-
     // So celli belongs to at least one cyclic patch.
     // Let us figure out which.
     if (globalNumbering_.isLocal(globalIdx)) // celli and globalIdx on same proc
     {
+        // Making list of cyclic patches to which celli belongs
+        List<label> celliCyclicPatches;
+        forAll(bMesh, patchi)
+        {
+            if (isA<cyclicPolyPatch>(bMesh[patchi]))
+            {
+                // Note: Probably not efficient due to use of found(celli) but
+                // typically only used for very cells (interface cells and their
+                // point neighbours on cyclic boundaries).
+                if (bMesh[patchi].faceCells().found(celli))
+                {
+                    celliCyclicPatches.append(patchi);
+                }
+            }
+        }
+
         // Get all local point neighbor cells of celli, i.e. all point
         // neighbours that are not on the other side of a cyclic patch.
         List<label> localPointNeiCells(0);
@@ -244,14 +243,14 @@ Foam::List<Foam::label> Foam::zoneDistribute::getCyclicPatches
         {
             for (const label patchi : celliCyclicPatches)
             {
-                // Find the corresponding cyclic neighbor patch ID
+                // find the corresponding cyclic neighbor patch ID
                 const cyclicPolyPatch& cpp =
                     static_cast<const cyclicPolyPatch&>(bMesh[patchi]);
 
                 const label neiPatch = cpp.neighbPatchID();
 
                 // Check if the cell globalIdx is on neiPatch.
-                // If it is, append neiPatch to list of patches to return
+                // If it is, append neiPatch to list of
                 if (bMesh[neiPatch].faceCells().found(localIdx))
                 {
                     patches.append(neiPatch);
@@ -268,71 +267,75 @@ Foam::List<Foam::label> Foam::zoneDistribute::getCyclicPatches
     }
     else // celli and globalIdx on differet processors
     {
-        List<label> cyclicID(3, -1);
-        List<vector> separationVectors(3, vector(0,0,0));
-        scalar distance = GREAT;
+        // Note: The following is needed if a celli is located at the interface
+        // (plicRDF), on a cyclic patch and a processor patch. In this case
+        // globalIdx may be on a different processor, but requires the
+        // transformation from a cyclic patch on the processor of celli.
 
-        forAll(celliCyclicPatches, cID)
+        const List<label>& faces = mesh_.cells()[celli];
+
+        // Loop over all faces of celli and find cyclic patches
+        for (const label facei : faces)
         {
-            cyclicID[cID] = celliCyclicPatches[cID];
+            if (mesh_.isInternalFace(facei)) continue;
 
-            const label& patchI = celliCyclicPatches[cID];
-            const cyclicPolyPatch& cpp =
-                static_cast<const cyclicPolyPatch&>(bMesh[patchI]);
+            const label patchi = bMesh.whichPatch(facei);
+            const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(bMesh[patchi]);
 
-            if(cpp.transform() == coupledPolyPatch::transformType::ROTATIONAL)
+            if (cpp)
             {
-                FatalErrorInFunction
-                    << "Rotational cyclic patches are not supported in parallel.\n"
-                    << "Try to decompose the domain so that the rotational cyclic patch "
-                    << "is not split in between processors."
-                    << exit(FatalError);
-            }
-            cpp.neighbPatch().transformPosition(separationVectors[cID], 0);
-        }
+                // Get the neighbor cell across the cyclic face
+                const label cycNeiPatch = cpp->neighbPatchID();
+                const label cycNeiCell =
+                    bMesh[cycNeiPatch].faceCells()[facei - cpp->start()];
+                const List<label>& cycNeiCellFaces = mesh_.cells()[cycNeiCell];
 
-        for(int i = 0; i < 2; i++)
-        {
-            for(int j = 0; j < 2; j++)
-            {
-                for(int k = 0; k < 2; k++)
+                // Loop over all the faces of the neighbor cell
+                for (const label cycNeiCellFace : cycNeiCellFaces)
                 {
-                    vector separation =   i*separationVectors[0]
-                                        + j*separationVectors[1]
-                                        + k*separationVectors[2];
-
-                    scalar testDistance = mag
-                    (
-                        (globalIdxCellCentre - separation)
-                        -
-                        mesh_.C()[celli]
-                    );
-                    if(debug) Info << "testDistance " << testDistance << endl;
-
-                    if( testDistance < distance )
+                    if (mesh_.isInternalFace(cycNeiCellFace))
                     {
-                        distance = testDistance;
-                        patches = List<label>(0);
-                        List<label> applyCyclic({i,j,k});
+                        continue;
+                    }
 
-                        if(debug) Info << "distance " << distance << endl;
-                        if(debug) Info << "separation " << separation << endl;
-                        if(debug) Info << "applyCyclic " << applyCyclic << endl;
+                    // Check if the neighbor cell has processor patches
+                    const label neiPatch = bMesh.whichPatch(cycNeiCellFace);
+                    const processorPolyPatch* ppp =
+                        isA<processorPolyPatch>(bMesh[neiPatch]);
 
-                        for(int n = 0; n < 3; n++)
+                    if (ppp)
+                    {
+                        // Avoid duplicate entries
+                        if (patches.found(cycNeiPatch))
                         {
-                            if(cyclicID[n] != -1 && applyCyclic[n] == 1)
+                            continue;
+                        }
+
+                        // Since we can not access any information on globalIdx
+                        // we use the cell centre map from the stencil to
+                        // identify the cell.
+                        const label localFaceID = cycNeiCellFace - ppp->start();
+                        const vector neiCentre =
+                            ppp->neighbFaceCellCentres()[localFaceID];
+
+                        forAll(cyclicCentres_()[celli], k)
+                        {
+                            if
+                            (
+                                (
+                                    mag(cyclicCentres_()[celli][k] - neiCentre)
+                                        < 100*SMALL
+                                )
+                                && (stencil_[celli][k] == globalIdx)
+                            )
                             {
-                                const cyclicPolyPatch& cpp =
-                                    static_cast<const cyclicPolyPatch&>
-                                    (
-                                        bMesh[cyclicID[n]]
-                                    );
-                                if(debug)
-                                {
-                                    Info << "cpp.name() " << cpp.name() << endl;
-                                }
-                                patches.append(cpp.neighbPatchID());
+                                patches.append(cycNeiPatch);
+                                // Here an alternative might be to append patchi
+                                // and do:
+                                //     cpp.transformPosition()
+                                // instead of
+                                //     cpp.neighbPatch().transformPosition()
+                                // in getValue()
                             }
                         }
                     }

src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H

@@ -186,7 +186,7 @@ public:
 
         //- Finds and returns list of all cyclic patch labels to which celli's
         //  point neighbour cell, globalIdx, belongs. celli and globalIdx touch
-        //  in at least one point on these patches. globalIdx typically belongs
+        //  in at least one point on these patces. globalIdx typically belongs
         //  to stencil_[celli]. The returned label list is used to transform
         //  positions across cyclic boundaries e.g. to be able to calculate
         //  distances between cell centres and interface centres in plicRDF
@@ -194,8 +194,7 @@ public:
         List<label> getCyclicPatches
         (
             const label celli,
-            const label globalIdx,
-            const vector globalIdxCellCentre
+            const label globalIdx
         ) const;
 
         //- Gives patchNumber and patchFaceNumber for a given

src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H

@@ -174,12 +174,7 @@ Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getPositionFields
                 List<label> cyclicPatches(0);
                 if(checkTransformation)
                 {
-                    cyclicPatches = getCyclicPatches
-                    (
-                        celli,
-                        gblIdx,
-                        getValue(phi, neiValues, gblIdx)
-                    );
+                    cyclicPatches = getCyclicPatches(celli, gblIdx);
                 }
 
                 tmpField.append

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

@@ -705,6 +705,8 @@ bool Foam::faceAreaWeightAMI::calculate
         srcCentroids_[i].transfer(srcCtr[i]);
     }
 
+    tgtAddress_.setSize(tgtAddr.size());
+    tgtWeights_.setSize(tgtWght.size());
     forAll(tgtAddr, i)
     {
         tgtAddress_[i].transfer(tgtAddr[i]);

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveUniformityCellZone/objectiveUniformityCellZone.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2025 PCOpt/NTUA
-    Copyright (C) 2013-2025 FOSS GP
+    Copyright (C) 2007-2023 PCOpt/NTUA
+    Copyright (C) 2013-2023 FOSS GP
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -118,11 +118,11 @@ scalar objectiveUniformityCellZone::J()
     const volVectorField& U = vars_.UInst();
     const scalarField& V = mesh_.V().field();
 
-    forAll(zones_, zI)
+    for (const label zI : zones_)
     {
-        const cellZone& cz = mesh_.cellZones()[zones_[zI]];
-        scalarField VZone(V, cz);
-        vectorField UZone(U.primitiveField(), cz);
+        const cellZone& zoneI = mesh_.cellZones()[zI];
+        scalarField VZone(V, zoneI);
+        vectorField UZone(U.primitiveField(), zoneI);
         volZone_[zI] = gSum(VZone);
         UMean_[zI] = gSum(UZone*VZone)/volZone_[zI];
         UVar_[zI] = gSum(magSqr(UZone - UMean_[zI])*VZone)/volZone_[zI];
@@ -138,10 +138,10 @@ void objectiveUniformityCellZone::update_dJdv()
 {
     const volVectorField& U = vars_.U();
 
-    forAll(zones_, zI)
+    for (const label zI : zones_)
     {
-        const cellZone& cz = mesh_.cellZones()[zones_[zI]];
-        for (const label cellI : cz)
+        const cellZone& zoneI = mesh_.cellZones()[zI];
+        for (const label cellI : zoneI)
         {
             dJdvPtr_()[cellI] = (U[cellI] - UMean_[zI])/volZone_[zI];
         }
@@ -154,10 +154,10 @@ void objectiveUniformityCellZone::update_divDxDbMultiplier()
     volScalarField& divDxDbMult = divDxDbMultPtr_();
     const volVectorField& U = vars_.U();
 
-    forAll(zones_, zI)
+    for (const label zI : zones_)
     {
-        const cellZone& cz = mesh_.cellZones()[zones_[zI]];
-        for (const label cellI : cz)
+        const cellZone& zoneI = mesh_.cellZones()[zI];
+        for (const label cellI : zoneI)
         {
             divDxDbMult[cellI] =
                 0.5*(magSqr(U[cellI] - UMean_[zI]) - UVar_[zI])/volZone_[zI];

src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C

@@ -260,7 +260,7 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
     const labelListList& stencil = distribute.getStencil();
 
 
-    forAll(nextToInterface, celli)
+    forAll(nextToInterface,celli)
     {
         if (nextToInterface[celli])
         {
@@ -289,17 +289,7 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                                 centre,
                                 mapCentres,
                                 gblIdx,
-                                distribute.getCyclicPatches
-                                (
-                                    celli,
-                                    gblIdx,
-                                    distribute.getValue
-                                    (
-                                        centre,
-                                        mapCentres,
-                                        gblIdx
-                                    )
-                                )
+                                distribute.getCyclicPatches(celli, gblIdx)
                             )
                         );
                         vector distanceToIntSeg(c - p);
@@ -351,12 +341,7 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                     forAll(stencil[pCellI], j)
                     {
                         const label gblIdx = stencil[pCellI][j];
-                        vector n = -distribute.getValue
-                        (
-                            normal,
-                            mapNormal,
-                            gblIdx
-                        );
+                        vector n = -distribute.getValue(normal, mapNormal, gblIdx);
                         if (mag(n) != 0)
                         {
                             n /= mag(n);
@@ -367,17 +352,7 @@ const Foam::volScalarField&  Foam::reconstructedDistanceFunction::constructRDF
                                     centre,
                                     mapCentres,
                                     gblIdx,
-                                    distribute.getCyclicPatches
-                                    (
-                                        pCellI,
-                                        gblIdx,
-                                        distribute.getValue
-                                        (
-                                            centre,
-                                            mapCentres,
-                                            gblIdx
-                                        )
-                                    )
+                                    distribute.getCyclicPatches(pCellI, gblIdx)
                                 )
                             );
                             vector distanceToIntSeg(c - p);

src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C

@@ -103,12 +103,7 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
                         centre_,
                         mapCentre,
                         gblIdx,
-                        exchangeFields.getCyclicPatches
-                        (
-                            celli,
-                            gblIdx,
-                            exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
-                        )
+                        exchangeFields.getCyclicPatches(celli, gblIdx)
                     )
                 );
                 vector distanceToIntSeg = (tensor::I- n*n) & (p - centre);
@@ -162,12 +157,7 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
                     exchangeFields.getPosition
                     (
                         mesh_.C(), mapCC, gblIdx,
-                        exchangeFields.getCyclicPatches
-                        (
-                            celli,
-                            gblIdx,
-                            exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
-                        )
+                        exchangeFields.getCyclicPatches(celli, gblIdx)
                     )
                 );
                 alphaValues.append
@@ -218,12 +208,7 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
                 exchangeFields.getPosition
                 (
                     mesh_.C(), mapCC, gblIdx,
-                    exchangeFields.getCyclicPatches
-                    (
-                        celli,
-                        gblIdx,
-                        exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
-                    )
+                    exchangeFields.getCyclicPatches(celli, gblIdx)
                 )
             );
             phiValues.append