diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e5d9f0a3c3..0d4f455a75 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -514,7 +514,7 @@ void Foam::conformalVoronoiMesh::createFeaturePoints()
 
     insertMixedFeaturePoints();
 
-    Info<< "    Inserted " << number_of_vertices() << " vertices" << endl;
+    Pout<< "    Inserted " << number_of_vertices() << " vertices" << endl;
 
     featureVertices_.setSize(number_of_vertices());
 
@@ -563,10 +563,15 @@ void Foam::conformalVoronoiMesh::insertConvexFeaturePoints()
             ptI++
         )
         {
-            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
-
             const Foam::point& featPt = feMesh.points()[ptI];
 
+            if (!geometryToConformTo_.positionOnThisProc(featPt))
+            {
+                continue;
+            }
+
+            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
+
             vector cornerNormal = sum(featPtNormals);
             cornerNormal /= mag(cornerNormal);
 
@@ -610,10 +615,15 @@ void Foam::conformalVoronoiMesh::insertConcaveFeaturePoints()
             ptI++
         )
         {
-            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
-
             const Foam::point& featPt = feMesh.points()[ptI];
 
+            if (!geometryToConformTo_.positionOnThisProc(featPt))
+            {
+                continue;
+            }
+
+            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
+
             vector cornerNormal = sum(featPtNormals);
             cornerNormal /= mag(cornerNormal);
 
@@ -703,6 +713,11 @@ void Foam::conformalVoronoiMesh::insertMixedFeaturePoints()
 
                 const Foam::point& pt(feMesh.points()[ptI]);
 
+                if (!geometryToConformTo_.positionOnThisProc(pt))
+                {
+                    continue;
+                }
+
                 scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
 
                 forAll(pEds, e)
@@ -1162,7 +1177,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
 :
     Delaunay(),
     runTime_(runTime),
-    rndGen_(7864293),
+    rndGen_(64293*Pstream::myProcNo()),
     allGeometry_
     (
         IOobject
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 417eba94e5..118de92f1e 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -416,6 +416,18 @@ private:
             const Delaunay::Finite_vertices_iterator& vit
         ) const;
 
+        //- Check to see if dual cell specified by given vertex iterator
+        //  intersects the parallel interface, hence needs referred.  Returns
+        //  destination processor(s).
+        List<label> parallelInterfaceIntersection
+        (
+            const Delaunay::Finite_vertices_iterator& vit
+        ) const;
+
+        //- Which processor should a vertex generating the given dual edge be
+        //  referred to
+        label targetProc(direction ptOnFace) const;
+
         //- Find the "worst" protrusion of a dual cell through the surface,
         //  subject to the maxSurfaceProtrusion tolerance
         void dualCellLargestSurfaceProtrusion
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 4d0e65a90c..10552a59e4 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -152,6 +152,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                     hitSurface
                 );
 
+                // Pout<< "vert " << vert << endl;
+
                 if (surfHit.hit())
                 {
                     vit->setNearBoundary();
@@ -176,10 +178,12 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         );
                     }
                 }
+
+                // Pout<< "    surfHit " << surfHit << endl;
             }
         }
 
-        Info<< nl << "Initial conformation" << nl
+        Pout<< nl << "Initial conformation" << nl
             << "    Number of vertices " << number_of_vertices() << nl
             << "    Number of surface hits " << surfaceHits.size() << nl
             << "    Number of edge hits " << featureEdgeHits.size()
@@ -204,6 +208,171 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
         initialTotalHits = surfaceHits.size() + featureEdgeHits.size();
     }
 
