Feature: comms update for immersed boundary solver

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C

@@ -1095,21 +1095,11 @@ void Foam::immersedBoundaryFvPatch::makeIbCellCells() const
     ibCellCellsPtr_ = new labelListList(ibc.size());
     labelListList& cellCells = *ibCellCellsPtr_;
 
-    ibProcCentresPtr_ = new FieldField<Field, vector>(Pstream::nProcs());
+    ibProcCentresPtr_ = new vectorListList(Pstream::nProcs());
-    FieldField<Field, vector>& procCentres = *ibProcCentresPtr_;
+    vectorListList& procCentres = *ibProcCentresPtr_;
 
-    forAll (procCentres, procI)
+    ibProcGammaPtr_ = new scalarListList(Pstream::nProcs());
-    {
+    scalarListList& procGamma = *ibProcGammaPtr_;
-        procCentres.set(procI, new vectorField(0));
-    }
-
-    ibProcGammaPtr_ = new FieldField<Field, scalar>(Pstream::nProcs());
-    FieldField<Field, scalar>& procGamma = *ibProcGammaPtr_;
-
-    forAll (procGamma, procI)
-    {
-        procGamma.set(procI, new scalarField(0));
-    }
 
 
     ibCellProcCellsPtr_ = new List<List<labelPair> >(ibc.size());
@@ -1250,113 +1240,29 @@ void Foam::immersedBoundaryFvPatch::makeIbCellCells() const
         labelList procIbCells = procIbCellsSet.toc();
         sort(procIbCells);
 
-        // Note: consider more sophisticated gather-scatter
+        // Note: new gather-scatter operations
-        // HJ, 18/Jun/2015
+        // HJ, 11/Aug/2016
-
-        // Send and receive number  of immersed boundary cells
-        // next to processor boundaries
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            if (procI != Pstream::myProcNo())
-            {
-                // Parallel data exchange
-                {
-                    OPstream toProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        sizeof(label)
-                    );
-
-                    toProc << procIbCells.size();
-                }
-            }
-        }
-
-        labelList sizes(Pstream::nProcs(), 0);
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            if (procI != Pstream::myProcNo())
-            {
-                // Parallel data exchange
-                {
-                    IPstream fromProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        sizeof(label)
-                    );
-
-                    fromProc >> sizes[procI];
-                }
-            }
-        }
 
         // Send and receive ibc centres and radii
-        vectorField centres(procIbCells.size(), vector::zero);
+        vectorListList ctrs(Pstream::nProcs());
+
+        ctrs[Pstream::myProcNo()].setSize(procIbCells.size());
+        vectorList& centres = ctrs[Pstream::myProcNo()];
+
         forAll (centres, cellI)
         {
             centres[cellI] = C[ibc[procIbCells[cellI]]];
         }
-        scalarField procRMax(rM, procIbCells);
 
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        Pstream::gatherList(ctrs);
-        {
+        Pstream::scatterList(ctrs);
-            if (procI != Pstream::myProcNo())
-            {
-                // Parallel data exchange
-                {
-                    OPstream toProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        centres.size()*sizeof(vector)
-                      + procRMax.size()*sizeof(scalar)
-                    );
 
-                    toProc << centres << procRMax;
+        scalarListList rMax(Pstream::nProcs());
-                }
-            }
-        }
 
-        FieldField<Field, vector> ctrs(Pstream::nProcs());
+        rMax[Pstream::myProcNo()] = scalarField(rM, procIbCells);
-        FieldField<Field, scalar> rMax(Pstream::nProcs());
 
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        Pstream::gatherList(rMax);
-        {
+        Pstream::scatterList(rMax);
-            if (procI != Pstream::myProcNo())
-            {
-                ctrs.set
-                (
-                    procI,
-                    new vectorField(sizes[procI], vector::zero)
-                );
-
-                rMax.set
-                (
-                    procI,
-                    new scalarField(sizes[procI], 0)
-                );
-
-                // Parallel data exchange
-                {
-                    IPstream fromProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        sizes[procI]*sizeof(vector)
-                      + sizes[procI]*sizeof(scalar)
-                    );
-
-                    fromProc >> ctrs[procI] >> rMax[procI];
-                }
-            }
-            else
-            {
-                ctrs.set(procI, new vectorField(0));
-                rMax.set(procI, new scalarField(0));
-            }
-        }
 
         // Find cells needed by other processors
         if (ibProcCellsPtr_)
