Parallelism and robustness improvements. Inno Gatin

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C

@@ -181,8 +181,32 @@ void Foam::immersedBoundaryFvPatch::updatePhi
     scalarField sGamma =
         mag(ibPolyPatch_.correctedFaceAreas())/mag(mesh.faceAreas());
 
+    // Scaling of internal mesh flux field should be done only for the current
+    // ib patch to avoid scaling multiple times in case of multiple Ib patches
+    // present. (IG 3/Dec/2018)
+
     // Scale internalField
-    phi.internalField() *= scalarField::subField(sGamma, mesh.nInternalFaces());
+    scalarField& phiIn = phi.internalField();
+
+    const labelList& deadFaces = ibPolyPatch_.deadFaces();
+    forAll(deadFaces, dfI)
+    {
+        const label faceI = deadFaces[dfI];
+        if (mesh.isInternalFace(faceI))
+        {
+            phiIn[faceI] *= sGamma[faceI];
+        }
+    }
+ 
+    const labelList& cutFaces = ibPolyPatch_.ibFaces();
+    forAll(cutFaces, cfI)
+    {
+        const label faceI = cutFaces[cfI];
+        if (mesh.isInternalFace(faceI))
+        {
+            phiIn[faceI] *= sGamma[faceI];
+        }
+    }
 
     // Scale all other patches
     forAll (mesh.boundary(), patchI)
@@ -216,7 +240,6 @@ void Foam::immersedBoundaryFvPatch::updatePhi
     const unallocLabelList& owner = mesh.owner();
     const unallocLabelList& neighbour = mesh.neighbour();
 
-    const scalarField& phiIn = phi.internalField();
     forAll(owner, faceI)
     {
         divPhi[owner[faceI]] += phiIn[faceI];

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/immersedBoundaryFieldBase/immersedBoundaryFieldBase.C

@@ -79,18 +79,75 @@ void Foam::immersedBoundaryFieldBase<Type>::writeField
     const fvPatchField<Type>& f
 ) const
 {
-    // Write VTK on master only
-//    if (Pstream::master())
+    // Write immersed boundary data as a vtk file
+    autoPtr<surfaceWriter> writerPtr = surfaceWriter::New("vtk");
+
+    // Get the intersected patch
+    const standAlonePatch& ts = ibPatch_.ibPolyPatch().ibPatch();
+
+    if (Pstream::parRun())
     {
-        // Add parallel reduction of all faces and data to proc 0
-        // and write the whole patch together.  To Do: HJ, 4/Dec/2018
+        // Gather points, fields and faces
 
-        // Write immersed boundary data as a vtk file
-        autoPtr<surfaceWriter> writerPtr = surfaceWriter::New("vtk");
+        // Points
+        List<pointField> procPoints(Pstream::nProcs());
+        procPoints[Pstream::myProcNo()] = ts.points();
+        Pstream::gatherList(procPoints);
 
-        // Get the intersected patch
-        const standAlonePatch& ts = ibPatch_.ibPolyPatch().ibPatch();
+        // Fields
+        List<Field<Type> > procFields(Pstream::nProcs());
+        procFields[Pstream::myProcNo()] = f;
+        Pstream::gatherList(procFields);
 
+        // Faces
+        List<faceList> procFaces(Pstream::nProcs());
+        procFaces[Pstream::myProcNo()] = ts;
+        Pstream::gatherList(procFaces);
+
+        if (Pstream::master())
+        {
+            // Assemble unique lists to correspond to a single surface
+            pointField allPoints(0);
+            Field<Type> completeField(0);
+            faceList allFaces(0);
+            label prevProcPatchSize = 0;
+
+            forAll(procPoints, procI)
+            {
+                allPoints.append(procPoints[procI]);
+                completeField.append(procFields[procI]);
+
+                // Point labels in faces need to be incremented with respect to
+                // the size of the size of the previous processore patch
+                forAll(procFaces[procI], faceI)
+                {
+                    face curFace = procFaces[procI][faceI];
+                    forAll(curFace, pointI)
+                    {
+                        curFace[pointI] += prevProcPatchSize;
+                    }
+                    allFaces.append(curFace);
+                }
+
+                // Increment the total number of points
+                prevProcPatchSize += procPoints[procI].size();
+            }
+
+            writerPtr->write
+            (
+                f.dimensionedInternalField().path(),
+                ibPatch_.name(),
+                allPoints,
+                allFaces,
+                f.dimensionedInternalField().name(),
+                completeField,
+                false, // FACE_DATA
+                false  // verbose
+            );
+        }
+    }
+    else
+    {
         writerPtr->write
         (
             f.dimensionedInternalField().path(),

src/immersedBoundary/immersedBoundary/immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.C

@@ -1282,6 +1282,8 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
         if
         (
             mag(sumSf[ccc]) > 1e-12
+            // The below criteria is still not firmly established, potential
+            // place for cutting errors (IG 19/Dec/2018)
          && mag(sumSf[ccc] + ibSf[cutCellI])/cutCellVolumes[cutCellI] > 1e-6
         )
         {