+    label totalInterfacePoints = 0;
+
+    label nIter = 0;
+
+    if (Pstream::parRun())
+    {
+        do
+        {
+            // Propagating vertices to other processors to form halo information
+
+            DynamicList<Foam::point> parallelInterfacePoints;
+            DynamicList<label> targetProcessor;
+
+            for
+            (
+                Delaunay::Finite_vertices_iterator vit =
+                    finite_vertices_begin();
+                vit != finite_vertices_end();
+                vit++
+            )
+            {
+                if (vit->internalOrBoundaryPoint())
+                {
+                    List<label> toProcs = parallelInterfaceIntersection(vit);
+
+                    forAll(toProcs, tPI)
+                    {
+                        label toProc = toProcs[tPI];
+
+                        if (toProc > -1)
+                        {
+                            parallelInterfacePoints.append
+                            (
+                                topoint(vit->point())
+                            );
+
+                            targetProcessor.append(toProc);
+
+                            std::list<Vertex_handle> incidentVertices;
+
+                            incident_vertices
+                            (
+                                vit,
+                                std::back_inserter(incidentVertices)
+                            );
+
+                            for
+                            (
+                                std::list<Vertex_handle>::iterator ivit =
+                                    incidentVertices.begin();
+                                ivit != incidentVertices.end();
+                                ++ivit
+                            )
+                            {
+                                parallelInterfacePoints.append
+                                (
+                                    topoint((*ivit)->point())
+                                );
+
+                                targetProcessor.append(toProc);
+                            }
+                        }
+                    }
+                }
+            }
+
+            writePoints("parallelInterfacePoints.obj", parallelInterfacePoints);
+
+            // Determine send map
+            // ~~~~~~~~~~~~~~~~~~
+
+            // 1. Count
+            labelList nSend(Pstream::nProcs(), 0);
+
+            forAll(targetProcessor, i)
+            {
+                label procI = targetProcessor[i];
+
+                nSend[procI]++;
+            }
+
+            // Send over how many I need to receive
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            labelListList sendSizes(Pstream::nProcs());
+
+            sendSizes[Pstream::myProcNo()] = nSend;
+
+            combineReduce(sendSizes, UPstream::listEq());
+
+            // 2. Size sendMap
+            labelListList sendMap(Pstream::nProcs());
+
+            forAll(nSend, procI)
+            {
+                sendMap[procI].setSize(nSend[procI]);
+
+                nSend[procI] = 0;
+            }
+
+            // 3. Fill sendMap
+            forAll(targetProcessor, i)
+            {
+                label procI = targetProcessor[i];
+
+                sendMap[procI][nSend[procI]++] = i;
+            }
+
+            // Determine receive map
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            labelListList constructMap(Pstream::nProcs());
+
+            // Local transfers first
+            constructMap[Pstream::myProcNo()] = identity
+            (
+                sendMap[Pstream::myProcNo()].size()
+            );
+
+            label constructSize = constructMap[Pstream::myProcNo()].size();
+
+            forAll(constructMap, procI)
+            {
+                if (procI != Pstream::myProcNo())
+                {
+                    label nRecv = sendSizes[procI][Pstream::myProcNo()];
+
+                    constructMap[procI].setSize(nRecv);
+
+                    for (label i = 0; i < nRecv; i++)
+                    {
+                        constructMap[procI][i] = constructSize++;
+                    }
+                }
+            }
+
+            mapDistribute pointMap
+            (
+                constructSize,
+                sendMap.xfer(),
+                constructMap.xfer()
+            );
+
+            totalInterfacePoints = parallelInterfacePoints.size();
+
+            reduce(totalInterfacePoints, sumOp<label>());
+
+            pointMap.distribute(parallelInterfacePoints);
+
+            forAll(parallelInterfacePoints, ptI)
+            {
+                insertPoint
+                (
+                    parallelInterfacePoints[ptI],
+                    Vb::ptFarPoint
+                );
+            }
+
+            Info<< "totalInterfacePoints " << totalInterfacePoints << endl;
+
+            nIter++;
+
+        } while (totalInterfacePoints > 0 && nIter < 3);
+    }
+
     label iterationNo = 0;
 
     label maxIterations =
@@ -214,7 +383,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
 
     label hitLimit = label(iterationToIntialHitRatioLimit*initialTotalHits);
 
-    Info<< nl << "Stopping iterations when: " << nl
+    Pout<< nl << "Stopping iterations when: " << nl
         <<"    total number of hits drops below "
         << iterationToIntialHitRatioLimit << " of initial hits ("
         << hitLimit << ")" << nl
@@ -469,7 +638,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
     //     }
     // }
 
-    storeSurfaceConformation();
+    Pout<< "NOT STORING SURFACE CONFORMATION" << endl;
+    //storeSurfaceConformation();
 }
 
 
@@ -498,21 +668,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
         }
 
         Foam::point dE0 = topoint(dual(fit->first));
-
-        // If edge end is outside bounding box then edge cuts boundary
-        if (!geometryToConformTo_.bounds().contains(dE0))
-        {
-            return true;
-        }
-
         Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
 