@@ -1431,131 +1337,27 @@ void Foam::immersedBoundaryFvPatch::makeIbCellCells() const
             procCells[procI] = procCellSet.toc();
         }
 
+        Pstream::gatherList(procCells);
+        Pstream::scatterList(procCells);
 
-        // Send and receive sizes
+        procCentres[Pstream::myProcNo()] =
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
+            vectorField
-        {
+            (
-            if (procI != Pstream::myProcNo())
+                C,
-            {
+                procCells[Pstream::myProcNo()]
-                // Parallel data exchange
+            );
-                {
-                    OPstream toProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        sizeof(label)
-                    );
 
-                    toProc << procCells[procI].size();
+        Pstream::gatherList(procCentres);
-                }
+        Pstream::scatterList(procCentres);
-            }
-        }
-
-        labelList procSizes(Pstream::nProcs(), 0);
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            if (procI != Pstream::myProcNo())
-            {
-                // Parallel data exchange
-                {
-                    IPstream fromProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        sizeof(label)
-                    );
-
-                    fromProc >> procSizes[procI];
-                }
-            }
-        }
-
-        // Send cell centres
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            if (procI != Pstream::myProcNo())
-            {
-                vectorField centres(C, procCells[procI]);
-
-                // Parallel data exchange
-                {
-                    OPstream toProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        centres.size()*sizeof(vector)
-                    );
-
-                    toProc << centres;
-                }
-            }
-        }
-
-        // Receive cell centres
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            if (procI != Pstream::myProcNo())
-            {
-                procCentres[procI].setSize(procSizes[procI]);
-
-                // Parallel data exchange
-                {
-                    IPstream fromProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        procSizes[procI]*sizeof(vector)
-                    );
-
-                    fromProc >> procCentres[procI];
-                }
-            }
-            // else: already set to zero-size field
-        }
 
         // Send cell gamma
-        const scalarField& gammaI = gamma().internalField();
+        procGamma[Pstream::myProcNo()] =
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
+            scalarField
-        {
+            (
-            if (procI != Pstream::myProcNo())
+                gamma().internalField(),
-            {
+                procCells[Pstream::myProcNo()]
-                scalarField gamma(gammaI, procCells[procI]);
+            );
-
+        // Reset own size to zero?  HJ, 11/Aug/2016
-                // Parallel data exchange
-                {
-                    OPstream toProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        gamma.size()*sizeof(scalar)
-                    );
-
-                    toProc << gamma;
-                }
-            }
-        }
-
-        // Receive cell gamma
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            if (procI != Pstream::myProcNo())
-            {
-                procGamma[procI].setSize(procSizes[procI]);
-
-                // Parallel data exchange
-                {
-                    IPstream fromProc
-                    (
-                        Pstream::blocking,
-                        procI,
-                        procSizes[procI]*sizeof(scalar)
-                    );
-
-                    fromProc >> procGamma[procI];
-                }
-            }
-            // else: already set to zero-size field
-        }
 
         // Cell-procCells addressing
         forAll (cellProcCells, cellI)
@@ -2410,7 +2212,7 @@ const Foam::labelListList& Foam::immersedBoundaryFvPatch::ibCellCells() const
 }
 
 
