From 4df782615ce0186250e76f4f7390938a5021f478 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 13 Mar 2018 17:58:02 +0100
Subject: [PATCH 1/5] BUG: surface proxy zone handling with dangling? reference
 (closes #757)

- using const reference to temporary was failing.  Remedy by using a
  direct copy, which is a reasonable solution since surfZone content
  is quite minimal.
---
 src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C     | 4 ++--
 src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C     | 2 +-
 src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C       | 2 +-
 src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C       | 2 +-
 src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C       | 2 +-
 src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C   | 2 +-
 src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C | 2 +-
 src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C       | 4 ++--
 src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C       | 2 +-
 src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C       | 2 +-
 src/surfMesh/surfaceFormats/vtp/VTPsurfaceFormat.C       | 2 +-
 src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C       | 2 +-
 12 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C b/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C
index 299e9a5333..e86bcc7962 100644
--- a/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C
@@ -312,7 +312,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
     const pointField& pointLst = surf.points();
     const UList<Face>& faceLst = surf.surfFaces();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().size()
       ? surf.surfZones()
@@ -393,7 +393,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
 
     if (zoneLst.size() <= 1)
     {
-        const List<surfZone>& zones =
+        const surfZoneList zones =
         (
             zoneLst.size()
           ? zoneLst
diff --git a/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C b/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C
index b834b4ae5f..7225a47cc4 100644
--- a/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C
@@ -169,7 +169,7 @@ void Foam::fileFormats::FLMAsurfaceFormat<Face>::write
     const UList<label>& faceMap = surf.faceMap();
 
     // for no zones, suppress the group name
-    const List<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst, word::null)
diff --git a/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C b/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C
index 352babafbb..ca460aa019 100644
--- a/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C
@@ -253,7 +253,7 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
     const UList<point>& pointLst = surf.points();
     const UList<Face>& faceLst  = surf.surfFaces();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().size()
       ? surf.surfZones()
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
index 03e66e3c75..503616989c 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
@@ -420,7 +420,7 @@ void Foam::fileFormats::NASsurfaceFormat<Face>::write
     const UList<label>& faceMap  = surf.faceMap();
 
     // for no zones, suppress the group name
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst, "")
diff --git a/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C b/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
index 77128da7ca..c2897bbac8 100644
--- a/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
@@ -217,7 +217,7 @@ void Foam::fileFormats::OBJsurfaceFormat<Face>::write
     const UList<label>& faceMap  = surf.faceMap();
 
     // for no zones, suppress the group name
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst, "")
diff --git a/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C b/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C
index 5f54451022..04552eb713 100644
--- a/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C
@@ -41,7 +41,7 @@ void Foam::fileFormats::SMESHsurfaceFormat<Face>::write
     const UList<Face>&  faceLst  = surf.surfFaces();
     const UList<label>& faceMap  = surf.faceMap();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst)
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
index b9f3f92afa..cc61fa0f29 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
@@ -241,7 +241,7 @@ void Foam::fileFormats::STARCDsurfaceFormat<Face>::write
     const UList<Face>&  faceLst  = surf.surfFaces();
     const UList<label>& faceMap  = surf.faceMap();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst)
diff --git a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
index 64ad423df8..4026f85e83 100644
--- a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
@@ -215,7 +215,7 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeAscii
     const UList<Face>&   faceLst = surf.surfFaces();
     const UList<label>&  faceMap = surf.faceMap();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst)
@@ -269,7 +269,7 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeBinary
     const UList<Face>&   faceLst = surf.surfFaces();
     const UList<label>&  faceMap = surf.faceMap();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().size() > 1
       ? surf.surfZones()
diff --git a/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C b/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C
index 898adcb833..a3b2046d27 100644
--- a/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C
@@ -161,7 +161,7 @@ void Foam::fileFormats::TRIsurfaceFormat<Face>::write
     const UList<Face>&   faceLst = surf.surfFaces();
     const UList<label>&  faceMap = surf.faceMap();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst)