-        // If other edge end is outside bounding box then edge cuts boundary
-        if (!geometryToConformTo_.bounds().contains(dE1))
-        {
-            return true;
-        }
-
         // Check for the edge passing through a surface
         if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
         {
@@ -524,6 +681,115 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
 }
 
 
+Foam::List<Foam::label>
+Foam::conformalVoronoiMesh::parallelInterfaceIntersection
+(
+    const Delaunay::Finite_vertices_iterator& vit
+) const
+{
+    std::list<Facet> facets;
+    incident_facets(vit, std::back_inserter(facets));
+
+    DynamicList<label> procs;
+
+    for
+    (
+        std::list<Facet>::iterator fit=facets.begin();
+        fit != facets.end();
+        ++fit
+    )
+    {
+        if
+        (
+            is_infinite(fit->first)
+         || is_infinite(fit->first->neighbor(fit->second))
+        )
+        {
+            return procs;
+        }
+
+        Foam::point dE0 = topoint(dual(fit->first));
+        Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
+
+        Foam::point boxPt(vector::one*GREAT);
+        direction ptOnFace = -1;
+
+        geometryToConformTo_.bounds().intersects
+        (
+            dE0,
+            dE1 - dE0,
+            dE0,
+            dE1,
+            boxPt,
+            ptOnFace
+        );
+
+        if
+        (
+            ptOnFace > 0
+         && !geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1)
+        )
+        {
+            // A surface penetration is not found, and a bounds penetration is
+            // found, so the dual edge has crossed the parallel interface.
+
+            label target = targetProc(ptOnFace);
+
+            if
+            (
+                target > -1
+             && findIndex(procs, target) == -1
+            )
+            {
+                procs.append(target);
+            }
+
+            // Pout<< "Parallel box penetration, face " << ptOnFace
+            //     << " proc " << target << endl;
+            // meshTools::writeOBJ(Pout, dE0);
+            // meshTools::writeOBJ(Pout, dE1);
+            // meshTools::writeOBJ(Pout, boxPt);
+            // Pout << "l 1 2" << endl;
+        }
+    }
+
+    return procs;
+}
+
+
+Foam::label Foam::conformalVoronoiMesh::targetProc
+(
+     direction ptOnFace
+) const
+{
+    // Hard coded to 8 proc, bound-box octant split:
+
+    label faceIndex = 0;
+
+    // From: http://graphics.stanford.edu/~seander/bithacks.html
+    while (ptOnFace >>= 1)
+    {
+        faceIndex++;
+    }
+
+    label procArray[8][6] =
+    {
+        {-1,  1, -1,  2, -1, 4},
+        { 0, -1, -1,  3, -1, 5},
+        {-1,  3,  0, -1, -1, 6},
+        { 2, -1,  1, -1, -1, 7},
+        {-1,  5, -1,  6,  0, -1},
+        { 4, -1, -1,  7,  1, -1},
+        {-1,  7,  4, -1,  2, -1},
+        { 6, -1,  5, -1,  3, -1}
+    };
+
+    return initListList<List<label>, label, 8, 6>(procArray)
+        [Pstream::myProcNo()][faceIndex];
+}
+
+
+
 void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
 (
     const Delaunay::Finite_vertices_iterator& vit,
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
index 3322afc077..928b964c87 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
@@ -94,6 +94,11 @@ bool Foam::conformalVoronoiMesh::insertSpecialisedFeaturePoint
 
     if (nExternal == 2 && nInternal == 1)
     {
+        // if (!geometryToConformTo_.positionOnThisProc(pt))
+        // {
+        //     continue;
+        // }
+
         // Info<< "nExternal == 2 && nInternal == 1" << endl;
 
         // const Foam::point& featPt = feMesh.points()[ptI];
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index ca0c67a21c..8580338991 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -240,14 +240,14 @@ inline void Foam::conformalVoronoiMesh::insertPoints
     const std::vector<Point>& points
 )
 {
-    Info<< "    " << points.size() << " points to insert..." << endl;
+    Pout<< "    " << points.size() << " points to insert..." << endl;
 
     label nVert = number_of_vertices();
 
     // using the range insert (faster than inserting points one by one)
     insert(points.begin(), points.end());
 
-    Info<< "    " << number_of_vertices() - startOfInternalPoints_
+    Pout<< "    " << number_of_vertices() - startOfInternalPoints_
         << " points inserted" << endl;
 
     for
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 91f0581291..a7480fe7a2 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -150,7 +150,7 @@ void Foam::conformalVoronoiMesh::writeInternalDelaunayVertices
         internalDelaunayVertices
     );
 
-    Info<< nl
+    Pout<< nl
         << "Writing " << internalDVs.name()
         << " to " << internalDVs.instance()
         << endl;
@@ -228,7 +228,7 @@ void Foam::conformalVoronoiMesh::writeMesh
             filterFaces
         );
 
