diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index c37e9c6f833827502e393000d1595b8a050bde0a..107c732bbfaf5f6e880f9ca933a0044fe4e6f59b 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -989,6 +989,14 @@ public:
                 const UList<T>&
             ) const;
 
+            //- Print list according to (collected and) sorted coordinate
+            template<class T>
+            static void collectAndPrint
+            (
+                const pointField& points,
+                const UList<T>& data
+            );
+
             //- Print some mesh stats.
             void printMeshInfo(const bool, const string&) const;
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
index 812f313abe402b2b8e080d7a2b6c985b0f8fe559..13a96eba7982d8ec5bbdf3af2836a998b6bc6db9 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -25,6 +25,7 @@ License
 
 #include "meshRefinement.H"
 #include "fvMesh.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -112,11 +113,11 @@ T meshRefinement::gAverage
 (
     const polyMesh& mesh,
     const PackedBoolList& isMasterElem,
-    const labelList& meshPoints,
+    const labelList& meshElems,
     const UList<T>& values
 )
 {
-    if (values.size() != meshPoints.size())
+    if (values.size() != meshElems.size())
     {
         FatalErrorIn
         (
@@ -128,8 +129,8 @@ T meshRefinement::gAverage
             "    const UList<T>& values\n"
             ")\n"
         )   << "Number of elements in list " << values.size()
-            << " does not correspond to number of elements in meshPoints "
-            << meshPoints.size()
+            << " does not correspond to number of elements in meshElems "
+            << meshElems.size()
             << exit(FatalError);
     }
 
@@ -138,7 +139,7 @@ T meshRefinement::gAverage
 
     forAll(values, i)
     {
-        if (isMasterElem[meshPoints[i]])
+        if (isMasterElem[meshElems[i]])
         {
             sum += values[i];
             n++;
@@ -218,6 +219,53 @@ void meshRefinement::testSyncBoundaryFaceList
 }
 
 
+// Print list sorted by coordinates. Used for comparing non-parallel v.s.
+// parallel operation
+template<class T>
+void meshRefinement::collectAndPrint
+(
+    const pointField& points,
+    const UList<T>& data
+)
+{
+    globalIndex globalPoints(points.size());
+
+    pointField allPoints;
+    globalPoints.gather
+    (
+        Pstream::worldComm,
+        identity(Pstream::nProcs()),
+        points,
+        allPoints,
+        UPstream::msgType(),
+        Pstream::blocking
+    );
+
+    Field<T> allData;
+    globalPoints.gather
+    (
+        Pstream::worldComm,
+        identity(Pstream::nProcs()),
+        data,
+        allData,
+        UPstream::msgType(),
+        Pstream::blocking
+    );
+
+
+    scalarField magAllPoints(mag(allPoints));
+
+    labelList visitOrder;
+    sortedOrder(magAllPoints, visitOrder);
+    forAll(visitOrder, i)
+    {
+        label allPointI = visitOrder[i];
+        Info<< allPoints[allPointI] << " : " << allData[allPointI]
+            << endl;
+    }
+}
+
+
 //template<class T, class Mesh>
 template<class GeoField>
 void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)