diff --git a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
index f75920ae9a..1bb3d3cdaf 100644
--- a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
@@ -265,7 +265,7 @@ void Foam::fileFormats::VTKsurfaceFormat<Face>::write
     const UList<Face>&   faceLst = surf.surfFaces();
     const UList<label>&  faceMap = surf.faceMap();
 
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst)
diff --git a/src/surfMesh/surfaceFormats/vtp/VTPsurfaceFormat.C b/src/surfMesh/surfaceFormats/vtp/VTPsurfaceFormat.C
index bf3d09911a..c1f22611bd 100644
--- a/src/surfMesh/surfaceFormats/vtp/VTPsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/vtp/VTPsurfaceFormat.C
@@ -118,7 +118,7 @@ void Foam::fileFormats::VTPsurfaceFormat<Face>::write
     const UList<Face>&   faceLst = surf.surfFaces();
     const UList<label>&  faceMap = surf.faceMap();
 
-    const List<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst)
diff --git a/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C b/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C
index fc02e1a085..03aee4006a 100644
--- a/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C
@@ -41,7 +41,7 @@ void Foam::fileFormats::X3DsurfaceFormat<Face>::write
     const UList<label>&  faceMap = surf.faceMap();
 
     // for no zones, suppress the group name
-    const UList<surfZone>& zones =
+    const surfZoneList zones =
     (
         surf.surfZones().empty()
       ? surfaceFormatsCore::oneZone(faceLst, word::null)

From fae1ba0f6377780ad5f146b5088907dab1996430 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Wed, 14 Mar 2018 15:34:53 +0000
Subject: [PATCH 2/5] BUG: Corrections to externalWallHeatFluxTemperature BC
 (ref EP624)

---
 ...nalWallHeatFluxTemperatureFvPatchScalarField.C | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
index dbd1291172..99ab2d9e55 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C
@@ -58,6 +58,8 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     temperatureCoupledBase(patch(), "undefined", "undefined", "undefined-K"),
     mode_(fixedHeatFlux),
     Q_(0),
+    q_(),
+    h_(),
     Ta_(),
     relaxation_(1),
     emissivity_(0),
@@ -84,7 +86,9 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     temperatureCoupledBase(patch(), dict),
     mode_(operationModeNames.lookup("mode", dict)),
     Q_(0),
-    Ta_(Function1<scalar>::New("Ta", dict)),
+    q_(),
+    h_(),
+    Ta_(),
     relaxation_(dict.lookupOrDefault<scalar>("relaxation", 1)),
     emissivity_(dict.lookupOrDefault<scalar>("emissivity", 0)),
     qrRelaxation_(dict.lookupOrDefault<scalar>("qrRelaxation", 1)),
@@ -109,6 +113,7 @@ externalWallHeatFluxTemperatureFvPatchScalarField
         case fixedHeatTransferCoeff:
         {
             h_ = scalarField("h", dict, p.size());
+            Ta_ = Function1<scalar>::New("Ta", dict);
 
             if (dict.found("thicknessLayers"))
             {
@@ -175,12 +180,12 @@ externalWallHeatFluxTemperatureFvPatchScalarField
     temperatureCoupledBase(patch(), ptf),
     mode_(ptf.mode_),
     Q_(ptf.Q_),
-    q_(ptf.q_, mapper),
-    h_(ptf.h_, mapper),
+    q_(),
+    h_(),
     Ta_(ptf.Ta_, false),
     relaxation_(ptf.relaxation_),
     emissivity_(ptf.emissivity_),
-    qrPrevious_(ptf.qrPrevious_, mapper),
+    qrPrevious_(),
     qrRelaxation_(ptf.qrRelaxation_),
     qrName_(ptf.qrName_),
     thicknessLayers_(ptf.thicknessLayers_),
@@ -194,12 +199,14 @@ externalWallHeatFluxTemperatureFvPatchScalarField
         }
         case fixedHeatFlux:
         {
+            q_.setSize(mapper.size());
             q_.map(ptf.q_, mapper);
 
             break;
         }
         case fixedHeatTransferCoeff:
         {
+            h_.setSize(mapper.size());
             h_.map(ptf.h_, mapper);
 
             break;

From d7b9cea87c4326d6e82659fb0f94584d528b173b Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Mon, 19 Mar 2018 13:56:32 +0000
Subject: [PATCH 3/5] BUG: IsoAdvector - updated use of pos/neg as reported by
 Johan Roenby.  See #775

---
 .../isoAdvection/isoAdvection/isoAdvection.C     | 16 ++++++++--------
 .../solvers/isoAdvection/isoCutFace/isoCutFace.C |  8 ++++----
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
index e039b0b572..d149e0c5ef 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
@@ -615,8 +615,8 @@ void Foam::isoAdvection::limitFluxes()
     const scalar aTol = 1.0e-12;          // Note: tolerances
     const scalar maxAlphaMinus1 = 1;      // max(alphaNew - 1);
     const scalar minAlpha = -1;           // min(alphaNew);
-    const label nUndershoots = 20;        // sum(neg(alphaNew + aTol));
-    const label nOvershoots = 20;         // sum(pos(alphaNew - 1 - aTol));
+    const label nUndershoots = 20;        // sum(neg0(alphaNew + aTol));
+    const label nOvershoots = 20;         // sum(pos0(alphaNew - 1 - aTol));
     cellIsBounded_ = false;
 
     // Loop number of bounding steps
@@ -682,8 +682,8 @@ void Foam::isoAdvection::limitFluxes()
             scalarField alphaNew(alpha1In_ - fvc::surfaceIntegrate(dVf_)());
             label maxAlphaMinus1 = max(alphaNew - 1);
             scalar minAlpha = min(alphaNew);
-            label nUndershoots = sum(neg(alphaNew + aTol));
-            label nOvershoots = sum(pos(alphaNew - 1 - aTol));
+            label nUndershoots = sum(neg0(alphaNew + aTol));
+            label nOvershoots = sum(pos0(alphaNew - 1 - aTol));
             Info<< "After bounding number " << n + 1 << " of time "
                 << mesh_.time().value() << ":" << endl;
             Info<< "nOvershoots = " << nOvershoots << " with max(alphaNew-1) = "
@@ -792,7 +792,7 @@ void Foam::isoAdvection::boundFromAbove
                         fluidToPassOn*mag(phi[fi]*dt)/dVftot;
 
                     nFacesToPassFluidThrough +=
-                        pos(dVfmax[fi] - fluidToPassThroughFace);
+                        pos0(dVfmax[fi] - fluidToPassThroughFace);
 
                     fluidToPassThroughFace =
                         min(fluidToPassThroughFace, dVfmax[fi]);
@@ -1014,9 +1014,9 @@ void Foam::isoAdvection::applyBruteForceBounding()
     {
         alpha1_ =
             alpha1_
-           *pos(alpha1_ - snapAlphaTol)
-           *neg(alpha1_ - (1.0 - snapAlphaTol))
-          + pos(alpha1_ - (1.0 - snapAlphaTol));
+           *pos0(alpha1_ - snapAlphaTol)
+           *neg0(alpha1_ - (1.0 - snapAlphaTol))
+          + pos0(alpha1_ - (1.0 - snapAlphaTol));
 
         alpha1Changed = true;
     }
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
index 33d4e7f412..83292012a0 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
@@ -362,7 +362,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
     {
         // If all face cuttings were in the past and cell is filling up (Un0>0)
         // then face must be full during whole time interval
-        tIntArea = magSf*dt*pos(Un0);
+        tIntArea = magSf*dt*pos0(Un0);
         return tIntArea;
     }
 
@@ -373,7 +373,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
         // If all cuttings are in the future but non of them within [0,dt] then
         // if cell is filling up (Un0 > 0) face must be empty during whole time
         // interval
-        tIntArea = magSf*dt*(1 - pos(Un0));
+        tIntArea = magSf*dt*(1 - pos0(Un0));
         return tIntArea;
     }
 
@@ -398,7 +398,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
         // If Un0 > 0 cell is filling up and it must initially be empty.
         // If Un0 < 0 cell must initially be full(y immersed in fluid A).
         time = firstTime;
-        initialArea = magSf*(1.0 - pos(Un0));
+        initialArea = magSf*(1.0 - pos0(Un0));
         tIntArea = initialArea*time;
         cutPoints(fPts, pTimes, time, FIIL);
     }
@@ -470,7 +470,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
     {
         // FIIL will leave the face at lastTime and face will be fully in fluid
         // A or fluid B in the time interval from lastTime to dt.
-        tIntArea += magSf*(dt - lastTime)*pos(Un0);
+        tIntArea += magSf*(dt - lastTime)*pos0(Un0);
     }
 
     return tIntArea;

From 546fdb6207da5baca1f423e02929987e7d7ac5ca Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 20 Mar 2018 10:17:33 +0000
Subject: [PATCH 4/5] BUG: corrected 'binned' distribution model

---
 .../distributionModels/binned/binned.C        | 22 +++++++------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/src/lagrangian/distributionModels/binned/binned.C b/src/lagrangian/distributionModels/binned/binned.C
index 57530fa649..a5cf7e1f2b 100644
--- a/src/lagrangian/distributionModels/binned/binned.C
+++ b/src/lagrangian/distributionModels/binned/binned.C
@@ -47,10 +47,6 @@ const char* Foam::distributionModels::binned::header =
 void Foam::distributionModels::binned::initialise()
 {
     const label nSample(xy_.size());
-    forAll(xy_, bini)
-    {
-        xy_[bini][1] /= scalar(nSample);
-    }
 
     // Convert values to integral values
     for (label bini = 1; bini < nSample; ++bini)
@@ -58,6 +54,13 @@ void Foam::distributionModels::binned::initialise()
         xy_[bini][1] += xy_[bini - 1][1];
     }
 
+    // Normalise
+    scalar sumProb = xy_.last()[1];
+    forAll(xy_, bini)
+    {
+        xy_[bini][1] /= sumProb;
+    }
+
     // Calculate the mean value
     label bini = 0;
     forAll(xy_, i)
@@ -188,16 +191,7 @@ Foam::scalar Foam::distributionModels::binned::sample() const
     {
         if (xy_[i][1] > y)
         {
-            scalar d1 = y - xy_[i][1];
-            scalar d2 = xy_[i+1][1] - y;
-            if (d1 < d2)
-            {
-                return xy_[i][0];
-            }
-            else
-            {
-                return xy_[i+1][0];
-            }
+            return xy_[i][0];
         }
     }
 

From 1edb2f7dd0dc085f0403202703d764711584a17a Mon Sep 17 00:00:00 2001
From: mattijs <mattis>
Date: Wed, 21 Mar 2018 21:54:40 +0100
Subject: [PATCH 5/5] BUG: mesh refinement crash with degenerate mesh
 distributions (closes #778)

- occurred when the initial mesh distribution was missing cells on
  some processors.
---
 .../algorithms/indexedOctree/indexedOctree.C  | 50 ++++++++++++-------
 1 file changed, 33 insertions(+), 17 deletions(-)

diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index ada7d06b5b..64346e1d25 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -2560,6 +2560,11 @@ Foam::labelBits Foam::indexedOctree<Type>::findNode
 template<class Type>
 Foam::label Foam::indexedOctree<Type>::findInside(const point& sample) const
 {
+    if (nodes_.empty())
+    {
+        return -1;
+    }
+
     labelBits index = findNode(0, sample);
 
     const node& nod = nodes_[getNode(index)];
@@ -2592,6 +2597,11 @@ const Foam::labelList& Foam::indexedOctree<Type>::findIndices
     const point& sample
 ) const
 {
+    if (nodes_.empty())
+    {
+        return emptyList<label>();
+    }
+
     labelBits index = findNode(0, sample);
 
     const node& nod = nodes_[getNode(index)];
@@ -2603,10 +2613,8 @@ const Foam::labelList& Foam::indexedOctree<Type>::findIndices
     {
         return contents_[getContent(contentIndex)];
     }
-    else
-    {
-        return emptyList<label>();
-    }
+
+    return emptyList<label>();
 }
 
 
@@ -2688,18 +2696,21 @@ void Foam::indexedOctree<Type>::findNear
     CompareOp& cop
 ) const
 {
-    findNear
-    (
-        nearDist,
-        true,
-        *this,
-        nodePlusOctant(0, 0),
-        bb(),
-        tree2,
-        nodePlusOctant(0, 0),
-        tree2.bb(),
-        cop
-    );
+    if (!nodes_.empty())
+    {
+        findNear
+        (
+            nearDist,
+            true,
+            *this,
+            nodePlusOctant(0, 0),
+            bb(),
+            tree2,
+            nodePlusOctant(0, 0),
+            tree2.bb(),
+            cop
+        );
+    }
 }
 
 
@@ -2711,6 +2722,11 @@ void Foam::indexedOctree<Type>::print
     const label nodeI
 ) const
 {
+    if (nodes_.empty())
+    {
+        return;
+    }
+
     const node& nod = nodes_[nodeI];
     const treeBoundBox& bb = nod.bb_;