-        Info<< nl << "Writing polyMesh to " << instance << endl;
+        Pout<< nl << "Writing polyMesh to " << instance << endl;
 
         writeMesh
         (
@@ -572,7 +572,7 @@ void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
 
     if (!protrudingCells.empty())
     {
-        Info<< nl << "Found " << protrudingCells.size()
+        Pout<< nl << "Found " << protrudingCells.size()
             << " cells protruding from the surface, writing cellSet "
             << protrudingCells.name()
             << endl;
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 870bd53658..3c8477bbf5 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -59,6 +59,7 @@ class indexedVertex
         //  ptNearBoundaryPoint  : internal near boundary point.
         //  ptInternalPoint      : internal point.
         //  ptFarPoint           : far-point.
+        //  > ptFarPoint, < 0    : referred point from processor -(type_ + 1)
         //  >= 0                 : part of point-pair. Index of other point.
         //                         Lowest numbered is inside one (master).
         int type_;
@@ -72,10 +73,9 @@ public:
 
     enum pointTypes
     {
-        ptNearBoundaryPoint  = -4,
-        ptInternalPoint      = -3,
-        ptMirrorPoint        = -2,
-        ptFarPoint           = -1
+        ptNearBoundaryPoint  = INT_MIN,
+        ptInternalPoint      = INT_MIN + 1,
+        ptFarPoint           = INT_MIN + 2
     };
 
     typedef typename Vb::Vertex_handle      Vertex_handle;
@@ -198,6 +198,27 @@ public:
     }
 
 
+    // is this a referred vertex
+    inline bool referred() const
+    {
+        return type_ < 0 && type_ > ptFarPoint;
+    }
+
+
+    // For referred vertices, what is the processor index
+    inline int procIndex() const
+    {
+        if (referred())
+        {
+            return -(type_ + 1);
+        }
+        else
+        {
+            return -1;
+        }
+    }
+
+
     //- Set the point to be internal
     inline void setInternal()
     {
@@ -218,12 +239,6 @@ public:
         type_ = ptNearBoundaryPoint;
     }
 
-    //- Is point a mirror point
-    inline bool mirrorPoint() const
-    {
-        return type_ == ptMirrorPoint;
-    }
-
 
     //- Either master or slave of pointPair.
     inline bool pairPoint() const
@@ -270,7 +285,7 @@ public:
     //- Is point near the boundary or part of the boundary definition
     inline bool nearOrOnBoundary() const
     {
-        return pairPoint() || mirrorPoint() || nearBoundary();
+        return pairPoint() || nearBoundary();
     }
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 13d60d4463..32d3e2832f 100644
--- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -361,6 +361,15 @@ bool Foam::conformationSurfaces::overlaps(const treeBoundBox& bb) const
 }
 
 