-const Foam::FieldField<Foam::Field, Foam::vector>&
+const Foam::vectorListList&
 Foam::immersedBoundaryFvPatch::ibProcCentres() const
 {
     if (!ibProcCentresPtr_)
@@ -2422,7 +2224,7 @@ Foam::immersedBoundaryFvPatch::ibProcCentres() const
 }
 
 
-const Foam::FieldField<Foam::Field, Foam::scalar>&
+const Foam::scalarListList&
 Foam::immersedBoundaryFvPatch::ibProcGamma() const
 {
     if (!ibProcGammaPtr_)

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.H

@@ -49,6 +49,8 @@ SourceFiles
 #include "surfaceFieldsFwd.H"
 #include "dynamicLabelList.H"
 #include "labelPair.H"
+#include "scalarList.H"
+#include "vectorList.H"
 #include "FieldFields.H"
 #include "scalarMatrices.H"
 
@@ -70,7 +72,7 @@ class immersedBoundaryFvPatch
 {
     // Private data
 
-        //- Reference to processor patch
+        //- Reference to immersed boundary patch
         const immersedBoundaryPolyPatch& ibPolyPatch_;
 
         //- Finite volume mesh reference
@@ -168,10 +170,10 @@ class immersedBoundaryFvPatch
             mutable labelListList* ibProcCellsPtr_;
 
             //- Centres of cells from neighbour processors
-            mutable FieldField<Field, vector>* ibProcCentresPtr_;
+            mutable vectorListList* ibProcCentresPtr_;
 
             //- Gamma of cells from neighbour processors
-            mutable FieldField<Field, scalar>* ibProcGammaPtr_;
+            mutable scalarListList* ibProcGammaPtr_;
 
             //- Cell-proc-cell addressing
             mutable List<List<labelPair> >* ibCellProcCellsPtr_;
@@ -447,12 +449,12 @@ public:
             //- Return neighbour proc centres
             //  These are centres from neighbouring processors the
             //  local processor needs to receive
-            const FieldField<Field, vector>& ibProcCentres() const;
+            const vectorListList& ibProcCentres() const;
 
             //- Return neighbour proc gamma
             //  These are gamma values from neighbouring processors the
             //  local processor needs to receive
-            const FieldField<Field, scalar>& ibProcGamma() const;
+            const scalarListList& ibProcGamma() const;
 
             //- Return neighbour cell addressing
             const List<List<labelPair> >& ibCellProcCells() const;

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchLeastSquaresFit.C

@@ -66,7 +66,7 @@ void Foam::immersedBoundaryFvPatch::makeInvDirichletMatrices() const
     // Initialize maxRowSum for debug
     scalar maxRowSum = 0.0;
 
-    const FieldField<Field, vector>& procC = ibProcCentres();
+    const vectorListList& procC = ibProcCentres();
 
     label nCoeffs = 5;
 
@@ -305,7 +305,7 @@ void Foam::immersedBoundaryFvPatch::makeInvNeumannMatrices() const
     // Initialize maxRowSum for debug
     scalar maxRowSum = 0.0;
 
-    const FieldField<Field, vector>& procC = ibProcCentres();
+    const vectorListList& procC = ibProcCentres();
 
     label nCoeffs = 6;
 

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchSamplingWeights.C

@@ -74,8 +74,8 @@ void Foam::immersedBoundaryFvPatch::makeIbSamplingWeights() const
     const scalarField& gammaIn = gamma().internalField();
     const vectorField& CIn = mesh_.C().internalField();
 
-    const FieldField<Field, scalar>& gammaProc = ibProcGamma();
+    const scalarListList& gammaProc = ibProcGamma();
-    const FieldField<Field, vector>& CProc = ibProcCentres();
+    const vectorListList& CProc = ibProcCentres();
 
     // Go through all cellCells and calculate inverse distance for
     // all live points

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchTemplates.C

@@ -41,6 +41,9 @@ Foam::immersedBoundaryFvPatch::sendAndReceive
     );
     FieldField<Field, Type>& procPsi = tprocPsi();
 
+    // This requires a rewrite useng mapDistribute
+    // HJ, 11/Aug/2016
+
     forAll (procPsi, procI)
     {
         procPsi.set

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchTriAddressing.C

@@ -61,7 +61,7 @@ void Foam::immersedBoundaryFvPatch::makeTriAddressing() const
     {
         hitTris[hf[hfI]] = hfI;
     }
-
+    Info<< "triangles: " << triPatch.size() << " hit: " << hf.size() << endl;
     // Allocate storage
     cellsToTriAddrPtr_ = new labelListList(triPatch.size());
     labelListList& addr = *cellsToTriAddrPtr_;
@@ -84,8 +84,16 @@ void Foam::immersedBoundaryFvPatch::makeTriAddressing() const
     label faceIndex = 0;
     label counter = 0;
 
-    forAll (triPatch, triI)
+    boolList visited(triPatch.size(), false);
+
+    register label curTri;
+
+    // Only search for tri faces in the mesh
+    forAll (triFacesInMesh, tfimI)
     {
+
+        const label triI = triFacesInMesh[tfimI];
+
         if (hitTris[triI] > -1)
         {
             // Triangle contains IB point
@@ -99,8 +107,8 @@ void Foam::immersedBoundaryFvPatch::makeTriAddressing() const
         {
             // No direct hit.  Start a neighbourhood search
 
-            // Record already visited faces
+            // Reset visited faces
-            labelHashSet visited;
+            visited = false;
 
             // Collect new faces to visit
             SLList<label> nextToVisit;
@@ -113,47 +121,53 @@ void Foam::immersedBoundaryFvPatch::makeTriAddressing() const
 
             do
             {
-                const label curTri = nextToVisit.removeHead();
+                // Pick next face that was not visited by skipping
+                // already visited faces
+                do
+                {
+                    curTri = nextToVisit.removeHead();
+                }
+                while (visited[curTri]);
+
                 // Discard tri if already visited
-                if (visited[curTri])
+                if (visited[curTri]) continue;
-                {
+
-                    continue;
+                visited.insert(curTri);
-                }
-                else
-                {
-                    visited.insert(curTri);
-                }
 
                 const triFace& curTriPoints = triPatch[curTri];
 
                 // For all current points of face, pick up neighbouring faces
                 forAll (curTriPoints, tpI)
                 {
-                    const labelList curNbrs = pf[curTriPoints[tpI]];
+                    const labelList& curNbrs = pf[curTriPoints[tpI]];
 
                     forAll (curNbrs, nbrI)
                     {
-                        if (!visited.found(curNbrs[nbrI]))
+                        if (visited[curNbrs[nbrI]])
+                        {
+                            continue;
+                        }
+                        else
                         {
-                            // Found a face which is not visited.  Add it to
-                            // the list of faces to visit
-                            nextToVisit.append(curNbrs[nbrI]);
-
                             if (hitTris[curNbrs[nbrI]] > -1)
                             {
                                 // Found a neighbour with a hit: use this
                                 // IB point
                                 ibPointsToUse.insert(hitTris[curNbrs[nbrI]]);
                             }
+
+                            // Found a face which is not visited.  Add it to
+                            // the list of faces to visit
+                            nextToVisit.append(curNbrs[nbrI]);
                         }
                     }
                 }
 
-                // If the search has gone wrong, escape with
+                // If the search has gone wrong eg. because of a discrepancy
+                // in the resolution between the mesh and the STL, escape with
                 // poorer interpolation.
                 if (nextToVisit.size() > 200 && !ibPointsToUse.empty())
                 {
-                    Info<< "ESCAPE " << ibPointsToUse.size() << endl;
                     break;
                 }
             } while

tutorials/immersedBoundary/sphereInChannel/system/controlDict

@@ -51,7 +51,6 @@ maxCo           0.2;
 
 libs
 (
-    "liblduSoLvers.so"
     "libimmersedBoundary.so"
 );