+bool Foam::conformationSurfaces::positionOnThisProc(const point& pt) const
+{
+    // This is likely to give problems when a point is on the boundary between
+    // two processors.
+
+    return bounds_.contains(pt);
+}
+
+
 Foam::Field<bool> Foam::conformationSurfaces::inside
 (
     const pointField& samplePts
diff --git a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index 9b9fd9b9ab..6271ddcfc4 100644
--- a/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/src/mesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -170,6 +170,9 @@ public:
             //  the surfaces
             bool overlaps(const treeBoundBox& bb) const;
 
+            //- Check if the point is in the domain handled by this processor
+            bool positionOnThisProc(const point& pt) const;
+
             //- Check if points are inside surfaces to conform to
             Field<bool> inside(const pointField& samplePts) const;
 
diff --git a/src/mesh/conformalVoronoiMesh/initialPointsMethod/hierarchicalDensityWeightedStochastic/hierarchicalDensityWeightedStochastic.C b/src/mesh/conformalVoronoiMesh/initialPointsMethod/hierarchicalDensityWeightedStochastic/hierarchicalDensityWeightedStochastic.C
index 53e1284810..8dcbf2c2ff 100644
--- a/src/mesh/conformalVoronoiMesh/initialPointsMethod/hierarchicalDensityWeightedStochastic/hierarchicalDensityWeightedStochastic.C
+++ b/src/mesh/conformalVoronoiMesh/initialPointsMethod/hierarchicalDensityWeightedStochastic/hierarchicalDensityWeightedStochastic.C
@@ -119,7 +119,7 @@ void Foam::hierarchicalDensityWeightedStochastic::recurseAndFill
 
                 if (debug)
                 {
-                    Info<< newName + "_overlap " << subBB << endl;
+                    Pout<< newName + "_overlap " << subBB << endl;
                 }
 
                 if (!fillBox(initialPoints, subBB, true))
@@ -144,7 +144,7 @@ void Foam::hierarchicalDensityWeightedStochastic::recurseAndFill
 
             if (debug)
             {
-                Info<< newName + "_inside " << subBB << endl;
+                Pout<< newName + "_inside " << subBB << endl;
             }
 
             if (!fillBox(initialPoints, subBB, false))
@@ -221,7 +221,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
         if (surfHit.hit())
         {
-            // Info<< "box wellInside, no need to sample surface." << endl;
+            // Pout<< "box wellInside, no need to sample surface." << endl;
 
             wellInside = true;
         }
@@ -248,7 +248,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
             minimumSurfaceDistanceCoeffSqr_*sqr(cornerSizes)
         );
 
-        // Info<< corners << nl << cornerSizes << nl << insideCorners << endl;
+        // Pout<< corners << nl << cornerSizes << nl << insideCorners << endl;
 
         forAll(insideCorners, i)
         {
@@ -267,7 +267,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
             if (maxCellSize/minCellSize > maxSizeRatio_)
             {
-                // Info<< "Abort fill at corner sample stage,"
+                // Pout<< "Abort fill at corner sample stage,"
                 //     << " minCellSize " << minCellSize
                 //     << " maxCellSize " << maxCellSize
                 //     << " maxSizeRatio " << maxCellSize/minCellSize
@@ -281,7 +281,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
                 // If one or more corners is not "wellInside", then treat this
                 // as an overlapping box.
 
-                // Info<< "Inside box found to have some non-wellInside "
+                // Pout<< "Inside box found to have some non-wellInside "
                 //     << "corners, using overlapping fill."
                 //     << endl;
 
@@ -378,7 +378,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
                     if (maxCellSize/minCellSize > maxSizeRatio_)
                     {
-                        // Info<< "Abort fill at surface sample stage, "
+                        // Pout<< "Abort fill at surface sample stage, "
                         //     << " minCellSize " << minCellSize
                         //     << " maxCellSize " << maxCellSize
                         //     << " maxSizeRatio " << maxCellSize/minCellSize
@@ -393,7 +393,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
                         // then treat this as an overlapping box.
                         overlapping = true;
 
-                        // Info<< "Inside box found to have some non-"
+                        // Pout<< "Inside box found to have some non-"
                         //     << "wellInside surface points, using "
                         //     << "overlapping fill."
                         //     << endl;
@@ -480,7 +480,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
                 if (maxCellSize/minCellSize > maxSizeRatio_)
                 {
-                    // Info<< "Abort fill at sample stage,"
+                    // Pout<< "Abort fill at sample stage,"
                     //     << " minCellSize " << minCellSize
                     //     << " maxCellSize " << maxCellSize
                     //     << " maxSizeRatio " << maxCellSize/minCellSize
@@ -493,17 +493,17 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
         if (nInside == 0)
         {
-            // Info<< "No sample points found inside box" << endl;
+            // Pout<< "No sample points found inside box" << endl;
 
             return true;
         }
 
-        // Info<< scalar(nInside)/scalar(samplePoints.size())
+        // Pout<< scalar(nInside)/scalar(samplePoints.size())
         //     << " full overlapping box" << endl;
 
         totalVolume *= scalar(nInside)/scalar(samplePoints.size());
 
-        // Info<< "Total volume to fill = " << totalVolume << endl;
+        // Pout<< "Total volume to fill = " << totalVolume << endl;
 
         // Using the sampledPoints as the first test locations as they are
         // randomly shuffled, but unfiormly sampling space and have wellInside
@@ -547,7 +547,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
                         scalar r = rnd.scalar01();
 
-                        // Info<< "totalVolume " << totalVolume << nl
+                        // Pout<< "totalVolume " << totalVolume << nl
                         //     << "volumeAdded " << volumeAdded << nl
                         //     << "localVolume " << localVolume << nl
                         //     << "addProbability " << addProbability << nl
@@ -559,7 +559,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
                             // Place this volume before finishing filling this
                             // box
 
-                            // Info<< "Final volume probability break accept"
+                            // Pout<< "Final volume probability break accept"
                             //     << endl;
 
                             initialPoints.push_back
@@ -583,7 +583,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
     if (volumeAdded < totalVolume)
     {
-        // Info<< "Adding random points, remaining volume "
+        // Pout<< "Adding random points, remaining volume "
         //     << totalVolume - volumeAdded
         //     << endl;
 
@@ -633,7 +633,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
                 if (maxCellSize/minCellSize > maxSizeRatio_)
                 {
-                    // Info<< "Abort fill at random fill stage,"
+                    // Pout<< "Abort fill at random fill stage,"
                     //     << " minCellSize " << minCellSize
                     //     << " maxCellSize " << maxCellSize
                     //     << " maxSizeRatio " << maxCellSize/minCellSize
@@ -665,7 +665,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
                         scalar r = rnd.scalar01();
 
-                        // Info<< "totalVolume " << totalVolume << nl
+                        // Pout<< "totalVolume " << totalVolume << nl
                         //     << "volumeAdded " << volumeAdded << nl
                         //     << "localVolume " << localVolume << nl
                         //     << "addProbability " << addProbability << nl
@@ -677,7 +677,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
                             // Place this volume before finishing filling this
                             // box
 
-                            // Info<< "Final volume probability break accept"
+                            // Pout<< "Final volume probability break accept"
                             //     << endl;
 
                             initialPoints.push_back
@@ -701,7 +701,7 @@ bool Foam::hierarchicalDensityWeightedStochastic::fillBox
 
     globalTrialPoints_ += trialPoints;
 
-    // Info<< trialPoints
+    // Pout<< trialPoints
     //     << " locations queried, " << initialPoints.size() - initialSize
     //     << " points placed, ("
     //     << scalar(initialPoints.size() - initialSize)
@@ -777,7 +777,7 @@ hierarchicalDensityWeightedStochastic::initialPoints() const
 
     if (debug)
     {
-        Info<< "    Filling box " << hierBB << endl;
+        Pout<< "    Filling box " << hierBB << endl;
     }
 
     recurseAndFill
@@ -788,10 +788,18 @@ hierarchicalDensityWeightedStochastic::initialPoints() const
         "recursionBox"
     );
 
-    Info<< "        " << initialPoints.size() << " points placed" << nl
+    label nInitialPoints = initialPoints.size();
+
+    if (Pstream::parRun())
+    {
+        reduce(nInitialPoints, sumOp<label>());
+        reduce(globalTrialPoints_, sumOp<label>());
+    }
+
+    Info<< "        " << nInitialPoints << " points placed" << nl
         << "        " << globalTrialPoints_ << " locations queried" << nl
         << "        "
-        << scalar(initialPoints.size())/scalar(max(globalTrialPoints_, 1))
+        << scalar(nInitialPoints)/scalar(max(globalTrialPoints_, 1))
         << " success rate"
         << endl;
 
diff --git a/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C b/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
index 8dfcaad9b0..f914cad1b3 100644
--- a/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
+++ b/src/mesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
@@ -74,6 +74,18 @@ std::vector<Vb::Point> pointFile::initialPoints() const
             << exit(FatalError) << endl;
     }
 
+    boolList procPt(points.size(), false);
+
+    forAll(points, ptI)
+    {
+        procPt[ptI] = cvMesh_.geometryToConformTo().positionOnThisProc
+        (
+            points[ptI]
+        );
+    }
+
+    inplaceSubset(procPt, points);
+
     std::vector<Vb::Point> initialPoints;
 
     Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside