ENH: PDRsetFields utility (#1216)

- the PDRsetFields utility processes a set of geometrical obstructions
  to determine the equivalent blockage effects.

  These fields are necessary inputs for PDRFoam calculations.

  After setting up the geometries, the -dry-run option can be used to
  generate a VTK file for diagnosis and post-processing purposes.

- this is an initial release, with improvements slated for the future.

NOTE
  - the field results may be less than fully reliable when run in
    single-precision. This howver does not represent a realistic
    restriction since the prepared fields target a combustion
    application which will invariably be double-precision.

applications/utilities/preProcessing/PDRsetFields/Make/files

@@ -0,0 +1,19 @@
+PDRsetFields.C
+
+PDRarrays.C
+PDRarraysAnalyse.C
+PDRarraysCalc.C
+PDRmeshArrays.C
+PDRparams.C
+PDRpatchDef.C
+PDRlegacyMeshSpec.C
+PDRutilsIntersect.C
+PDRutilsOverlap.C
+
+obstacles/PDRobstacle.C
+obstacles/PDRobstacleIO.C
+obstacles/PDRobstacleTypes.C
+obstacles/PDRobstacleLegacyIO.C
+obstacles/PDRobstacleLegacyRead.C
+
+EXE = $(FOAM_APPBIN)/PDRsetFields

applications/utilities/preProcessing/PDRsetFields/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -Iobstacles \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfileFormats \
+    -lsurfMesh \
+    -lmeshTools \
+    -lblockMesh

applications/utilities/preProcessing/PDRsetFields/PDRarrays.C

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRarrays.H"
+#include "PDRblock.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Smaller helper to resize matrix and assign to Zero
+template<class T>
+inline void resizeMatrix(SquareMatrix<T>& mat, const label n)
+{
+    mat.setSize(n);
+    mat = Zero;
+}
+
+
+// Smaller helper to resize i-j-k field and assign to uniform value,
+// normally Zero
+template<class T>
+inline void resizeField
+(
+    IjkField<T>& fld,
+    const labelVector& ijk,
+    const T& val = T(Zero)
+)
+{
+    fld.resize(ijk);
+    fld = val;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PDRarrays::PDRarrays()
+:
+    pdrBlock_(std::cref<PDRblock>(PDRblock::null()))
+{}
+
+
+Foam::PDRarrays::PDRarrays(const PDRblock& pdrBlock)
+:
+    PDRarrays()
+{
+    reset(pdrBlock);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PDRarrays::reset(const PDRblock& pdrBlock)
+{
+    pdrBlock_ = std::cref<PDRblock>(pdrBlock);
+
+    // Resize all the major arrays, which are grouped in the structure arrp
+    // All the relevant dimensions are in PDRblock
+
+    // Cell-based addressing
+    const labelVector cellDims = pdrBlock.sizes();
+
+    // Face or point-based addressing
+    const labelVector faceDims(cellDims + labelVector::one);
+
+    // Max addressing dimensions for 2D arrays, with some extra space
+    // These will be used for any combination of x,y,z,
+    // so need to be dimensioned to the maximum size in both directions
+    const label maxDim = cmptMax(pdrBlock.sizes()) + 2;
+
+    resizeField(v_block, cellDims);
+    resizeField(surf, cellDims);
+
+    resizeField(area_block_s, cellDims);
+    resizeField(area_block_r, cellDims);
+    resizeField(dirn_block, cellDims);
+
+    resizeField(face_block, faceDims);
+
+    resizeField(along_block, cellDims);
+
+    resizeField(betai_inv1, cellDims);
+
+    resizeField(obs_count, cellDims);
+    resizeField(sub_count, cellDims);
+    resizeField(grating_count, cellDims);
+
+    resizeField(drag_s, cellDims);
+    resizeField(drag_r, cellDims);
+
+    resizeField(obs_size, cellDims);
+
+    for (auto& list : overlap_1d)
+    {
+        list.resize(maxDim);
+        list = Zero;
+    }
+
+    resizeMatrix(aboverlap, maxDim);
+    resizeMatrix(abperim, maxDim);
+    resizeMatrix(a_lblock, maxDim);
+    resizeMatrix(b_lblock, maxDim);
+    resizeMatrix(ac_lblock, maxDim);
+    resizeMatrix(bc_lblock, maxDim);
+    resizeMatrix(c_count, maxDim);
+    resizeMatrix(c_drag, maxDim);
+
+    resizeField(face_patch, faceDims, labelVector::uniform(-1));
+    resizeField(hole_in_face, faceDims);
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRarrays.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PDRarrays
+
+Description
+    Work array definitions for PDR fields
+
+SourceFiles
+    PDRarrays.C
+    PDRarraysCalc.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRarrays_H
+#define PDRarrays_H
+
+#include "symmTensor.H"
+#include "symmTensor2D.H"
+#include "SquareMatrix.H"
+#include "IjkField.H"
+#include <functional>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class PDRblock;
+class PDRmeshArrays;
+class PDRobstacle;
+class PDRpatchDef;
+
+
+/*---------------------------------------------------------------------------*\
+                          Class PDRarrays Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRarrays
+{
+    //- Reference to PDRblock
+    std::reference_wrapper<const PDRblock> pdrBlock_;
+
+public:
+
+    // Data Members
+    // Entries used for analysis and when writing fields
+
+    //- Volume blockage
+    IjkField<scalar> v_block;
+
+    //- Surface area in cell
+    IjkField<scalar> surf;
+
+    //- Obstacle size in cell
+    IjkField<scalar> obs_size;
+
+    //- Summed area blockage (directional) from sharp obstacles
+    IjkField<vector> area_block_s;
+
+    //- Summed area blockage (directional) from round obstacles
+    IjkField<vector> area_block_r;
+
+    //- A total directional blockage in the cell
+    IjkField<Vector<bool>> dirn_block;
+
+    //- Face area blockage for face,
+    //- summed from cell centre-plane to cell centre-plane
+    IjkField<vector> face_block;
+
+    //- Longitudinal area blockage from obstacles that extend all the way
+    //- through the cell in a given direction.
+    IjkField<vector> along_block;
+
+    IjkField<vector> betai_inv1;
+
+    //- Number of obstacles in cell.
+    //  Can be non-integer if an obstacle does not pass all way through cell
+    IjkField<scalar> obs_count;
+
+    //- Number of obstacles parallel to specified direction
+    IjkField<vector> sub_count;
+
+    //- Addition to count to account for grating comprises many bars
+    //- (to get Lobs right)
+    IjkField<vector> grating_count;
+
+    //- Tensorial drag from sharp obstacles
+    IjkField<symmTensor> drag_s;
+
+    //- Directional drag from round obstacles
+    IjkField<vector> drag_r;
+
+
+    // Next arrays are for 2D calculations of intersection
+
+    // One-dimensional scratch areas for cell overlaps
+    Vector<List<scalar>> overlap_1d;
+
+    // In two dimensions, area of cell covered by circle
+    SquareMatrix<scalar> aboverlap;
+
+    // In two dimensions, length of perimeter of circle witthin cell
+    SquareMatrix<scalar> abperim;
+
+    // For offset cells, i.e. face blockage
+    SquareMatrix<scalar> a_lblock, b_lblock;
+
+    // For centred cells
+    SquareMatrix<scalar> ac_lblock, bc_lblock;
+
+    // The count in the cells
+    SquareMatrix<scalar> c_count;
+
+    //- Cell-centred drag
+    SquareMatrix<symmTensor2D> c_drag;
+
+    //- Face field for (directional) for patch Id
+    IjkField<labelVector> face_patch;
+
+    //- Face field for (directional) hole in face
+    IjkField<Vector<bool>> hole_in_face;
+
+
+    // Constructors
+
+        //- Construct null
+        PDRarrays();
+
+        //- Construct and reset
+        explicit PDRarrays(const PDRblock& pdrBlock);
+
+
+    //- Destructor
+    ~PDRarrays() = default;
+
+
+    // Member Functions
+
+        //- Reset PDRblock reference, resize and zero arrays
+        void reset(const PDRblock& pdrBlock);
+
+        //- Reference to PDRblock
+        const PDRblock& block() const
+        {
+            return pdrBlock_.get();
+        }
+
+        //- Summary of the blockages
+        //  For diagnostics and general overview
+        void blockageSummary() const;
+
+        //- Add cylinder blockage
+        void addCylinder(const PDRobstacle& obs);
+
+        //- Add general (non-cylinder) blockage
+        void addBlockage
+        (
+            const PDRobstacle& obs,
+            DynamicList<PDRpatchDef>& patches,
+            const int volumeSign
+        );
+
+
+        static void calculateAndWrite
+        (
+            PDRarrays& arr,
+            const PDRmeshArrays& meshIndexing,
+            const fileName& casepath,
+            const UList<PDRpatchDef>& patches
+        );
+
+        void calculateAndWrite
+        (
+            const fileName& casepath,
+            const PDRmeshArrays& meshIndexing,
+            const UList<PDRpatchDef>& patches
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRlegacy.H

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::PDRlegacy
+
+Description
+    Legacy and transitional routines
+
+SourceFiles
+    PDRlegacyMeshSpec.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRlegacy_H
+#define PDRlegacy_H
+
+#include "PDRblock.H"
+#include "fileName.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class ISstream;
+class PDRblock;
+
+/*---------------------------------------------------------------------------*\
+                              Namespace PDRlegacy
+\*---------------------------------------------------------------------------*/
+
+namespace PDRlegacy
+{
+    void print_info(const PDRblock& block);
+
+    void read_mesh_spec(const fileName& casepath, PDRblock& pdrBlock);
+
+    void read_mesh_spec(ISstream& is, PDRblock& pdrBlock);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRlegacyMeshSpec.C

@@ -0,0 +1,340 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// Read spec that resemble this:
+//
+//   xmesh
+//   (
+//       ( -8.18  14  0.83 )
+//       ( -0.69  11  1.00 )
+//       (  0.69  14  1.20 )
+//       (  8.18  0 )
+//   )
+//
+//   ymesh
+//   (
+//       ...
+//   )
+
+#include "PDRlegacy.H"
+
+// OpenFOAM includes
+#include "error.H"
+#include "IFstream.H"
+#include "stringOps.H"
+#include "OSspecific.H"
+
+#define XMESH_TAG "xmesh"
+#define YMESH_TAG "ymesh"
+#define ZMESH_TAG "zmesh"
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PDRlegacy
+{
+namespace Detail
+{
+
+//- Simple structure for processing a mesh spec entry
+//  Handles an entry with (float), (float int), or (float int float)
+//
+//  This is for handling a sub-spec that resembles this:
+//
+//  (
+//      ( -8.18  14  0.83 )
+//      ( -0.69  11  1.00 )
+//      (  0.69  14  1.20 )
+//      (  8.18  0 )
+//  )
+struct pdrMeshSpecLine
+{
+    scalar knot;
+    label  ndiv;
+    scalar factor;
+
+    pdrMeshSpecLine() : knot(0), ndiv(0), factor(0) {}
+
+    //- Cheap means to avoid uniform list output
+    bool operator!=(const pdrMeshSpecLine&) const
+    {
+        return false;
+    }
+};
+
+
+// Read mesh-spec entry
+Istream& operator>>(Istream& is, pdrMeshSpecLine& spec)
+{
+    spec.knot = 0;
+    spec.ndiv = 0;
+    spec.factor = 0;
+
+    is.readBegin("pdrMeshSpecLine");
+
+    // Must have a point
+    is >> spec.knot;
+
+    token tok(is);
+    if (tok.isLabel())
+    {
+        spec.ndiv = tok.labelToken();
+
+        if (spec.ndiv)
+        {
+            is >> spec.factor;
+        }
+    }
+    else
+    {
+        is.putBack(tok);
+    }
+
+    is.readEnd("pdrMeshSpecLine");
+
+    is.check(FUNCTION_NAME);
+    return is;
+}
+
+#ifdef FULLDEBUG
+// Write mesh-spec entry
+Ostream& operator<<(Ostream& os, const pdrMeshSpecLine& spec)
+{
+    os << token::BEGIN_LIST << spec.knot;
+
+    if (spec.ndiv)
+    {
+        os  << token::SPACE << spec.ndiv
+            << token::SPACE << spec.factor;
+    }
+
+    os << token::END_LIST;
+
+    return os;
+}
+#endif
+
+
+void read_spec(ISstream& is, const direction cmpt, List<scalar>& gridPoint)
+{
+    if (!gridPoint.empty())
+    {
+        FatalErrorInFunction
+            << "Duplicate specification of "
+            << vector::componentNames[cmpt]
+            << " grid"
+            << exit(FatalError);
+    }
+
+    List<pdrMeshSpecLine> specs(is);
+
+    if (specs.size() < 2)
+    {
+        FatalErrorInFunction
+            << "Grid specification for " << vector::componentNames[cmpt]
+            << " is too small. Need at least two points!" << nl
+            << exit(FatalError);
+    }
+
+    specs.last().ndiv = 0; // safety
+
+
+    DynamicList<scalar> knots;
+    DynamicList<label>  divisions;
+    DynamicList<scalar> factors;
+
+    for (const auto& spec : specs)
+    {
+        knots.append(spec.knot);
+
+        if (spec.ndiv < 1)
+        {
+            break;
+        }
+        divisions.append(spec.ndiv);
+        factors.append(spec.factor);
+    }
+
+    label nPoints = 1;
+    for (const label nDiv : divisions)
+    {
+        nPoints += nDiv;
+    }
+
+    if (nPoints < 2)
+    {
+        FatalErrorInFunction
+            << "No cells defined for direction "
+            << vector::componentNames[cmpt] << nl
+            << exit(FatalError);
+    }
+
+
+    // Define the grid points
+    gridPoint.resize(nPoints);
+
+    const label nSegments = divisions.size();
+
+    label start = 0;
+
+    for (label segmenti=0; segmenti < nSegments; ++segmenti)
+    {
+        const label nDiv = divisions[segmenti];
+        const scalar factor = factors[segmenti];
+
+        SubList<scalar> subPoint(gridPoint, nDiv+1, start);
+        start += nDiv;
+
+        subPoint[0] = knots[segmenti];
+        subPoint[nDiv] = knots[segmenti+1];
+
+        const scalar dist = (subPoint.last() - subPoint.first());
+
+        if (equal(factor, scalar(1)))
+        {
+            for (label i=1; i < nDiv; ++i)
+            {
+                subPoint[i] = (subPoint[0] + (dist * i)/nDiv);
+            }
+        }
+        else
+        {
+            scalar delta = dist * (1.0 - factor) / (1.0 - ::pow(factor, nDiv));
+
+            scalar xyz = subPoint[0];
+
+            for (label i=0; i < nDiv; ++i)
+            {
+                subPoint[i] = xyz;
+                xyz += delta;
+                delta *= factor;
+            }
+        }
+    }
+}
+
+
+} // End namespace Detail
+} // End namespace PDRlegacy
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+void Foam::PDRlegacy::print_info(const PDRblock& block)
+{
+    Info<< "PDRblock" << nl
+        << "    nCells: " << block.sizes() << nl
+        << "    Box: " << block.bounds() << nl
+        << "x " << flatOutput(block.grid().x()) << nl
+        << "y " << flatOutput(block.grid().y()) << nl
+        << "z " << flatOutput(block.grid().z()) << nl
+        << endl;
+}
+
+
+void Foam::PDRlegacy::read_mesh_spec(const fileName& casepath, PDRblock& block)
+{
+    Info<< "Reading pdrMeshSpec (legacy format)" << nl;
+
+    bool processed = false;
+
+    for (const fileName dirName : { "system", "constant/polyMesh" })
+    {
+        fileName path
+        (
+            casepath / dirName / "pdrMeshSpec"
+        );
+
+        if (Foam::isFile(path))
+        {
+            IFstream is(path);
+
+            read_mesh_spec(is, block);
+            processed = true;
+            break;
+        }
+    }
+
+    if (!processed)
+    {
+        FatalErrorInFunction
+            << "Did not process pdrMeshSpec" << nl
+            << exit(FatalError);
+    }
+}
+
+
+void Foam::PDRlegacy::read_mesh_spec(ISstream& is, PDRblock& block)
+{
+    Vector<scalarList> grid;
+
+    string line;
+
+    while (is.good())
+    {
+        is.getLine(line);
+        stringOps::inplaceTrim(line);
+
+        if (line == XMESH_TAG)
+        {
+            Detail::read_spec(is, vector::X, grid.x());
+        }
+        else if (line == YMESH_TAG)
+        {
+            Detail::read_spec(is, vector::Y, grid.y());
+        }
+        else if (line == ZMESH_TAG)
+        {
+            Detail::read_spec(is, vector::Z, grid.z());
+        }
+    }
+
+    for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+    {
+        if (grid[cmpt].empty())
+        {
+            FatalErrorInFunction
+                << "No specification for "
+                << vector::componentNames[cmpt]
+                << " grid" << nl
+                << exit(FatalError);
+        }
+    }
+
+    block.reset(grid.x(), grid.y(), grid.z());
+
+    #ifdef FULLDEBUG
+    print_info(block);
+    #endif
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRmeshArrays.C

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRmeshArrays.H"
+#include "PDRblock.H"
+#include "polyMesh.H"
+#include "Time.H"
+#include "IjkField.H"
+
+// Notes
+//
+// Determines the face and cell numbers of all faces and cells in the
+// central rectangular region where CAD_PDR operates. First,
+// "points" is read and the coordinates (by which I mean here the
+// indices in the x, y and z coordinate arrays) are determined. Then
+// "faces" is read and for each the coordinates of the lower- x,y,z
+// corner are determioned, also the orientation (X, Y or Z).
+// (Orientation in the sense of e.g. + or -x is not noted.) The files
+// "owner" and "neighbour" specify the six faces around each cell, so
+// from these the coordinates of the cells are determined.
+//
+// Full checks are made that the mesh in the central region is consistent
+// with CAD_PDR's mesh specified by the PDRmeshSpec file.
+//
+// Eventually, when writing out results, we shall work through the
+// full list of cells, writing default values for any cells that are
+// not in the central regtion.
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::scalar Foam::PDRmeshArrays::gridPointRelTol = 0.02;
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PDRmeshArrays::classify
+(
+    const polyMesh& mesh,
+    const PDRblock& pdrBlock
+)
+{
+    // Additional copy of i-j-k addressing
+    cellDims = pdrBlock.sizes();
+    faceDims = (cellDims + labelVector::one);
+
+    const label maxPointId = cmptMax(pdrBlock.sizes())+1;
+
+    Info<< "Mesh" << nl
+        << "  nPoints:" << mesh.nPoints()
+        << "  nCells:" << mesh.nCells()
+        << "  nFaces:" << mesh.nFaces() << nl;
+
+    Info<< "PDRblock" << nl
+        << "  minEdgeLen:" << pdrBlock.minEdgeLen() << nl;
+
+
+    // Bin points into i-j-k locations
+    List<labelVector> pointIndex(mesh.nPoints());
+
+    for (label pointi=0; pointi < mesh.nPoints(); ++pointi)
+    {
+        const point& pt = mesh.points()[pointi];
+        pointIndex[pointi] = pdrBlock.gridIndex(pt, gridPointRelTol);
+    }
+
+    // Min x,y,z index
+    const labelMinMax invertedLimits(maxPointId, -maxPointId);
+    Vector<labelMinMax> faceLimits;
+
+    const Vector<direction> faceBits
+    (
+        boundBox::XDIR,
+        boundBox::YDIR,
+        boundBox::ZDIR
+    );
+
+    faceIndex.resize(mesh.nFaces());
+    faceOrient.resize(mesh.nFaces());
+
+    for (label facei=0; facei < mesh.nFaces(); ++facei)
+    {
+        faceLimits.x() = faceLimits.y() = faceLimits.z() = invertedLimits;
+
+        for (const label pointi : mesh.faces()[facei])
+        {
+            for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
+            {
+                faceLimits[cmpt].add(pointIndex[pointi][cmpt]);
+            }
+        }
+
+        direction inPlane(0u);
+
+        for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
+        {
+            const auto& limits = faceLimits[cmpt];
+
+            if (!limits.valid())
+            {
+                // This should be impossible
+                FatalErrorInFunction
+                    << "Unexpected search failure for " << facei << " in "
+                    << vector::componentNames[cmpt] << "-direction" << nl
+                    << exit(FatalError);
+            }
+
+            if (limits.min() < 0)
+            {
+                FatalErrorInFunction
+                    << "Face " << facei << " contains non-grid point in "
+                    << vector::componentNames[cmpt] << "-direction" << nl
+                    << exit(FatalError);
+            }
+            else if (limits.min() == limits.max())
+            {
+                // In plane
+                inPlane |= faceBits[cmpt];
+            }
+            else if (limits.min() + 1 != limits.max())
+            {
+                FatalErrorInFunction
+                    << "Face " << facei
+                    << " not in " << vector::componentNames[cmpt] << "-plane" << nl
+                    << exit(FatalError);
+            }
+        }
+
+        switch (inPlane)
+        {
+            case boundBox::XDIR:
+                faceOrient[facei] = vector::X;
+                break;
+
+            case boundBox::YDIR:
+                faceOrient[facei] = vector::Y;
+                break;
+
+            case boundBox::ZDIR:
+                faceOrient[facei] = vector::Z;
+                break;
+
+            default:
+                FatalErrorInFunction
+                    << "Face " << facei << " not in an x/y/z plane?" << nl
+                    << exit(FatalError);
+                break;
+        }
+
+        faceIndex[facei] =
+            labelVector
+            (
+                faceLimits.x().min(),
+                faceLimits.y().min(),
+                faceLimits.z().min()
+            );
+    }
+
+
+    // Bin cells into i-j-k locations
+    cellIndex = std::move(pointIndex);
+    cellIndex = labelVector::uniform(maxPointId);
+    cellIndex.resize(mesh.nCells(), labelVector::uniform(maxPointId));
+
+    // Option 1: use PDRblock.findCell() method
+    if (true)
+    {
+        const pointField& cc = mesh.cellCentres();
+
+        for (label celli=0; celli < mesh.nCells(); ++celli)
+        {
+            cellIndex[celli] = pdrBlock.findCell(cc[celli]);
+        }
+    }
+
+    // Option 2: walk cell faces and use faceIndex information
+    if (false)
+    {
+        for (label celli=0; celli < mesh.nCells(); ++celli)
+        {
+            labelVector& cellIdx = cellIndex[celli];
+
+            for (const label facei : mesh.cells()[celli])
+            {
+                cellIdx.x() = min(cellIdx.x(), faceIndex[facei].x());
+                cellIdx.y() = min(cellIdx.y(), faceIndex[facei].y());
+                cellIdx.z() = min(cellIdx.z(), faceIndex[facei].z());
+            }
+
+            if (cmptMin(cellIdx) < 0)
+            {
+                cellIdx = labelVector(-1,-1,-1);
+            }
+        }
+    }
+
+
+    // Verify that all i-j-k cells were found
+    {
+        // This could be more efficient - but we want to be picky
+        IjkField<bool> cellFound(pdrBlock.sizes(), false);
+
+        for (label celli=0; celli < cellIndex.size(); ++celli)
+        {
+            const labelVector& cellIdx = cellIndex[celli];
+
+            if (cmptMin(cellIdx) >= 0)
+            {
+                cellFound(cellIdx) = true;
+            }
+        }
+
+        label firstMissing = cellFound.find(false);
+
+        if (firstMissing >= 0)
+        {
+            FatalErrorInFunction
+                << "No cell found for " << pdrBlock.index(firstMissing)
+                << " indexing"
+                << exit(FatalError);
+        }
+    }
+}
+
+
+void Foam::PDRmeshArrays::read
+(
+    const Time& runTime,
+    const PDRblock& pdrBlock
+)
+{
+    #include "createPolyMesh.H"
+    classify(mesh, pdrBlock);
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRmeshArrays.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PDRmeshArrays
+
+Description
+    OpenFOAM/PDRblock addressing information
+
+    Provides mapping for a rectilinear OpenFOAM mesh in terms of
+    i-j-k indices for faces and cells.
+
+    The mesh points are first binned according to their i-j-k locations.
+    Next the faces are classified according to their lowest x/y/z
+    coordinates and the face orientation as x/y/z.
+    Orientation in the sense +x or -x is not noted.
+    The cell faces are then examined to determine the appropriate i-j-k
+    location.
+
+SourceFiles
+    PDRmeshmeshArraysIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRmeshArrays_H
+#define PDRmeshArrays_H
+
+#include "labelVector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class PDRblock;
+class polyMesh;
+class Time;
+
+/*---------------------------------------------------------------------------*\
+                          Class PDRmeshArrays Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRmeshArrays
+{
+public:
+
+    //- Relative tolerance when matching grid points. Default = 0.02
+    static scalar gridPointRelTol;
+
+    //- The cell i-j-k addressing range
+    labelVector cellDims;
+
+    //- The face i-j-k addressing range
+    labelVector faceDims;
+
+    //- For each cell, the corresponding i-j-k address.
+    List<labelVector> cellIndex;
+
+    //- For each face, the corresponding i-j-k address.
+    List<labelVector> faceIndex;
+
+    //- For each face, the x/y/z orientation
+    List<direction> faceOrient;
+
+
+    // Constructors
+
+        //- Construct null
+        PDRmeshArrays() = default;
+
+
+    //- Destructor
+    ~PDRmeshArrays() = default;
+
+
+    // Member Functions
+
+        //- The number of cells
+        label nCells() const
+        {
+            return cellIndex.size();
+        }
+
+        //- The number of faces
+        label nFaces() const
+        {
+            return faceIndex.size();
+        }
+
+
+        //- Determine i-j-k indices for faces/cells
+        void classify(const polyMesh& mesh, const PDRblock& pdrBlock);
+
+        //- Read OpenFOAM mesh and determine i-j-k indices for faces/cells
+        void read(const Time& runTime, const PDRblock& pdrBlock);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRparams.C

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRparams.H"
+#include "stringOps.H"
+
+// * * * * * * * * * * * * * * * * Global Data * * * * * * * * * * * * * * * //
+
+// Global parameter settings
+Foam::PDRparams Foam::pars;
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PDRparams::readDefaults(const dictionary& dict)
+{
+    dict.readIfPresent("legacyMeshSpec", legacyMeshSpec);
+    dict.readIfPresent("legacyObsSpec",  legacyObsSpec);
+
+    dict.readIfPresent("two_d", two_d);
+    dict.readIfPresent("yCyclic", yCyclic);
+    dict.readIfPresent("ySymmetry", ySymmetry);
+    dict.readIfPresent("deluge", deluge);
+
+    dict.readIfPresent("newFields", new_fields);
+    dict.readIfPresent("noIntersectN", noIntersectN);
+
+    dict.readIfPresent("blockedFacesWallFn", blockedFacesWallFn);
+    dict.readIfPresent("ignoreGratings", ignoreGratings);
+
+    outer_orthog = dict.found("outer_orthog");
+
+    dict.readIfPresent("debugLevel", debugLevel);
+    dict.readIfPresent("nFacesToBlockC", nFacesToBlockC);
+    dict.readIfPresent("nPairsToBlockC", nPairsToBlockC);
+    dict.readIfPresent("overlaps", overlaps);
+
+    dict.readIfPresent("gridPointTol", gridPointTol);
+
+    dict.readIfPresent("Cb_r", cb_r);
+    dict.readIfPresent("Cb_s", cb_s);
+
+    dict.readIfPresent("Cd_r", cd_r);
+    dict.readIfPresent("Cd_s", cd_s);
+
+    dict.readIfPresent("congRegionMaxBetav", cong_max_betav);
+
+    dict.readIfPresent("min_overlap_vol", min_overlap_vol);
+    dict.readIfPresent("min_overlap_area", min_overlap_area);
+    dict.readIfPresent("min_width", min_width);
+    dict.readIfPresent("empty_lobs_fac", empty_lobs_fac);
+    dict.readIfPresent("outerCombFac", outerCombFac);
+    dict.readIfPresent("obs_expand", obs_expand);
+
+    dict.readIfPresent("def_grating_slat_w", def_grating_slat_w);
+    dict.readIfPresent("blockedCellPoros", blockedCellPoros);
+    dict.readIfPresent("blockedFacePar", blockedFacePar);
+    dict.readIfPresent("maxCR", maxCR);
+
+    dict.readIfPresent("blockageNoCT", blockageNoCT);
+    dict.readIfPresent("scale", scale);
+
+    UPatchBc = "fixedValue;value uniform (0 0 0)";
+    if (dict.readIfPresent("UPatchBc", UPatchBc))
+    {
+        stringOps::inplaceTrim(UPatchBc);
+    }
+}
+
+
+void Foam::PDRparams::read(const dictionary& dict)
+{
+    readDefaults(dict);
+
+    dict.readEntry("obsFileDir", obsfile_dir);
+    dict.readEntry("obsFileNames", obsfile_names);
+
+    stringOps::inplaceExpand(obsfile_dir);
+
+    for (auto& f : obsfile_names)
+    {
+        stringOps::inplaceExpand(f);
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRparams.H

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2018-2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PDRparams
+
+Description
+    Parameters for PDRsetFields
+
+SourceFiles
+    PDRparams.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRparams_H
+#define PDRparams_H
+
+#include "labelList.H"
+#include "scalarList.H"
+#include "wordList.H"
+#include "fileNameList.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class PDRparams Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRparams
+{
+public:
+
+    // Data Members
+
+        fileName obsfile_dir;
+        wordList obsfile_names;
+        word    timeName;
+        string  UPatchBc;   //!< "fixedValue;value uniform (0 0 0)"
+
+        bool legacyMeshSpec{false};
+        bool legacyObsSpec{false};
+
+        bool two_d{false};
+        bool yCyclic{false};
+        bool ySymmetry{false};
+        bool deluge{false};
+        bool new_fields{true};
+        bool noIntersectN{true};
+        bool blockedFacesWallFn{false};
+        bool ignoreGratings{false};
+        bool outer_orthog{false};
+
+        int debugLevel{0};
+
+        //- Min number of blocked cell faces
+        //- for a cell to be marked as blocked
+        int nFacesToBlockC{6};
+
+        //- Min number of blocked cell face pairs (on opposite faces of a cell)
+        //- for a cell to be marked as blocked
+        int nPairsToBlockC{3};
+
+        //- Flag to control which overlap calculations are performed
+        int overlaps{0x7};
+
+        scalar gridPointTol{0.02};
+
+        scalar cb_r{0.035};
+        scalar cb_s{0.08};
+
+        scalar cd_r{1.2};
+        scalar cd_s{2.0};
+
+        scalar cong_max_betav{1.0};
+
+        scalar min_overlap_vol{0};
+        scalar min_overlap_area{0};
+
+        //- Ignore obstacles with second dimension (or diameter) less than this
+        scalar min_width{0.001};
+
+        //- Lobs in empty cell is this * cube root of cell volume
+        scalar empty_lobs_fac{1.0};
+
+        //- Value for outer region
+        scalar outerCombFac{1.0};
+
+        scalar obs_expand{0};
+
+        //- Default slat thickness grating
+        scalar def_grating_slat_w{0.005};
+
+        //- Cells with porosity less than this are blocked
+        scalar blockedCellPoros{0.05};
+
+        //- Faces with area blockage greater than this are blocked
+        scalar blockedFacePar{0.95};
+
+        //- Upper limit on CR (CT also gets limited)
+        scalar maxCR{1e30};
+
+        //- If a single obstacle blocks a cell by more than this,
+        //- then no CT in that direction
+        scalar blockageNoCT{0.95};
+
+        //- Overall scale factor
+        scalar scale{1.0};
+
+
+    // Constructors
+
+        //- Construct null
+        PDRparams() = default;
+
+
+    // Member Functions
+
+        //- Set or read defaults from dictionary.
+        //  Can also be used with an empty dictionary
+        void readDefaults(const dictionary& dict);
+
+        //- Read program parameters from dictionary
+        void read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+extern Foam::PDRparams pars;  //!< Globals for program parameters (ugly hack)
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRpatchDef.C

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRpatchDef.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::Enum
+<
+    Foam::PDRpatchDef::predefined
+>
+Foam::PDRpatchDef::names
+({
+    { predefined::BLOCKED_FACE, "blockedFaces" },
+    { predefined::MERGING_PATCH, "mergingFaces" },
+    { predefined::WALL_PATCH, "wallFaces" },
+});
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::PDRpatchDef::operator=(const std::string& newName)
+{
+    patchName = word::validate(newName);
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRpatchDef.H

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PDRpatchDef
+
+Description
+    Bookkeeping for patch definitions
+
+SourceFiles
+    PDRpatchDef.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRpatchDef_H
+#define PDRpatchDef_H
+
+#include "string.H"
+#include "scalar.H"
+#include "Enum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class PDRpatchDef Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRpatchDef
+{
+public:
+
+    //- Patch predefines
+    enum predefined
+    {
+        BLOCKED_FACE = 0,
+        MERGING_PATCH = 1,
+        WALL_PATCH = 2,
+        LAST_PREDEFINED = 2, // First user patch number will be 1 more
+        NUM_PREDEFINED = 3
+    };
+
+    //- Names for predefined types
+    static const Enum<predefined> names;
+
+
+    // Data Members
+
+    word patchName;
+
+    label patchType;
+
+    scalar blowoffPress;
+
+    scalar blowoffTime;
+
+
+    // Constructors
+
+        //- Construct null
+        PDRpatchDef()
+        :
+             patchName(),
+             patchType(0),
+             blowoffPress(0),
+             blowoffTime(0)
+        {}
+
+        //- Construct with given patch name
+        explicit PDRpatchDef(const word& name)
+        :
+             patchName(name),
+             patchType(0),
+             blowoffPress(0),
+             blowoffTime(0)
+        {}
+
+
+        //- Construct with given patch name
+        PDRpatchDef& operator=(const PDRpatchDef&) = default;
+        PDRpatchDef& operator=(PDRpatchDef&&) = default;
+
+        //- Assign new patch name
+        void operator=(const std::string& newName);
+};
+
+
+typedef PDRpatchDef PATCH;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRsetFields.C

@@ -0,0 +1,351 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Applications
+    PDRsetFields
+
+Description
+    Preparation of fields for PDRFoam
+
+SourceFiles
+    PDRsetFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "IOdictionary.H"
+
+#include "PDRsetFields.H"
+#include "PDRlegacy.H"
+#include "PDRutils.H"
+#include "IOmanip.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char* argv[])
+{
+    argList::addNote
+    (
+        "Processes a set of geometrical obstructions to determine the"
+        " equivalent blockage effects when setting cases for PDRFoam"
+    );
+    argList::noParallel();
+    argList::noFunctionObjects();
+
+    argList::addOption
+    (
+        "time",
+        "time",
+        "Specify a time"
+    );
+
+    argList::addOption("dict", "file", "Use alternative PDRsetFieldsDict");
+
+    argList::addBoolOption
+    (
+        "legacy",
+        "Force use of legacy obstacles table"
+    );
+
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Read obstacles and write VTK only"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    const word dictName("PDRsetFieldsDict");
+    #include "setSystemRunTimeDictionaryIO.H"
+
+    Info<< "Reading " << dictName << "\n" << endl;
+
+    IOdictionary setFieldsDict(dictIO);
+
+    const bool dryrun = args.found("dry-run");
+
+    const fileName& casepath = runTime.globalPath();
+
+    pars.timeName = "0";
+    args.readIfPresent("time", pars.timeName);
+
+    // Program parameters (globals)
+    pars.read(setFieldsDict);
+
+    if (args.found("legacy"))
+    {
+        pars.legacyObsSpec = true;
+    }
+
+    // Always have the following:
+    // 0 = blockedFaces patch (no wall functions)
+    // 1 = mergingFaces patch
+    // 2 = wallFaces patch
+
+    DynamicList<PDRpatchDef> patches;
+    patches.resize(PDRpatchDef::NUM_PREDEFINED);
+
+    for
+    (
+        PDRpatchDef::predefined predef :
+        {
+            PDRpatchDef::BLOCKED_FACE,
+            PDRpatchDef::MERGING_PATCH,
+            PDRpatchDef::WALL_PATCH,
+        }
+    )
+    {
+        patches[predef] = PDRpatchDef::names[predef];
+    }
+
+
+    // Dimensions and grid points for i-j-k domain
+    PDRblock pdrBlock;
+
+    if (pars.legacyMeshSpec)
+    {
+        PDRlegacy::read_mesh_spec(casepath, pdrBlock);
+    }
+    else
+    {
+        IOdictionary iodict
+        (
+            IOobject
+            (
+                "PDRblockMeshDict",
+                runTime.system(),
+                runTime,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            )
+        );
+
+        pdrBlock.read(iodict);
+
+        #ifdef FULLDEBUG
+        PDRlegacy::print_info(pdrBlock);
+        #endif
+    }
+
+    // Storage for obstacles and cylinder-like obstacles
+    DynamicList<PDRobstacle> obstacles, cylinders;
+
+    // Read in obstacles
+    const scalar volObstacles =
+    (
+        pars.legacyObsSpec
+      ? PDRobstacle::legacyReadFiles
+        (
+            pars.obsfile_dir, pars.obsfile_names,
+            pdrBlock.bounds(),
+            obstacles,
+            cylinders
+        )
+      : PDRobstacle::readFiles
+        (
+            pars.obsfile_dir, pars.obsfile_names,
+            pdrBlock.bounds(),
+            obstacles,
+            cylinders
+        )
+    );
+
+
+    PDRobstacle::generateVtk(casepath/"VTK", obstacles, cylinders);
+
+    if (dryrun)
+    {
+        Info<< nl
+            << "dry-run: stopping after reading/writing obstacles" << nl
+            << "\nEnd\n" << nl;
+        return 0;
+    }
+
+
+    // Bookkeeping of the ranges within the obstacle list
+
+    // Original blockage at the start
+    const labelRange origBlocks(0, obstacles.size());
+
+    // Intersection blockage
+    labelRange interBlocks(origBlocks.after(), 0);
+
+    scalar volSubtract = 0;
+
+    // Do binary intersections between blocks and cylinders (or diag-beam)
+    // by creating -ve blocks at the overlap
+
+    labelRange int1Blocks(origBlocks.after(), 0);
+
+    if (pars.overlaps % 2 > 0)
+    {
+        Info<< "    block/cylinder intersections" << endl;
+
+        label nblocked = obstacles.size();
+
+        volSubtract += block_cylinder_overlap(obstacles, origBlocks, cylinders);
+
+        nblocked = (obstacles.size() - nblocked);
+
+        interBlocks += nblocked;
+        int1Blocks += nblocked;
+    }
+
+    // Do binary intersections between blocks
+    // by creating -ve blocks at the overlap
+
+    labelRange int2Blocks(int1Blocks.after(), 0);
+    if (pars.overlaps % 4 > 1)
+    {
+        Info<< "    block/block intersections" << endl;
+
+        label nblocked = obstacles.size();
+
+        volSubtract += block_overlap(obstacles, origBlocks, 1.0);
+
+        nblocked = (obstacles.size() - nblocked);
+
+        interBlocks += nblocked;
+        int2Blocks += nblocked;
+    }
+
+    // Correct for triple intersections
+    // by looking for overlaps between the -ve blocks just created
+
+    labelRange int3Blocks(int2Blocks.after(), 0);
+    if (pars.overlaps % 8 > 3)
+    {
+        Info<< "    triple intersections" << endl;
+
+        label nblocked = obstacles.size();
+
+        volSubtract += block_overlap(obstacles, interBlocks, 1.0/3.0);
+
+        nblocked = (obstacles.size() - nblocked);
+
+        interBlocks += nblocked;
+        int3Blocks += nblocked;
+    }
+
+
+    // The field arrays, in one structure pass around easily
+    PDRarrays arr(pdrBlock);
+
+    Info<< "Apply blockage" << endl;
+
+    // Create blockage and count arrays by working through
+    // real and extra blocks and cylinders
+
+    // User-defined negative blocks. Use "sign" to distinguish
+    if (origBlocks.size())
+    {
+        Info<< "    negative blocks: " << origBlocks.size() << nl;
+
+        for (const PDRobstacle& obs : obstacles[origBlocks])
+        {
+            arr.addBlockage(obs, patches, -1);
+        }
+    }
+
+    // Do the intersection blocks positive and negative
+    // These are done first so that negative area blockage cancels positive
+
+    if (interBlocks.size())
+    {
+        Info<< "    blocks " << interBlocks.size() << nl;
+
+        for (const PDRobstacle& obs : obstacles[interBlocks])
+        {
+            arr.addBlockage(obs, patches, 0);
+        }
+    }
+
+    // The positive real bocks
+    if (origBlocks.size())
+    {
+        Info<< "    positive blocks: " << origBlocks.size() << nl;
+
+        for (const PDRobstacle& obs : obstacles[origBlocks])
+        {
+            arr.addBlockage(obs, patches, 1);
+        }
+    }
+
+    // The cylinders
+    if (cylinders.size())
+    {
+        Info<< "    cylinders: " << cylinders.size() << nl;
+
+        for (const PDRobstacle& obs : cylinders)
+        {
+            arr.addCylinder(obs);
+        }
+    }
+
+    // Calculation of the fields of drag, turbulence
+    // generation and combustion enhancement
+
+    arr.blockageSummary();
+
+    // Mapping of OpenFOAM cells/faces to i-j-k indices
+    PDRmeshArrays meshIdx;
+    meshIdx.gridPointRelTol = pars.gridPointTol;
+
+    meshIdx.read(runTime, pdrBlock);
+
+    PDRarrays::calculateAndWrite(arr, meshIdx, casepath, patches);
+
+    Info<< nl
+        << setw(6) << origBlocks.size() << " blocks and "
+        << cylinders.size() << " cylinders/diagonal blocks" << nl;
+
+    Info<< setw(6) << int2Blocks.size()
+        << " intersections amongst blocks" << nl;
+
+    Info<< setw(6) << int1Blocks.size()
+        << " intersections between blocks and cyl/beams" << nl;
+
+    Info<< setw(6) << int1Blocks.size()
+        << "/3 triple intersections" << nl;
+
+    Info<< "Volume of obstacles read in: " << volObstacles
+        << ", volume of intersections: " << volSubtract << nl;
+
+    Info<< nl << "After corrections:" << nl;
+    arr.blockageSummary();
+
+    Info<< nl << "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRsetFields.H

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Preparation of fields for PDRFoam
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRsetFields_H
+#define PDRsetFields_H
+
+#include "PDRarrays.H"
+#include "PDRblock.H"
+#include "PDRmeshArrays.H"
+#include "PDRobstacle.H"
+#include "PDRpatchDef.H"
+#include "PDRparams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam;
+
+
+//YCYCLIC is set to 1 for running the test cases with a cyclic boundry condition in the y direction
+//TWO_D is set to 1 for running the 2D test cases (no z direction) - usually same case as YCYCLIC
+//Now specified in CAD_PDRDict and read in as globals.
+
+// The program also labels intersection obstacles as types 96 and 86, but not then use these values
+// Obstacles in group definitions have a multiple of 100 added to the type number
+
+#define floatSMALL  1.e-10
+
+// Default initial values of field variables, used outside congested area,
+//and everywhere for uniform fields. They atre strings because same routines
+//are used to create b.c.s for scalars and tensors.
+#define DEFAULT_K       0.00015
+#define DEFAULT_EPS     1e-5
+#define DEFAULT_T       300
+#define DEFAULT_P       100000
+#define DEFAULT_SU      0.5
+#define DEFAULT_LOBS    0.1   // Does not matter what it is outside congestion
+                              // but zero would cause problems with Ep
+#define DEFAULT_EP      0.01  // Gives length scale 0.1, calc. as (Xp-0.999)/Ep with Xp=1
+
+// Boundary conditions on walls for all variables where it is not "zero_gradient"
+#define K_WALL_FN       "kqRWallFunction"
+#define EPS_WALL_FN     "epsilonWallFunction"
+#define ALPHAT_WALL     "nutkWallFunction"
+#define MUT_WALL_FN     "mutkWallFunction"
+#define NUT_WALL_FN     "nutkWallFunction"
+
+#define K_WALL_FN_LEGACY        "compressible::kqRWallFunction"
+#define EPS_WALL_FN_LEGACY      "compressible::epsilonWallFunction"
+#define ALPHAT_WALL_FN_LEGACY   "alphatWallFunction;\n\t\tPrt\t0.85"
+
+
+// The following parameters are used to decide when there arMAX_Ne sufficient (parts of)
+// obstacles ina cell for them to define the length scale of the generated turbulence.
+#define MIN_AB_FOR_SIZE 0.002
+#define MAX_VB_FOR_SIZE 0.9
+#define COUNT_FOR_SIZE  0.1
+#define MIN_COUNT_FOR_SIZE 0.05
+
+// These define how blocked a face or cell has to be for removal from the mesh
+//#define BLOCKED_CELL_PAR  0.05  //<-  Now pars.blockedCellPoros
+//#define BLOCKED_FACE_PAR  0.95  //<-  Now pars.blockedFacePar
+
+
+//- Calculate block/block overlaps
+//
+// Binary self-intersections are to be checked for blocks.
+// Resulting negative blocks are appended to blocks.
+// These new blocks have the opposite sign from input blocks, and
+// blockage multiplied by multiplier.
+//
+// If the number of newly generated blocks is required, check the size
+// of blocks on output vs input to see how many have been added.
+//
+// \param[in,out] blocks
+// \param[in] range - the range within blocks to be examined
+//
+// \return overlap volume
+scalar block_overlap
+(
+    DynamicList<PDRobstacle>& blocks,
+    const labelRange& range,
+    const scalar multiplier = 1.0
+);
+
+
+//- Calculate block/cylinder overlaps
+//
+// Binary intersections are to be checked for blocks and cylinders.
+// Resulting negative blocks are appended to blocks.
+// These new blocks have the opposite sign from input blocks, and
+// blockage multiplied by multiplier.
+//
+// If the number of newly generated blocks is required, check the size
+// of blocks on output vs input to see how many have been added.
+//
+// \param[in,out] arrp
+// \param[in,out] blocks
+// \param[in] range - the range within blocks to be examined
+// \param[in] cylinders - the cylinders to be examined
+//
+//  \return overlap volume
+scalar block_cylinder_overlap
+(
+    DynamicList<PDRobstacle>& blocks,
+    const labelRange& range,
+    const UList<PDRobstacle>& cylinders
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRutils.H

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::PDRutils
+
+Description
+    Utilities for PDR (eg, for setFields)
+
+SourceFiles
+    PDRUtils.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRutils_H
+#define PDRutils_H
+
+#include "PDRarrays.H"
+#include "PDRblock.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class PDRobstacle;
+
+namespace PDRutils
+{
+
+
+} // End namespace PDRutils
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRutilsInternal.H

@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::PDRutils
+
+Description
+    Utilities for PDR (eg, for setFields). Internal usage only.
+
+    The C lineage of the original code is still evident in the use of
+    pointers instead of references.
+    This will be addressed in later versions of the code (2019-12).
+
+SourceFiles
+    PDRUtils.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRutilsInternal_H
+#define PDRutilsInternal_H
+
+#include "PDRutils.H"
+#include "PDRarrays.H"
+#include "PDRblock.H"
+#include "symmTensor2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PDRutils
+{
+
+//- Determine 1-D overlap locations for a geometric entity
+//
+// \param[in] xmin - min position of the geometric entity
+// \param[in] xmax - max position of the geometric entity
+// \param[in] grid - grid point information
+// \param[out] olap - Fraction of cell-width with overlap
+//     0 for no overlap, 1 for full overlap.
+// \param[out] cmin - first cell index (inclusive) with overlap,
+//     values in the range \c [0,nCells]
+// \param[out] cmax - last cell index (inclusive) with overlap,
+//     values in the range \c [0,nCells]
+// \param[out] cfmin - first cell index (inclusive) with staggered face,
+//     values in the range \c [0,nCells]
+// \param[out] cfmax - last cell index (inclusive) with staggered face,
+//     values in the range \c [0,nCells]
+void one_d_overlap
+(
+    scalar xmin,
+    scalar xmax,
+    const PDRblock::location& grid,
+    List<scalar>& olap,
+    int *cmin, int *cmax,
+    int *cfmin, int *cfmax
+);
+
+
+//- Combine two 1D overlaps.
+//  Multiplying the two 1-d overlaps yields the proportion of each (2D) cell
+//  that is covered.
+//
+//  \note We go one over the relevant min/max limits since these
+//  values might be used.
+//  The 1D arrays will have bee initially zeroed throughout.
+void two_d_overlap
+(
+    const UList<scalar>& a_olap, label amin, label amax,
+    const UList<scalar>& b_olap, label bmin, label bmax,
+    SquareMatrix<scalar>& ab_olap
+);
+
+
+//- Calculate the proportion of each (two-dimensional) grid cell
+//- overlapped by the circle or angled rectangle.
+//
+//  Coordinates are labelled a and b.
+//
+// \param[in] ac, bc  coordinates of centre of circle or rectangle
+// \param[in] dia      diameter of circle (zero for rectangle)
+// \param[in] theta, wa, wb     parameters for rectangle
+// \param[in] amin, amax first and last cells in a-grid overlapped by object
+// \param[in] agrid               locations of grid lines of a-grid
+// \param[in] amin, amax first and last cells in b-grid overlapped by object
+// \param[in] bgrid               locations of grid lines of b-grid
+//
+// \param[out] abolap
+//     2-D array of (proportionate) area blockage by grid cell
+// \param[out] a_lblock
+//     2-D array of (proportionate) blockage to a-direction flow
+//     (This will be area blockage when extruded in the third
+//     coordinate).
+//
+// \param[out] a_count
+//   2-D array The contribution of this object to the count of
+// obstacles blocking a-direction flow. This is only non-zero if the
+// object is inside the lateral boundaries of the cell. It is large
+// negative if the cell is totally blocked in this direction.
+//
+//
+// \param[out] c_drag
+//
+// 2-D array of tensor that will give tensor drag in each cell (when
+// multiplied Cd, cylinder length, and 0.5 rho*U^2) Dimension: L.
+//
+// \note this routine does not zero array elements outside the amin
+//  to amax, bmin to bmax area.
+void circle_overlap
+(
+    scalar ac, scalar bc, scalar dia,
+    scalar theta, scalar wa, scalar wb,
+    const PDRblock::location& agrid, label amin, label amax,
+    const PDRblock::location& bgrid, label bmin, label bmax,
+    SquareMatrix<scalar>& ab_olap,
+    SquareMatrix<scalar>& ab_perim,
+    SquareMatrix<scalar>& a_lblock,
+    SquareMatrix<scalar>& ac_lblock,
+    SquareMatrix<scalar>& c_count,
+    SquareMatrix<symmTensor2D>& c_drag,
+    SquareMatrix<scalar>& b_lblock,
+    SquareMatrix<scalar>& bc_lblock
+);
+
+
+//- Area of intersection between circle and rectangle.
+//
+// Calculates the area of intersection between the circle, centre (xc, yc), radius rad,
+// and the rectangle with sides at x = x1 & x2, and y = y1 and y2.
+//
+// The return value is the fraction of the rectangle's area covered by the circle.
+double inters_cy
+(
+    double xc,   //!< circle centre (x)
+    double yc,   //!< circle centre (y)
+    double rad,  //!< circle radius
+    double x1, double x2,
+    double y1, double y2,
+    scalar* perim_p,
+    scalar* x_proj_edge_p,
+    scalar* y_proj_edge_p,
+    scalar* x_overlap_p,
+    scalar* y_overlap_p
+);
+
+
+//- The area overlap in the plane of a diagonal block and a cell.
+//
+// Calculates the overlap, in the plane of a diagonal block and a cell,
+// plus blockage and drag parameters.
+// Note that x and y herein may be any two of the three coordinates - would have been better not to label them x and y.
+//
+// On entry:
+//   xc, yc                       Coordinates of axis of d.b.
+//   theta, wa, wb        Angle and widths
+//
+//   The returned parameters will be multipled by the length of the obstacle's intersection with
+//   the third dimension of the 3-D cell to give this obstacle's contribution to the count, drag
+//   and area blockages.
+//   The return value is the area of intersection, which will multiply to volume blockage.
+//
+double inters_db
+(
+    double xc, double yc, double theta,
+    double wa, double wb,
+    double x1, double x2,
+    double y1, double y2,
+    scalar* count_p,
+    symmTensor2D& vdrag, scalar* perim_p,
+    scalar* x_lblk, scalar* y_lblk,
+    scalar* x_centre_p, scalar* y_centre_p
+);
+
+
+/* Calculates the blockage to x-direction flow presented by the specified circle on
+ the specified rectangle.
+ Becomes the area blockage when extruded to in the third dimension.
+ In other words, it is the projection on the y axis of the intersection between the
+ circle and the rectangle.
+ Returns fraction blocked
+ Note that x and y in this routine may in fact be any two of the three dimensions.
+ */
+double l_blockage
+(
+    double xc, double yc, double rad,
+    double x1, double x2,
+    double y1, double y2,
+    scalar* count_p, scalar* drag_p, scalar* centre_p
+);
+
+
+} // End namespace PDRutils
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRutilsIntersect.C

@@ -0,0 +1,712 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRsetFields.H"
+#include "PDRutilsInternal.H"
+#include "mathematicalConstants.H"
+
+#ifndef FULLDEBUG
+#define NDEBUG
+#endif
+#include <cassert>
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Calculate the area of the sector of a circle whose ends are at
+// (dxa, dya) and (dxb, dyb) relative to the centre. radsqu is radius
+// squared.
+//
+// (We trust that this is consistent with the other parameters..)
+inline static scalar sector
+(
+    scalar dxa, scalar dya,
+    scalar dxb, scalar dyb
+)
+{
+    scalar angle = (::atan2(dyb, dxb) - ::atan2(dya, dxa));
+
+    if (angle < -1.0E-10)
+    {
+        angle += mathematical::twoPi;
+    }
+
+    return angle;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+double Foam::PDRutils::inters_cy
+(
+    double xc, double yc, double rad,
+    double x1, double x2,
+    double y1, double y2,
+    scalar* perim_p,
+    scalar* x_proj_edge_p, scalar* y_proj_edge_p,
+    scalar* x_overlap_p, scalar* y_overlap_p
+)
+{
+    double angle, area, del;
+    double x_int[6][2], y_int[6][2]; // Coordinates of intersections between the circle and sides of the rectangle.
+    double x_arc[6][2], y_arc[6][2]; // Coordinates of end orc (within cell) for each quadrant
+    double dx[6], dy[6];
+
+    double x_olap_min = GREAT;
+    double x_olap_max = -GREAT;
+    double y_olap_min = GREAT;
+    double y_olap_max = -GREAT;
+    int n_vert, n_oppv;
+    int no_intersection;
+
+    const double dx1 = (x1 - xc);
+    const double dx2 = (x2 - xc);
+    const double dy1 = (y1 - yc);
+    const double dy2 = (y2 - yc);
+
+    const double dx1squ = dx1 * dx1;
+    const double dx2squ = dx2 * dx2;
+    const double dy1squ = dy1 * dy1;
+    const double dy2squ = dy2 * dy2;
+    const double radsqu = rad * rad;
+
+
+    /* Going clockwise from (x1, y1), vertices are labelled 1,2,3,4, with 0 the same as 4
+     and 5 the same as 1 (so that we can find the vertices on either side of any of 1 to 4).*/
+
+    dx[1] = dx1; dy[1] = dy1;
+    dx[2] = dx1; dy[2] = dy2;
+    dx[3] = dx2; dy[3] = dy2;
+    dx[4] = dx2; dy[4] = dy1;
+    dx[0] = dx2; dy[0] = dy1;
+    dx[5] = dx1; dy[5] = dy1;
+
+    // The positions of the ends of the arcs, if these points are
+    // inside the cell, they will be changed, if necessary, below.
+
+    x_arc[2][0] = x_arc[1][1] = -rad; y_arc[2][0] = y_arc[1][1] =  0.0;
+    x_arc[3][0] = x_arc[2][1] =  0.0; y_arc[3][0] = y_arc[2][1] =  rad;
+    x_arc[4][0] = x_arc[3][1] =  rad; y_arc[4][0] = y_arc[3][1] =  0.0;
+    x_arc[1][0] = x_arc[4][1] =  0.0; y_arc[1][0] = y_arc[4][1] = -rad;
+
+    // We catch arcs that are entirely inside the rectangle
+    // Note: this is wrong for a circle completely outside, but that
+    // will be dealt with separately
+
+    int arc_in[6] = { /* zero-initialied */ };
+    arc_in[1] = (dx1 < -rad && dy1 < -rad) ? 1 : 0;
+    arc_in[2] = (dx1 < -rad && dy2 >  rad) ? 1 : 0;
+    arc_in[3] = (dx2 >  rad && dy2 >  rad) ? 1 : 0;
+    arc_in[4] = (dx2 >  rad && dy1 < -rad) ? 1 : 0;
+
+    // Work out which vertices are in the circle
+
+    int vert_in[6];
+    vert_in[1] = (dx1squ + dy1squ <= radsqu);
+    vert_in[2] = (dx1squ + dy2squ <= radsqu);
+    vert_in[3] = (dx2squ + dy2squ <= radsqu);
+    vert_in[4] = (dx2squ + dy1squ <= radsqu);
+    vert_in[0] = vert_in[4];
+    vert_in[5] = vert_in[1];
+
+    int n_in = 0;
+    if (vert_in[1]) ++n_in;
+    if (vert_in[2]) ++n_in;
+    if (vert_in[3]) ++n_in;
+    if (vert_in[4]) ++n_in;
+
+
+    /* We now calculate the points of intersection of the circle with, successively,
+     x=x1, y=y2, x=x2. y=y1.
+
+     Where there are two intersections with one side, need to be careful about
+     the order of the two points (i.e. clockwise round the rectangle) so that
+     later on we get the right sector (short or long way round the circumference) */
+
+    int n_int[6] = { /* zero-initialied */ };
+    n_int[1] = 0;
+    if ( dx1squ <= radsqu)
+    {
+        del = std::sqrt( radsqu - dx1squ);
+        if ( ( ( -del ) <= dy2 ) && ( del >= dy1 ) )
+        {
+            x_int[1][0] = x_int[1][1] = dx1;
+            if ( (-del ) > dy1 )
+            {
+                y_int[1][0] = -del;
+                n_int[1]++;
+                // This intersection will be an end of the 3rd- or 4th-quadrant arc
+                if ( dx1 > 0.0 ) { x_arc[4][1] = dx1; y_arc[4][1] = -del; arc_in[4] = 1; }
+                else             { x_arc[1][1] = dx1; y_arc[1][1] = -del; arc_in[1] = 1; }
+            }
+            if ( ( del ) < dy2 )
+            {
+                y_int[1][n_int[1]] =  del;
+                n_int[1]++;
+                if ( dx1 > 0.0 ) { x_arc[3][0] = dx1; y_arc[3][0] = del; arc_in[3] = 1; }
+                else             { x_arc[2][0] = dx1; y_arc[2][0] = del; arc_in[2] = 1; }
+            }
+        }
+    }
+
+    n_int[2] = 0;
+    if ( dy2squ <= radsqu)
+    {
+        del = std::sqrt( radsqu - dy2squ);
+        if ( ( ( -del ) <= dx2 ) && ( del >= dx1 ) )
+        {
+            y_int[2][0] = y_int[2][1] = dy2;
+            if ( (-del ) > dx1 )
+            {
+                x_int[2][0] = -del;
+                n_int[2]++;
+                if ( dy2 > 0.0 ) { x_arc[2][1] = -del; y_arc[2][1] = dy2; arc_in[2] = 1; }
+                else             { x_arc[1][1] = -del; y_arc[1][1] = dy2; arc_in[1] = 1; }
+            }
+            if ( ( del ) < dx2 )
+            {
+                x_int[2][n_int[2]] =  del;
+                n_int[2]++;
+                if ( dy2 > 0.0 ) { x_arc[3][0] = del; y_arc[3][0] = dy2; arc_in[3] = 1; }
+                else             { x_arc[4][0] = del; y_arc[4][0] = dy2; arc_in[4] = 1; }
+            }
+        }
+    }
+
+    n_int[3] = 0;
+    if ( dx2squ <= radsqu)
+    {
+        del = std::sqrt( radsqu - dx2squ);
+        if ( ( ( -del ) <= dy2 ) && ( del >= dy1 ) )
+        {
+            x_int[3][0] = x_int[3][1] = dx2;
+            if ( ( del ) < dy2 )
+            {
+                y_int[3][0] = del;
+                n_int[3]++;
+                if ( dx2 > 0.0 ) { x_arc[3][1] = dx2; y_arc[3][1] = del; arc_in[3] = 1; }
+                else             { x_arc[2][1] = dx2; y_arc[2][1] = del; arc_in[2] = 1; }
+            }
+            if ( (-del ) > dy1 )
+            {
+                y_int[3][n_int[3]] = -del;
+                n_int[3]++;
+                if ( dx2 > 0.0 ) { x_arc[4][0] = dx2; y_arc[4][0] = -del; arc_in[4] = 1; }
+                else             { x_arc[1][0] = dx2; y_arc[1][0] = -del; arc_in[1] = 1; }
+            }
+        }
+    }
+
+    n_int[4] = 0;
+    if ( dy1squ <= radsqu)
+    {
+        del = std::sqrt( radsqu - dy1squ);
+        if ( ( ( -del ) <= dx2 ) && ( del >= dx1 ) )
+        {
+            y_int[4][0] = y_int[4][1] = dy1;
+            if ( ( del ) < dx2 )
+            {
+                x_int[4][0] = del;
+                n_int[4]++;
+                if ( dy1 > 0.0 ) { x_arc[3][1] = del; y_arc[3][1] = dy1; arc_in[3] = 1; }
+                else             { x_arc[4][1] = del; y_arc[4][1] = dy1; arc_in[4] = 1; }
+            }
+            if ( (-del ) > dx1 )
+            {
+                x_int[4][n_int[4]] = -del;
+                n_int[4]++;
+                if ( dy1 > 0.0 ) { x_arc[2][0] = -del; y_arc[2][0] = dy1; arc_in[2] = 1; }
+                else             { x_arc[1][0] = -del; y_arc[1][0] = dy1; arc_in[1] = 1; }
+            }
+        }
+    }
+
+    n_int[0] = n_int[4];
+    n_int[5] = n_int[1];
+
+    y_int[0][0] = y_int[0][1] = dy1;
+    x_int[0][0] = x_int[4][0];
+    x_int[0][1] = x_int[4][1];
+    x_int[5][0] = x_int[5][1] = dx1;
+    y_int[5][0] = y_int[1][0];
+    y_int[5][1] = y_int[1][1];
+
+    /* There are five separate cases, depending of the number of vertices inside the circle */
+    switch ( n_in )
+    {
+        case 0:
+        {
+            /* We start with the whole area of the circle, and then subtract any bits that stick out. */
+            area =  radsqu * mathematical::pi;
+            *perim_p = mathematical::twoPi * rad;
+            no_intersection = true;
+            for (n_vert = 1; n_vert < 5; n_vert++)
+            {
+                assert(n_int[n_vert] != 1);
+                if (n_int[n_vert] == 2)
+                {
+                    /* The area of the bit to be subtracted is a sector minus a triangle. */
+                    no_intersection = false;
+                    angle = sector( x_int[n_vert][1], y_int[n_vert][1], x_int[n_vert][0], y_int[n_vert][0]);
+                    area -= angle * radsqu * 0.5;
+                    *perim_p -= angle * rad;
+                    /* Two trinagles specified here, but one has zero area. */
+                    area += ( - ( y_int[n_vert][1] - y_int[n_vert][0] ) * x_int[n_vert][0]
+                              + ( x_int[n_vert][1] - x_int[n_vert][0] ) * y_int[n_vert][0] ) / 2.0;
+                }
+            }
+            /* Need to allow for when the circle is completely out side the rectanglle
+             by checking if the centre is outside the rectangle                                           */
+            if ( no_intersection )
+            {
+                if ( (dx1>0) ||(dx2<0) || (dy1>0) || (dy2<0) )
+                {
+                    *perim_p = *x_proj_edge_p = *y_proj_edge_p = 0.0;
+                    area = *x_overlap_p = *y_overlap_p = 0.0;
+                    return area;
+                }
+            }
+
+            break;
+        }
+
+        case 1:
+        {
+            /* Find which vertex is inside */
+            n_vert = 1;
+            while ( !vert_in[n_vert] ) { n_vert++; assert( n_vert < 5 ); }
+            assert( n_int[n_vert-1] == 1 );
+            if ( n_int[n_vert] != 1 )
+            {
+                assert( n_int[n_vert] == 1 );
+            }
+            angle = sector( x_int[n_vert-1][0], y_int[n_vert-1][0], x_int[n_vert][0], y_int[n_vert][0]);
+            area = angle * radsqu * 0.5;
+            *perim_p = angle * rad;
+            /* We subtract (or add) two triangles; the other two evaluate to zero */
+            area -= ( - ( x_int[n_vert][0]   - dx[n_vert] ) * dy[n_vert]
+                      + ( x_int[n_vert-1][0] - dx[n_vert] ) * dy[n_vert]
+                      + ( y_int[n_vert][0]   - dy[n_vert] ) * dx[n_vert]
+                      - ( y_int[n_vert-1][0] - dy[n_vert] ) * dx[n_vert] ) / 2.0;
+
+            break;
+        }
+
+        case 2:
+        {
+            /* This time n_vert is the number of the side which is completely inside the circle */
+            n_vert = 1;
+            while ( !(vert_in[n_vert] && vert_in[n_vert+1]) ) { n_vert++; assert( n_vert < 5 ); }
+            assert( n_int[n_vert-1] == 1 );
+            assert( n_int[n_vert+1] == 1 );
+            angle = sector( x_int[n_vert-1][0], y_int[n_vert-1][0], x_int[n_vert+1][0], y_int[n_vert+1][0]);
+            area = angle * radsqu * 0.5;
+            *perim_p = angle * rad;
+            /* We subtract (or add) three triangles; the other three evaluate to zero */
+            area += ( ( x_int[n_vert+1][0] - dx[n_vert+1] ) * dy[n_vert+1]
+                      - ( x_int[n_vert-1][0] - dx[n_vert] ) * dy[n_vert]
+                      - ( y_int[n_vert+1][0] - dy[n_vert+1] ) * dx[n_vert+1]
+                      + ( y_int[n_vert-1][0] - dy[n_vert] ) * dx[n_vert]
+                      + ( dx[n_vert+1] -dx[n_vert] ) * dy[n_vert]
+                      - ( dy[n_vert+1] -dy[n_vert] ) * dx[n_vert] ) / 2.0;
+
+            switch ( n_vert )
+            {
+                case 1: x_olap_min = dx1; break;
+                case 2: y_olap_max = dy2; break;
+                case 3: x_olap_max = dx2; break;
+                case 4: y_olap_min = dy1; break;
+            }
+
+            break;
+        }
+
+        case 3:
+        {
+            /* Find which vertex is NOT inside */
+            n_vert = 1;
+            while ( vert_in[n_vert] ) { n_vert++; assert( n_vert < 5 ); }
+            assert( n_int[n_vert-1] == 1 );
+            assert( n_int[n_vert] == 1 );
+            n_oppv = (n_vert + 2) % 4;
+            angle = sector( x_int[n_vert][0], y_int[n_vert][0], x_int[n_vert-1][0], y_int[n_vert-1][0]);
+            area = angle * radsqu * 0.5;
+            *perim_p = angle * rad;
+            /* We subtract (or add) four triangles; the other four evaluate to zero */
+            area += ( - ( x_int[n_vert][0]   - dx[n_vert+1] ) * dy[n_vert+1]
+                      + ( x_int[n_vert-1][0] - dx[n_vert-1] ) * dy[n_vert-1]
+                      + ( y_int[n_vert][0]   - dy[n_vert+1] ) * dx[n_vert+1]
+                      - ( y_int[n_vert-1][0] - dy[n_vert-1] ) * dx[n_vert-1]
+                      + ( dx[n_oppv] -dx[n_vert+1] ) * dy[n_oppv]
+                      - ( dx[n_oppv] -dx[n_vert-1] ) * dy[n_oppv]
+                      - ( dy[n_oppv] -dy[n_vert+1] ) * dx[n_oppv]
+                      + ( dy[n_oppv] -dy[n_vert-1] ) * dx[n_oppv] ) / 2.0;
+
+            x_olap_min = dx1;
+            y_olap_max = dy2;
+            x_olap_max = dx2;
+            y_olap_min = dy1;
+
+            break;
+        }
+
+        case 4:
+        {
+            /* Easy! We have the whole rectangle.  */
+            area = *x_overlap_p = *y_overlap_p = 1.0; // Normalised
+            *perim_p = *x_proj_edge_p = *y_proj_edge_p = 0.0;
+            return area;
+
+            break;
+        }
+    }
+
+    // The area may be very small negative by rounding errors
+    assert(area >=-1.0E-4);
+    if (area < 0.0) area = 0.0;
+    /* Return the overlap as a fraction of the rectangle's area. */
+    area /= ( (x2 - x1 ) * ( y2 - y1 ) );
+
+    // Sum the parts of the circumference that are inside the circle, projected onto axes
+    *x_proj_edge_p =
+    (
+        (y_arc[1][1] - y_arc[1][0]) * arc_in[1]
+      + (y_arc[2][1] - y_arc[2][0]) * arc_in[2]
+      + (y_arc[3][0] - y_arc[3][1]) * arc_in[3]
+      + (y_arc[4][0] - y_arc[4][1]) * arc_in[4]
+    );
+
+    *y_proj_edge_p =
+    (
+        (x_arc[1][0] - x_arc[1][1]) * arc_in[1]
+      + (x_arc[2][1] - x_arc[2][0]) * arc_in[2]
+      + (x_arc[3][1] - x_arc[3][0]) * arc_in[3]
+      + (x_arc[4][0] - x_arc[4][1]) * arc_in[4]
+    );
+
+    if (arc_in[1])
+    {
+        x_olap_min = min(x_olap_min, x_arc[1][1]);
+        x_olap_max = max(x_olap_max, x_arc[1][0]);
+        y_olap_min = min(y_olap_min, y_arc[1][0]);
+        y_olap_max = max(y_olap_max, y_arc[1][1]);
+    }
+    if (arc_in[2])
+    {
+        x_olap_min = min(x_olap_min, x_arc[2][0]);
+        x_olap_max = max(x_olap_max, x_arc[2][1]);
+        y_olap_min = min(y_olap_min, y_arc[2][0]);
+        y_olap_max = max(y_olap_max, y_arc[2][1]);
+    }
+    if (arc_in[3])
+    {
+        x_olap_min = min(x_olap_min, x_arc[3][0]);
+        x_olap_max = max(x_olap_max, x_arc[3][1]);
+        y_olap_min = min(y_olap_min, y_arc[3][1]);
+        y_olap_max = max(y_olap_max, y_arc[3][0]);
+    }
+    if (arc_in[4])
+    {
+        x_olap_min = min(x_olap_min, x_arc[4][1]);
+        x_olap_max = max(x_olap_max, x_arc[4][0]);
+        y_olap_min = min(y_olap_min, y_arc[4][1]);
+        y_olap_max = max(y_olap_max, y_arc[4][0]);
+    }
+
+    *x_overlap_p = ( x_olap_max - x_olap_min ) / ( x2 - x1 );
+    *y_overlap_p = ( y_olap_max - y_olap_min ) / ( y2 - y1 );
+    assert ( *x_overlap_p >= -floatSMALL );
+    assert ( *y_overlap_p >= -floatSMALL );
+
+    return area;
+} // End intersect
+
+
+// ************************************************************************* //
+
+double Foam::PDRutils::l_blockage
+(
+    double xc, double yc, double rad,
+    double x1, double x2,
+    double y1, double y2,
+    scalar* count_p, scalar* drag_p, scalar* centre_p
+)
+{
+    double xi = 0.0, lb, lb1, lb2, del;
+    bool within = true; // Indicates that the the intersection does not overlap the ends of the line
+
+    /* xi is the side we need to calc. intersections with */
+    if      ( xc < x1 ) { xi = x1; }
+    else if ( xc > x2 ) { xi = x2; }
+
+    if ( xi == 0.0 )
+    {
+        del = rad;     // The relevant lowest ( or highest) point is at end of vertical radius
+    }
+    else                // The relevant lowest ( or highest) point at intersection with x = xi
+    {
+        del = rad*rad - ( xi - xc ) * ( xi - xc );
+        if ( del < 0.0 ) { del = 0.0; }                  // No intersection
+        else             { del = std::sqrt(del); }
+    }
+
+    if ( ( yc + del ) > y2 ) { lb2  = y2; within = false; } else { lb2  = yc + del; }
+    if ( ( yc - del ) < y1 ) { lb1 = y1; within = false; } else { lb1 = yc - del; }
+
+    lb = (lb2 - lb1) / (y2 - y1);
+    *centre_p = (lb2 + lb1) * 0.5;
+
+    if ( lb < 0.0 )  lb = 0.0;
+
+    /* *count_p is 0 if the circle overlaps either y-side of the rectangle,
+     1 if the circle is entirely inside the rectangle
+     reduced if it overlaps x-sides.
+     A negative value indicates total blockage*/
+    if ( within && (lb > 0.0) )
+    {
+        *count_p = 1.0;
+        if ( ( xc - rad ) < x1 ) *count_p -= 0.5;
+        if ( ( xc + rad ) > x2 ) *count_p -= 0.5;
+    }
+    else
+    {
+        *count_p = 0.0;
+    }
+    *drag_p = lb * 1.2; //*drag_p = lb * CD_ROUND;
+    if ( lb > 0.99 ) { *count_p = -1000.0; *drag_p = 1000.0; }
+    assert(lb >-100.0);
+    return lb;
+}// End l_blockage
+
+
+// ************************************************************************* //
+
+double Foam::PDRutils::inters_db
+(
+    double xc, double yc, double theta,
+    double wa, double wb,
+    double x1, double x2,
+    double y1, double y2,
+    scalar* count_p,
+    symmTensor2D& vdrag, scalar* perim_p,
+    scalar* x_lblk_p, scalar* y_lblk_p,
+    scalar* x_centre_p, scalar* y_centre_p
+)
+{
+    double  x_int[6][2], y_int[6][2]; // Coordinates of intersections between the circle and sides of the rectangle.
+    double  area, lpa, lpb, len;
+
+    double  m = ::tan( theta );
+    double  cth = ::cos( theta );
+    double  sth = ::sin( theta );
+
+    double  was  = wa * sth * 0.5;
+    double  wac  = wa * cth * 0.5;
+    double  wbs  = wb * sth * 0.5;
+    double  wbc  = wb * cth * 0.5;
+    double  waos = wa / sth * 0.5;
+    double  waoc = wa / cth * 0.5;
+    double  wbos = wb / sth * 0.5;
+    double  wboc = wb / cth * 0.5;
+
+    double xb[6], yb[6], xp1, xp2, yp1, yp2;
+
+    double  dx1 = (x1 - xc);
+    double  dx2 = (x2 - xc);
+    double  dy1 = (y1 - yc);
+    double  dy2 = (y2 - yc);
+
+    *count_p = 0;
+
+// The vertices of the rectangle (all coordinates relative to centre of rectangle)
+    xb[1] = -wac - wbs;
+    yb[1] = -was + wbc;
+    xb[3] =  wac + wbs;
+    yb[3] =  was - wbc;
+    xb[2] =  wac - wbs;
+    yb[2] =  was + wbc;
+    xb[4] = -wac + wbs;
+    yb[4] = -was - wbc;
+
+    // First parameter of x_int or y_int determines which side of the cell we intersecting with
+    //  Second parameter 0 is first intersection, 1 is second, going clockwise
+
+    if ( xb[1] < dx1 )  // i.e. if left corner of block is to the left of x1
+    {
+        // Where one of lower sides of block intersects with x=x1
+        // Innermost max determines which intersection is the genuine one
+        // (not if whole block is to left of x1)
+        y_int[1][0] = min(max(max(dx1 * m - wboc, -dx1 / m - waos), dy1), dy2);
+        // Upper intersection
+        y_int[1][1] = min(max(min(dx1 * m + wboc, -dx1 / m + waos), dy1), dy2);
+    }
+    else
+    {
+        y_int[1][1] = dy1;
+        y_int[1][0] = dy2;
+        // We add a quarter to count for each vertex inside the cell
+        if ( (yb[1] > dy1) && (yb[1] < dy2) ) // ?? Seems inefficient ??
+        { *count_p += 0.25; }
+    }
+    if ( xb[3] > dx2 )
+    {
+        y_int[3][1] = min(max(max(dx2 * m - wboc, -dx2 / m - waos), dy1), dy2);
+        y_int[3][0] = min(max(min(dx2 * m + wboc, -dx2 / m + waos), dy1), dy2);
+    }
+    else
+    {
+        y_int[3][0] = dy1;
+        y_int[3][1] = dy2;
+        if (yb[3] > dy1 && yb[3] < dy2)
+        {
+            *count_p += 0.25;
+        }
+    }
+    if (yb[2] > dy2)
+    {
+        x_int[2][0] = min(max(max(dy2 / m - wbos, -dy2 * m - waoc), dx1), dx2);
+        x_int[2][1] = min(max(min(dy2 / m + wbos, -dy2 * m + waoc), dx1), dx2);
+    }
+    else
+    {
+        x_int[2][0] = dx2;
+        x_int[2][1] = dx1;
+        if ( (xb[2] > dx1) && (xb[2] < dx2) )
+        { *count_p += 0.25; }
+    }
+    if ( yb[4] < dy1 )
+    {
+        x_int[4][1] = min(max(max(dy1 / m - wbos, -dy1 * m - waoc ), dx1), dx2);
+        x_int[4][0] = min(max(min(dy1 / m + wbos, -dy1 * m + waoc ), dx1), dx2);
+    }
+    else
+    {
+        x_int[4][1] = dx2;
+        x_int[4][0] = dx1;
+        if ( (xb[4] > dx1) && (xb[4] < dx2) )
+        { *count_p += 0.25; }
+    }
+
+    y_int[0][0] = y_int[0][1] = dy1;
+    x_int[0][0] = x_int[4][0];
+    x_int[0][1] = x_int[4][1];
+    x_int[5][0] = x_int[5][1] = dx1;
+    y_int[5][0] = y_int[1][0];
+    y_int[5][1] = y_int[1][1];
+
+
+// We can now define a smaller enclosing rectangle
+
+    xp1 = min(x_int[2][0], x_int[4][1]);  // Leftmost of the intersections with top and bottom of cell
+    if ( yb[1] > dy1 && yb[1] < dy2 )       xp1 = min(xp1, xb[1] ); // left corner of block
+    xp1 = max(xp1, dx1);  // Make sure it is not to the left of x1
+
+    yp2 = max(y_int[1][1], y_int[3][0] );
+    if ( xb[2] > dx1 && xb[2] < dx2 )  yp2 = max(yp2, yb[2] );
+    yp2 = min(yp2, dy2);
+
+    xp2 = max(x_int[2][1], x_int[4][0] );
+    if ( yb[3] > dy1 && yb[3] < dy2 )  xp2 = max(xp2, xb[3] );
+    xp2 = min(xp2, dx2);
+
+    yp1 = min(y_int[1][0], y_int[3][1]);
+    if ( xb[4] > dx1 && xb[4] < dx2 )  yp1 = min(yp1, yb[4] );
+    yp1 = max(yp1, dy1 );
+
+    // Conveniently, the dimensions of the enclosing rectangle give us the line blockages
+    *x_lblk_p = (xp2 - xp1 ) / (x2 - x1 );
+    if ( *x_lblk_p < 0.0 ) { *x_lblk_p = 0.0; *count_p = 0.0; };  // ?? Better to trap no intersection earlier??
+    *y_lblk_p = (yp2 - yp1 ) / (y2 - y1 );
+    if ( *y_lblk_p < 0.0 ) { *y_lblk_p = 0.0; *count_p = 0.0; };
+
+    *x_centre_p = xc + (xp2 + xp1 ) * 0.5;
+    *y_centre_p = yc + (yp2 + yp1 ) * 0.5;
+
+    *perim_p = lpa = lpb = 0.0;;
+    area = (xp2 - xp1 ) * ( yp2 - yp1 );
+    {
+        double  dxx, dyy;
+        // Lower left
+        dyy = max(0.0, min(yb[1], y_int[1][0]) - yp1);
+        dxx = min(xb[4], x_int[0][1] ) - xp1;
+
+        if ( ( dxx * dyy) > 0.0 )
+        {
+            area -= dxx * dyy * 0.5;
+            len = std::hypot(dxx, dyy);
+            lpa += len * 0.5;
+            *perim_p += len;
+        }
+        // Upper left
+        dxx = max(0.0, min(xb[2], x_int[2][0]) - xp1);
+        dyy = yp2 - max(yb[1], y_int[1][1] );
+        if ( ( dxx * dyy) > 0.0 )
+        {
+            area -= dxx * dyy * 0.5;
+            len = std::hypot(dxx, dyy);
+            lpb += len * 0.5;
+            *perim_p += len;
+        }
+        // Upper right
+        dyy = max(0.0, yp2 - max(yb[3], y_int[3][0]));
+        dxx = xp2 - max(xb[2], x_int[2][1] );
+        if ( ( dxx * dyy) > 0.0 )
+        {
+            area -= dxx * dyy * 0.5;
+            len = std::hypot(dxx, dyy);
+            lpa += len * 0.5;
+            *perim_p += len;
+        }
+        // Lower right
+        dxx = max(0.0, xp2 - max(xb[4], x_int[4][0]));
+        dyy = min(yb[3], y_int[3][1] ) - yp1;
+        if ( ( dxx * dyy) > 0.0 )
+        {
+            area -= dxx * dyy * 0.5;
+            len = std::hypot(dxx, dyy);
+            lpb += len * 0.5;
+            *perim_p += len;
+        }
+
+    }
+
+    vdrag.xx() = lpa * cth * cth + lpb * sth * sth;
+    vdrag.xy() = lpa * cth * sth - lpb * sth * cth;
+    vdrag.yy() = lpa * sth * sth + lpb * cth * cth;
+
+    return area / ( (x2 - x1 ) * ( y2 - y1 ) );
+} // End inters_db
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/PDRutilsOverlap.C

@@ -0,0 +1,767 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRsetFields.H"
+#include "PDRutilsInternal.H"
+#include "mathematicalConstants.H"
+
+#ifndef FULLDEBUG
+#define NDEBUG
+#endif
+#include <cassert>
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    // A sign-corrected multiply
+    // This is used for porosity of obstacle intersections
+    inline static scalar COMBLK(const scalar a, const scalar b)
+    {
+        if (a < 0)
+        {
+            return -a * b;
+        }
+
+        return a * b;
+    }
+
+
+    // Obstacle satisfies some minimum size checks.
+    // A volume check misses thin plates, so use area.
+    // Thin sheet overlaps can be produced by touching objects
+    // if the obs_extend parameter is > 0.
+    inline static bool obsHasMinSize(const vector& span, const PDRparams& tol)
+    {
+        return
+        (
+            (cmptProduct(span) > tol.min_overlap_vol)
+         &&
+            (
+                (span.x() * span.y() > tol.min_overlap_area)
+             || (span.y() * span.z() > tol.min_overlap_area)
+             || (span.z() * span.x() > tol.min_overlap_area)
+            )
+        );
+    }
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+void Foam::PDRutils::one_d_overlap
+(
+    scalar xmin,
+    scalar xmax,
+    const PDRblock::location& grid,
+    List<scalar>& olap,
+    int *cmin, int *cmax,
+    int *cfmin, int *cfmax
+)
+{
+    // Looking at one coordinate direction, called x here, for something
+    // that extends from xmin to xmax, calculate how much it overlaps
+    // each cell in this direction. Result returned in 'olap' List is
+    // the proportion of the grid step overlapped, i.e dimensionless.
+    // First and last steps overlapped given by *cmin, *cmax
+    // Ditto for shifted grid given by *cfmin, *cfmax.
+
+    // Initially zero everwhere
+    olap = Zero;
+
+    if (olap.size() < grid.nPoints())
+    {
+        FatalErrorInFunction
+            << "The overlap scratch array is too small, has "
+            << olap.size() << " but needs " << grid.nPoints() << nl
+            << exit(FatalError);
+    }
+
+
+    // No intersection with the box
+    if (xmax <= grid.first() || grid.last() <= xmin)
+    {
+        // Mark as bad range, cannot iterate
+        *cmin = 0;
+        *cmax = -1;
+
+        // Another bad range (cannot iterate) but for extra safety ensure
+        // that (cfmin -> cmin) and (cmax -> cfmax) cannot iterate either
+        *cfmin = 1;
+        *cfmax = -2;
+        return;
+    }
+
+    // Ensure search is within the (point) bounds
+    xmin = grid.clip(xmin);
+    xmax = grid.clip(xmax);
+
+    // The begin/end of the obstacle
+    *cmin = grid.findCell(xmin);
+    *cmax = grid.findCell(xmax);
+
+    for (label ix = *cmin; ix <= *cmax; ++ix)
+    {
+        olap[ix] = 1.0;
+    }
+
+    // Fixup ends
+    if (*cmax == *cmin)
+    {
+        olap[*cmax] = (xmax - xmin) / grid.width(*cmax);
+    }
+    else
+    {
+        if (grid[*cmin] < xmin)
+        {
+            olap[*cmin] = (grid[*cmin+1] - xmin) / grid.width(*cmin);
+        }
+
+        if (xmax < grid[*cmax+1])
+        {
+            olap[*cmax] = (xmax - grid[*cmax]) / grid.width(*cmax);
+        }
+    }
+    assert(olap[*cmax] >= 0.0);
+
+
+    // Is xmin below/above the cell-centre (for virtual staggered-grid) ?
+    *cfmin =
+    (
+        xmin < grid.C(*cmin)
+      ? *cmin
+      : Foam::min(*cmin+1, grid.nCells()-1)
+    );
+
+    // Is xmax below/above the cell-centre (for virtual staggered-grid) ?
+    *cfmax =
+    (
+        xmax < grid.C(*cmax)
+      ? *cmax
+      : Foam::min(*cmax+1, grid.nCells()-1)
+    );
+}
+
+
+/**************************************************************************************************/
+
+void Foam::PDRutils::two_d_overlap
+(
+    const UList<scalar>& a_olap, label amin, label amax,
+    const UList<scalar>& b_olap, label bmin, label bmax,
+    SquareMatrix<scalar>& ab_olap
+)
+{
+    // We go one over the relevant min/max limits since these values might be
+    // used. If not, they would have been zeroed in one_d_overlap
+
+    amin = Foam::max(0, amin-1);
+    bmin = Foam::max(0, bmin-1);
+    amax = Foam::min(a_olap.size()-1, amax+1);
+    bmax = Foam::min(b_olap.size()-1, bmax+1);
+
+    for (label ia = amin; ia <= amax; ++ia)
+    {
+        for (label ib = bmin; ib <= bmax; ++ib)
+        {
+            ab_olap(ia,ib) = a_olap[ia] * b_olap[ib];
+        }
+    }
+}
+
+
+/**************************************************************************************************/
+
+void Foam::PDRutils::circle_overlap
+(
+    scalar ac, scalar bc, scalar dia,
+    scalar theta, scalar wa, scalar wb,
+    const PDRblock::location& agrid, label amin, label amax,
+    const PDRblock::location& bgrid, label bmin, label bmax,
+    SquareMatrix<scalar>& ab_olap,
+    SquareMatrix<scalar>& ab_perim,
+    SquareMatrix<scalar>& a_lblock,
+    SquareMatrix<scalar>& ac_lblock,
+    SquareMatrix<scalar>& c_count,
+    SquareMatrix<symmTensor2D>& c_drag,
+    SquareMatrix<scalar>& b_lblock,
+    SquareMatrix<scalar>& bc_lblock
+)
+{
+    /* This routine calculates the proportion of each (two-dimensional) grid cell
+       overlapped by the circle or angled rectangle. Coordinates are labelled a and b.
+       On entry:
+              ac, bc                coordinates of centre of circle or rectangle
+              dia                   diameter of circle (zero for rectangle)
+              theta, wa, wb parameters for rectangle
+              agrid[]               locations of grid lines of a-grid
+              amin, amax    first and last cells in a-grid overlapped by object
+              (similarly for b)
+            On exit:
+              abolap                2-D array of (proportionate) area blockage by grid cell
+              a_lblock              2-D array of (proportionate) blockage to a-direction flow
+                                            (This will be area blockage when extruded in the third coordinate).
+              a_count               (2-D array)The contribution of this object to the count of obstacles blocking
+                                            a-direction flow. This is only non-zero if the object is inside the
+                                            lateral boundaries of the cell. It is large negative if the cell is
+                                            totally blocked in this direction.
+              (similarly for b)
+              c_drag                2-D array of tensor that will give tensor drag in each cell (when multiplied
+                                            Cd, cylinder length, and 0.5 rho*U^2) Dimension: L.
+
+       Note that this routine does not zero array elements outside the amin to amax, bmin to bmax area.
+    */
+
+    scalar count, a_lblk, b_lblk, perim, dummy;
+
+    symmTensor2D vdrag(Zero);
+
+    // Prevent stepping outside of the array when the obstacle is on the
+    // upper boundary
+
+    // Upper limit of inclusive range is nCells-1
+    amin = Foam::max(0, amin);
+    bmin = Foam::max(0, bmin);
+    amax = Foam::min(amax, agrid.nCells()-1);
+    bmax = Foam::min(bmax, bgrid.nCells()-1);
+
+    for (label ia = amin; ia <= amax; ++ia)
+    {
+        // Cell-centred grid
+        const scalar a1 = agrid[ia];
+        const scalar a2 = agrid[ia+1];
+
+        // Left-shifted staggered face grid (-1 addressing is OK)
+        const scalar af1 = agrid.C(ia-1);
+        const scalar af2 = agrid.C(ia);
+
+        for (label ib = bmin; ib <= bmax; ++ib)
+        {
+            // Cell-centred grid
+            const scalar b1 = bgrid[ib];
+            const scalar b2 = bgrid[ib+1];
+
+            // Left-shifted staggered face grid (-1 addressing is OK)
+            const scalar bf1 = bgrid.C(ib-1);
+            const scalar bf2 = bgrid.C(ib);
+
+            // Do the centred cell
+            if ( dia > 0.0 )
+            {
+                ab_olap(ia,ib) = inters_cy
+                (
+                    ac, bc, 0.5*dia, a1, a2, b1, b2, &perim,
+                    &dummy, &dummy, &b_lblk, &a_lblk
+                );
+/* The last two arguments of the above call appear to be reversed, but the inters_cy routine returns
+   the amount of overlap in the a and b direcvtions, which are the blockage to the b and a directions. */
+
+/* abolap * cell area is area of cylinder in this cell. Divide by PI%D^2/4 to get proportion of cylinder in cell
+   For whole cylinger c_drag should be = D, so multiply by D.                                                   */
+
+                c_drag(ia,ib).xx() = c_drag(ia,ib).yy() = 4.0 * ab_olap(ia,ib) * (a2 - a1) * (b2 - b1) / dia / mathematical::pi;
+                c_drag(ia,ib).xy() = Zero;
+                c_count(ia,ib) = perim / (mathematical::pi * dia);
+
+//******?????
+                scalar area = (a2 - a1) * (b2 - b1);
+                scalar rat = dia * dia / area - 1.5;
+                if (rat > 0.0)
+                {
+                    scalar da = ac - 0.5 * (a1 + a2);
+                    scalar db = bc - 0.5 * (b1 + b2);
+                    scalar dc = std::hypot(da, db);
+                    scalar rat1 = min(max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
+                    scalar drg0 = c_drag(ia,ib).xx();
+                    scalar drg1 = c_drag(ia,ib).yy();
+                    scalar drg = std::hypot(drg0, drg1);
+                    c_drag(ia,ib).xx() = drg * ( 1.0 - rat1 ) + drg * da*da/dc/dc * rat1;
+                    c_drag(ia,ib).yy() = drg * ( 1.0 - rat1 ) + drg * db*db/dc/dc * rat1;
+                    c_drag(ia,ib).xy() = drg * da*db/dc/dc *rat1;
+                }
+            }
+            else
+            {
+                ab_olap(ia,ib) = inters_db( ac, bc, theta, wa, wb, a1, a2, b1, b2, &count, c_drag(ia,ib), &perim, &a_lblk, &b_lblk, &dummy, &dummy );
+                c_count(ia,ib) = perim / ( wa + wb ) * 0.5;
+            }
+            ac_lblock(ia,ib) = a_lblk;
+            bc_lblock(ia,ib) = b_lblk;
+            ab_perim(ia,ib) = perim;
+
+            // Do the a-shifted cell
+            if ( dia > 0.0 ) // I.e. a cylinder, not a d.b.
+            {
+                if (ac >= af1 && ac < af2)
+                {
+                    // Only want to block one layer of faces
+                    a_lblock(ia,ib) = l_blockage
+                    (
+                        ac, bc, 0.5*dia,
+                        af1, af2, b1, b2, &count, &dummy, &dummy
+                    );
+                }
+                inters_cy
+                (
+                    ac, bc, 0.5*dia,
+                    af1, af2, b1, b2,
+                    &perim, &count, &dummy, &dummy, &dummy
+                );
+            }
+            else
+            {
+                inters_db
+                (
+                    ac, bc, theta, wa, wb, af1, af2, b1, b2,
+                    &count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
+                );
+                a_lblock(ia,ib) = a_lblk;
+            }
+
+            // Do the b-shifted cell
+            if ( dia > 0.0 )
+            {
+                if (bc >= bf1 && bc < bf2)
+                {
+                    // Only want to block one layer of faces
+                    b_lblock(ia,ib) = l_blockage
+                    (
+                        bc, ac, 0.5*dia, bf1, bf2, a1, a2,
+                        &count, &(vdrag.yy()), &dummy
+                    );
+                }
+
+                inters_cy
+                (
+                    ac, bc, 0.5*dia,
+                    a1, a2, bf1, bf2,
+                    &perim, &dummy, &count, &dummy, &dummy
+                );
+            }
+            else
+            {
+                inters_db
+                (
+                    ac, bc, theta, wa, wb,
+                    a1, a2, bf1, bf2,
+                    &count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
+                );
+                b_lblock(ia,ib) =  b_lblk;
+            }
+        }
+    }
+
+}  // End circle_overlap
+
+
+/**************************************************************************************************/
+
+scalar block_overlap
+(
+    DynamicList<PDRobstacle>& blocks,
+    const labelRange& range,
+    const scalar multiplier
+)
+{
+    // Size information
+    const label nBlock = range.size();
+
+    // The return value
+    scalar totVolume = 0;
+
+    if (nBlock < 2) return 0;
+
+
+    // Sort blocks by their x-position (with sortBias)
+    labelList blkOrder;
+    sortedOrder(blocks[range], blkOrder);
+
+    DynamicList<PDRobstacle> newBlocks;
+
+    // Work through the sorted blocks
+    for (label i1 = 0; i1 < nBlock-1; ++i1)
+    {
+        const PDRobstacle& blk1 = blocks[range[blkOrder[i1]]];
+
+        // Upper coordinates
+        const vector max1 = blk1.pt + blk1.span;
+
+        // For second block start with the next one on the list, and
+        // stop when we find the first one whose biased x-position
+        // is beyond the end of the block1
+
+        for (label i2 = i1 + 1; i2 < nBlock; ++i2)
+        {
+            const PDRobstacle& blk2 = blocks[range[blkOrder[i2]]];
+
+            // Upper coordinates
+            const vector max2 = blk2.pt + blk2.span;
+
+            if (max1.x() <= blk2.x())
+            {
+                break;
+            }
+
+            if
+            (
+                max1.y() <= blk2.y()
+             || max1.z() <= blk2.z()
+             || max2.y() <= blk1.y()
+             || max2.z() <= blk1.z()
+             || (blk1.vbkge * blk2.vbkge <= 0)
+            )
+            {
+                continue;
+            }
+
+
+            {
+                PDRobstacle over;
+
+                over.pt = max(blk1.pt, blk2.pt);
+                over.span = min(max1, max2) - over.pt;
+
+                assert(cmptProduct(over.span) > 0.0);
+
+                // This routine should only have been called for all +ve o r all -ve obstacles
+                assert(blk1.vbkge * blk2.vbkge > 0);
+                /* At the first level of intersection, we create an obstacle of blockage -1 (if both objects solid)
+                 to cancel out the double counting. (multiplier is 1).
+                 ?? COMBLK does a (sign corrected) multiply; is this corrrect for porous obstacles?
+                 Depends on how blockages were summed in the first place. In fact this -ve obstacle
+                 concept only works if the blockages are summed??*/
+                over.vbkge = - COMBLK( blk1.vbkge, blk2.vbkge ) * multiplier;
+                over.xbkge = - COMBLK( blk1.xbkge, blk2.xbkge ) * multiplier;
+                over.ybkge = - COMBLK( blk1.ybkge, blk2.ybkge ) * multiplier;
+                over.zbkge = - COMBLK( blk1.zbkge, blk2.zbkge ) * multiplier;
+                over.typeId = 81 + int(15 * multiplier);  // Not subsequently used
+
+                if (obsHasMinSize(over.span, pars))
+                {
+                    // Obstacle satisfies some minimum size checks
+                    totVolume -= over.volume();
+
+                    newBlocks.append(over);
+                }
+            }
+        }
+    }
+
+    blocks.append(std::move(newBlocks));
+
+    return totVolume;
+}
+
+
+/**************************************************************************************************/
+
+using namespace Foam::PDRutils;
+
+scalar block_cylinder_overlap
+(
+    DynamicList<PDRobstacle>& blocks,
+    const labelRange& range,
+    const UList<PDRobstacle>& cylinders
+)
+{
+    // Size information
+    const label nBlock = range.size();
+    const label nCyl = cylinders.size();
+
+    // The return value
+    scalar totVolume = 0;
+
+    if (!nBlock || !nCyl) return 0;
+
+    scalar area, a_lblk, b_lblk, dummy, a_centre, b_centre;
+    symmTensor2D dum2;
+
+
+    // Sort blocks and cylinders by their x-position (with sortBias)
+    labelList blkOrder;
+    sortedOrder(blocks[range], blkOrder);
+
+    labelList cylOrder;
+    sortedOrder(cylinders, cylOrder);
+
+    DynamicList<PDRobstacle> newBlocks;
+
+    // Work through the sorted blocks
+    for (label i1 = 0; i1 < nBlock; i1++)
+    {
+        const PDRobstacle& blk1 = blocks[range[blkOrder[i1]]];
+
+        // Upper coordinates
+        const vector max1 = blk1.pt + blk1.span;
+
+        // Cyls whose end is before start of this block no longer
+        // need to be considered
+
+        label i2 = 0;
+        while (i2 < nCyl-1 && cylinders[cylOrder[i2]] < blk1)
+        {
+            ++i2;
+        }
+
+        for (/*nil*/; i2 < nCyl; ++i2)
+        {
+            const PDRobstacle& cyl2 = cylinders[cylOrder[i2]];
+
+            // Calculate overlap in axis direction; if zero continue.
+            // Calculate 2-d overlap and c 0f g; if area zero continue.
+
+            PDRobstacle over;
+
+
+            switch (cyl2.orient)
+            {
+                case vector::Z:
+                {
+                    const scalar zm2 = cyl2.z() + cyl2.len();
+                    if (blk1.z() > zm2 || cyl2.z() > max1.z()) continue;
+
+                    if ( cyl2.dia() == 0.0 )
+                    {
+                        area = inters_db
+                        (
+                            cyl2.x(), cyl2.y(), cyl2.theta(), cyl2.wa, cyl2.wb,
+                            blk1.x(), max1.x(),
+                            blk1.y(), max1.y(),
+                            &dummy, dum2, &dummy, &a_lblk, &b_lblk,
+                            &a_centre, &b_centre
+                        );
+                    }
+                    else
+                    {
+                        area = inters_cy
+                        (
+                            cyl2.x(), cyl2.y(), 0.5*cyl2.dia(),
+                            blk1.x(), max1.x(),
+                            blk1.y(), max1.y(),
+                            &dummy, &dummy, &dummy, &dummy, &dummy
+                        );
+                        b_lblk = l_blockage
+                        (
+                            cyl2.x(), cyl2.y(), 0.5*cyl2.dia(),
+                            blk1.x(), max1.x(),
+                            blk1.y(), max1.y(),
+                            &dummy, &dummy, &b_centre
+                        );
+                        a_lblk = l_blockage
+                        (
+                            cyl2.y(), cyl2.x(), 0.5*cyl2.dia(),
+                            blk1.y(), max1.y(),
+                            blk1.x(), max1.x(),
+                            &dummy, &dummy, &a_centre
+                        );
+                    }
+                    if (equal(area, 0)) continue;
+                    assert(a_lblk >0.0);
+                    assert(b_lblk >0.0);
+
+                    // The intersection between a circle and a rectangle  can be an odd shape.
+                    // We have its area. a_lblk and b_lblk are dimensions of enclosing rectangle
+                    // and a_centre and b_centre its centre. We scale this rectangle down to
+                    // the corect areacorrect area, as a rectangular approximation to the intersection.
+                    const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
+
+                    a_lblk *= blk1.span.x() * ratio;
+                    b_lblk *= blk1.span.y() * ratio;
+                    assert(b_lblk >0.0);
+                    assert(a_lblk >0.0);
+
+                    over.x() = a_centre - 0.5 * a_lblk;
+                    over.y() = b_centre - 0.5 * b_lblk;
+                    over.z() = max(blk1.z(), cyl2.z());
+
+                    over.span.x() = a_lblk;
+                    over.span.y() = b_lblk;
+                    over.span.z() = min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
+                    assert(over.x() > -200.0);
+                    assert(over.x() < 2000.0);
+                }
+                break;
+
+                case vector::Y:
+                {
+                    const scalar ym2 = cyl2.y() + cyl2.len();
+                    if (blk1.y() > ym2 || cyl2.y() > max1.y()) continue;
+
+                    if ( cyl2.dia() == 0.0 )
+                    {
+                        area = inters_db
+                        (
+                            cyl2.z(), cyl2.x(), cyl2.theta(), cyl2.wa, cyl2.wb,
+                            blk1.z(), max1.z(),
+                            blk1.x(), max1.x(),
+                            &dummy, dum2, &dummy, &a_lblk, &b_lblk,
+                            &a_centre, &b_centre
+                        );
+                    }
+                    else
+                    {
+                        area = inters_cy
+                        (
+                            cyl2.z(), cyl2.x(), 0.5*cyl2.dia(),
+                            blk1.z(), max1.z(),
+                            blk1.x(), max1.x(),
+                            &dummy, &dummy, &dummy, &dummy, &dummy
+                        );
+
+                        b_lblk = l_blockage
+                        (
+                            cyl2.z(), cyl2.x(), 0.5*cyl2.dia(),
+                            blk1.z(), max1.z(),
+                            blk1.x(), max1.x(),
+                            &dummy, &dummy, &b_centre
+                        );
+
+                        a_lblk = l_blockage
+                        (
+                            cyl2.x(), cyl2.z(), 0.5*cyl2.dia(),
+                            blk1.x(), max1.x(),
+                            blk1.z(), max1.z(),
+                            &dummy, &dummy, &a_centre
+                        );
+                    }
+
+                    if (equal(area, 0)) continue;
+                    assert(a_lblk >0.0);
+                    assert(b_lblk >0.0);
+
+                    // a_lblk and b_lblk are dimensions of enclosing rectangle.
+                    // Need to scale to correct area
+                    const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
+                    a_lblk *= blk1.span.z() * ratio;
+                    b_lblk *= blk1.span.x() * ratio;
+
+                    over.z() = a_centre - a_lblk * 0.5;
+                    over.x() = b_centre - b_lblk * 0.5;
+                    over.y() = max(blk1.y(), cyl2.y());
+
+                    over.span.z() = a_lblk;
+                    over.span.x() = b_lblk;
+                    over.span.y() = min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
+                }
+                break;
+
+                case vector::X:
+                {
+                    const scalar xm2 = cyl2.x() + cyl2.len();
+                    if (blk1.x() > xm2 || cyl2.x() > max1.x()) continue;
+
+                    if ( cyl2.dia() == 0.0 )
+                    {
+                        area = inters_db
+                        (
+                            cyl2.y(), cyl2.z(), cyl2.theta(), cyl2.wa, cyl2.wb,
+                            blk1.y(), max1.y(),
+                            blk1.z(), max1.z(),
+                            &dummy, dum2, &dummy, &a_lblk, &b_lblk,
+                            &a_centre, &b_centre
+                        );
+                    }
+                    else
+                    {
+                        area = inters_cy
+                        (
+                            cyl2.y(), cyl2.z(), 0.5*cyl2.dia(),
+                            blk1.y(), max1.y(),
+                            blk1.z(), max1.z(),
+                            &dummy, &dummy, &dummy, &dummy, &dummy
+                        );
+
+                        b_lblk = l_blockage
+                        (
+                            cyl2.y(), cyl2.z(), 0.5*cyl2.dia(),
+                            blk1.y(), max1.y(),
+                            blk1.z(), max1.z(),
+                            &dummy, &dummy, &b_centre
+                        );
+
+                        a_lblk = l_blockage
+                        (
+                            cyl2.z(), cyl2.y(), 0.5*cyl2.dia(),
+                            blk1.z(), max1.z(),
+                            blk1.y(), max1.y(),
+                            &dummy, &dummy, &a_centre
+                        );
+
+                    }
+
+                    if (equal(area, 0)) continue;
+                    assert(a_lblk >0.0);
+                    assert(b_lblk >0.0);
+
+                    // a_lblk and b_lblk are dimensions of enclosing rectangle.
+                    // Need to scale to correct area
+                    const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
+                    assert(ratio >-10000.0);
+                    assert(ratio <10000.0);
+                    a_lblk *= blk1.span.y() * ratio;
+                    b_lblk *= blk1.span.z() * ratio;
+
+                    over.y() = a_centre - a_lblk * 0.5;
+                    over.z() = b_centre - b_lblk * 0.5;
+                    over.x() = max(blk1.x(), cyl2.x());
+
+                    over.span.y() = a_lblk;
+                    over.span.z() = b_lblk;
+                    over.span.x() = min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
+                }
+                break;
+            }
+            over.vbkge = over.xbkge = over.ybkge = over.zbkge = -1.0;
+            over.typeId = PDRobstacle::IGNORE;
+
+            assert(cmptProduct(over.span) > 0.0);
+            assert(b_lblk >0.0);
+            assert(a_lblk >0.0);
+            assert(over.x() > -10000.0);
+
+            if (obsHasMinSize(over.span, pars))
+            {
+                // Obstacle satisfies some minimum size checks
+                totVolume -= over.volume();
+
+                newBlocks.append(over);
+            }
+        }
+    }
+
+    blocks.append(std::move(newBlocks));
+
+    return totVolume;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/ObstaclesDict

@@ -0,0 +1,179 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      obstaclesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+group1
+{
+    positions
+    (
+        (0 0 0)
+    );
+
+    obstacles
+    (
+        box { point (0 0 0); span (0.05 0.05 2); porosity 0; }
+        box { point (1 0 0); span (0.05 0.05 2); porosity 0; }
+        box { point (1 0 0); span (0.05 0.05 2); porosity 0; }
+        box { point (2 0 0); span (0.05 0.05 2); porosity 0; }
+        box { point (3 0 0); span (0.05 0.05 2); porosity 0; }
+        box { point (0 1 0); span (0.05 0.05 2); porosity 0; }
+        box { point (1 1 0); span (0.05 0.05 2); porosity 0; }
+        box { point (2 1 0); span (0.05 0.05 2); porosity 0; }
+        box { point (3 1 0); span (0.05 0.05 2); porosity 0; }
+        box { point (0 2 0); span (0.05 0.05 2); porosity 0; }
+        box { point (1 2 0); span (0.05 0.05 2); porosity 0; }
+        box { point (2 2 0); span (0.05 0.05 2); porosity 0; }
+        box { point (3 2 0); span (0.05 0.05 2); porosity 0; }
+        box { point (0 3 0); span (0.05 0.05 2); porosity 0; }
+        box { point (1 3 0); span (0.05 0.05 2); porosity 0; }
+        box { point (2 3 0); span (0.05 0.05 2); porosity 0; }
+        box { point (3 3 0); span (0.05 0.05 2); porosity 0; }
+        box { point (0 0 0); span (0.05 3.05 0.05); }
+        box { point (1 0 0); span (0.05 3.05 0.05); }
+        box { point (2 0 0); span (0.05 3.05 0.05); }
+        box { point (3 0 0); span (0.05 3.05 0.05); }
+        box { point (0 0 1); span (0.05 3.05 0.05); }
+        box { point (1 0 1); span (0.05 3.05 0.05); }
+        box { point (2 0 1); span (0.05 3.05 0.05); }
+        box { point (3 0 1); span (0.05 3.05 0.05); }
+        box { point (0 0 2); span (0.05 3.05 0.05); }
+        box { point (1 0 2); span (0.05 3.05 0.05); }
+        box { point (2 0 2); span (0.05 3.05 0.05); }
+        box { point (3 0 2); span (0.05 3.05 0.05); }
+        box { point (0 0 0); span (3.05 0.05 0.05); }
+        box { point (0 1 0); span (3.05 0.05 0.05); }
+        box { point (0 2 0); span (3.05 0.05 0.05); }
+        box { point (0 3 0); span (3.05 0.05 0.05); }
+        box { point (0 0 1); span (3.05 0.05 0.05); }
+        box { point (0 1 1); span (3.05 0.05 0.05); }
+        box { point (0 2 1); span (3.05 0.05 0.05); }
+        box { point (0 3 1); span (3.05 0.05 0.05); }
+        box { point (0 0 2); span (3.05 0.05 0.05); }
+        box { point (0 1 2); span (3.05 0.05 0.05); }
+        box { point (0 2 2); span (3.05 0.05 0.05); }
+        box { point (0 3 2); span (3.05 0.05 0.05); }
+
+        patch { point (1 1 1); span (2 2 2); inlet x; outlet y; name xpatch; }
+    );
+}
+
+group11
+{
+    positions
+    (
+        (0 0 0)
+    );
+
+}
+
+group14
+{
+    positions
+    (
+        (0 0 0.15)
+        (0 0 0.45)
+    );
+
+    obstacles
+    (
+        box { point (0.05   0  1.05); span(0.006 3.05 0.05); }
+        box { point (0.997  0  1.05); span(0.006 3.05 0.05); }
+        box { point (1.05   0  1.05); span(0.006 3.05 0.05); }
+        box { point (1.997  0  1.05); span(0.006 3.05 0.05); }
+        box { point (2.05   0  1.05); span(0.006 3.05 0.05); }
+        box { point (2.997  0  1.05); span(0.006 3.05 0.05); }
+
+        cyl { point (0.05  0.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  0.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  0.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  0.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  1.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  1.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  1.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  1.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  2.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  2.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  2.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  2.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (0.05  3.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  0.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  0.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  0.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  0.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  1.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  1.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  1.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  1.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  2.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  2.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  2.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  2.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (1.05  3.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  0.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  0.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  0.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  0.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  1.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  1.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  1.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  1.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  2.025  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  2.275  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  2.525  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  2.775  1.075); length 0.947; diameter 0.026; direction x; }
+        cyl { point (2.05  3.025  1.075); length 0.947; diameter 0.026; direction x; }
+    );
+}
+
+
+group17
+{
+    positions
+    (
+        (0  0  0)
+        (0  0  0.30)
+    );
+
+}
+
+group21
+{
+    positions
+    (
+        (0  0  0)
+    );
+
+}
+
+group31
+{
+    positions
+    (
+        (0  0  0)
+    );
+
+}
+
+group41
+{
+    positions
+    (
+        (0  0  0)
+    );
+
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacle.C

@@ -0,0 +1,737 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRobstacle.H"
+#include "boundBox.H"
+#include "meshedSurf.H"
+#include "axisAngleRotation.H"
+#include "coordinateSystem.H"
+#include "foamVtkSurfaceWriter.H"
+#include "unitConversion.H"
+#include "addToMemberFunctionSelectionTable.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineMemberFunctionSelectionTable(PDRobstacle, read, dictRead);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PDRobstacle::PDRobstacle()
+:
+    groupId(0),
+    typeId(0),
+    orient(vector::X),
+    sortBias(0),
+    pt(Zero),
+    span(Zero),
+    wa(0),
+    wb(0),
+    vbkge(0),
+    xbkge(0),
+    ybkge(0),
+    zbkge(0),
+    blowoff_type(0),
+    identifier()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PDRobstacle::clear()
+{
+    groupId = 0;
+    typeId = 0;
+    orient = vector::X;
+    sortBias = 0;
+    pt = Zero;
+    span = Zero;
+    wa = 0;
+    wb = 0;
+    vbkge = 0;
+    xbkge = 0;
+    ybkge = 0;
+    zbkge = 0;
+    blowoff_type = 0;
+    identifier.clear();
+}
+
+
+void Foam::PDRobstacle::readProperties(const dictionary& dict)
+{
+    PDRobstacle::clear();
+
+    // Read as word, which handles quoted or unquoted entries
+    word obsName;
+
+    if (dict.readIfPresent("name", obsName))
+    {
+        identifier = std::move(obsName);
+    }
+}
+
+
+void Foam::PDRobstacle::scale(const scalar factor)
+{
+    if (factor <= 0)
+    {
+        return;
+    }
+
+    sortBias *= factor;
+
+    switch (typeId)
+    {
+        case PDRobstacle::CYLINDER:
+        {
+            pt *= factor;
+
+            dia() *= factor;
+            len() *= factor;
+            break;
+        }
+
+        case PDRobstacle::DIAG_BEAM:
+        {
+            pt *= factor;
+
+            len() *= factor;
+            wa *= factor;
+            wb *= factor;
+            break;
+        }
+
+        case PDRobstacle::CUBOID_1:
+        case PDRobstacle::LOUVRE_BLOWOFF:
+        case PDRobstacle::CUBOID:
+        case PDRobstacle::WALL_BEAM:
+        case PDRobstacle::GRATING:
+        case PDRobstacle::RECT_PATCH:
+        {
+            pt *= factor;
+            span *= factor;
+
+            if (typeId == PDRobstacle::GRATING)
+            {
+                slat_width *= factor;
+            }
+            break;
+        }
+    }
+}
+
+
+Foam::scalar Foam::PDRobstacle::volume() const
+{
+    scalar vol = 0;
+
+    switch (typeId)
+    {
+        case PDRobstacle::CYLINDER:
+            vol = 0.25 * mathematical::pi * sqr(dia()) * len();
+            break;
+
+        case PDRobstacle::DIAG_BEAM:
+            vol = wa * wb * len();
+            break;
+
+        case PDRobstacle::CUBOID_1:
+        case PDRobstacle::LOUVRE_BLOWOFF:
+        case PDRobstacle::CUBOID:
+        case PDRobstacle::WALL_BEAM:
+        case PDRobstacle::GRATING:
+            vol = cmptProduct(span) * vbkge;
+            break;
+    }
+
+    return vol;
+}
+
+
+bool Foam::PDRobstacle::tooSmall(const scalar minWidth) const
+{
+    if (minWidth <= 0)
+    {
+        return false;
+    }
+
+    switch (typeId)
+    {
+        case PDRobstacle::CYLINDER:
+        {
+            // The effective half-width
+            if ((0.25 * dia() * sqrt(mathematical::pi)) <= minWidth)
+            {
+                return true;
+            }
+            break;
+        }
+
+        case PDRobstacle::DIAG_BEAM:
+        {
+            if
+            (
+                (len() <= minWidth && wa <= minWidth)
+             || (len() <= minWidth && wb <= minWidth)
+             || (wa  <= minWidth && wb <= minWidth)
+            )
+            {
+                return true;
+            }
+            break;
+        }
+
+        case PDRobstacle::CUBOID_1:
+        case PDRobstacle::LOUVRE_BLOWOFF:
+        case PDRobstacle::CUBOID:
+        case PDRobstacle::WALL_BEAM:
+        case PDRobstacle::GRATING:
+        case PDRobstacle::RECT_PATCH:
+        {
+            if
+            (
+                (span.x() <= minWidth && span.y() <= minWidth)
+             || (span.y() <= minWidth && span.z() <= minWidth)
+             || (span.z() <= minWidth && span.x() <= minWidth)
+            )
+            {
+                return true;
+            }
+
+            break;
+        }
+    }
+
+    return false;
+}
+
+
+Foam::volumeType Foam::PDRobstacle::trim(const boundBox& bb)
+{
+    volumeType::type vt = volumeType::UNKNOWN;
+
+    if (!bb.valid() || !typeId)
+    {
+        return vt;
+    }
+
+    switch (typeId)
+    {
+        case PDRobstacle::CYLINDER:
+        {
+            const scalar rad = 0.5*dia();
+
+            direction e1 = vector::X;
+            direction e2 = vector::Y;
+            direction e3 = vector::Z;
+
+            if (orient == vector::X)
+            {
+                e1 = vector::Y;
+                e2 = vector::Z;
+                e3 = vector::X;
+            }
+            else if (orient == vector::Y)
+            {
+                e1 = vector::Z;
+                e2 = vector::X;
+                e3 = vector::Y;
+            }
+            else
+            {
+                orient = vector::Z;  // extra safety?
+            }
+
+            if
+            (
+                (pt[e1] + rad <= bb.min()[e1])
+             || (pt[e2] + rad <= bb.min()[e2])
+             || (pt[e3] + len() <= bb.min()[e3])
+             || (pt[e1] - rad >= bb.max()[e1])
+             || (pt[e2] - rad >= bb.max()[e2])
+             || (pt[e3] >= bb.max()[e3])
+            )
+            {
+                // No overlap
+                return volumeType::OUTSIDE;
+            }
+
+            vt = volumeType::INSIDE;
+
+            // Trim obstacle length as required
+            if (pt[e3] < bb.min()[e3])
+            {
+                vt = volumeType::MIXED;
+                len() -= bb.min()[e3] - pt[e3];
+                pt[e3] = bb.min()[e3];
+            }
+
+            if (pt[e3] + len() > bb.max()[e3])
+            {
+                vt = volumeType::MIXED;
+                len() = bb.max()[e3] - pt[e3];
+            }
+
+            // Cannot trim diameter very well, so just mark as protruding
+            if
+            (
+                (pt[e1] - rad < bb.min()[e1]) || (pt[e1] + rad > bb.max()[e1])
+             || (pt[e2] - rad < bb.min()[e2]) || (pt[e2] + rad > bb.max()[e2])
+            )
+            {
+                vt = volumeType::MIXED;
+            }
+
+            break;
+        }
+
+
+        case PDRobstacle::DIAG_BEAM:
+        {
+            // Not implemented
+            break;
+        }
+
+
+        case PDRobstacle::CUBOID_1:
+        case PDRobstacle::LOUVRE_BLOWOFF:
+        case PDRobstacle::CUBOID:
+        case PDRobstacle::WALL_BEAM:
+        case PDRobstacle::GRATING:
+        case PDRobstacle::RECT_PATCH:
+        {
+            for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+            {
+                if
+                (
+                    ((pt[cmpt] + span[cmpt]) < bb.min()[cmpt])
+                 || (pt[cmpt] > bb.max()[cmpt])
+                )
+                {
+                    // No overlap
+                    return volumeType::OUTSIDE;
+                }
+            }
+
+
+            vt = volumeType::INSIDE;
+
+            // Trim obstacle as required
+
+            for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+            {
+                if (pt[cmpt] < bb.min()[cmpt])
+                {
+                    vt = volumeType::MIXED;
+                    if (span[cmpt] > 0)
+                    {
+                        span[cmpt] -= bb.min()[cmpt] - pt[cmpt];
+                    }
+                    pt[cmpt] = bb.min()[cmpt];
+                }
+
+
+                if (pt[cmpt] + span[cmpt] > bb.max()[cmpt])
+                {
+                    vt = volumeType::MIXED;
+                    span[cmpt] -= bb.max()[cmpt] - pt[cmpt];
+                }
+            }
+
+            break;
+        }
+    }
+
+    return vt;
+}
+
+
+Foam::meshedSurface Foam::PDRobstacle::surface() const
+{
+    meshedSurface surf;
+
+    const PDRobstacle& obs = *this;
+
+    switch (obs.typeId)
+    {
+        case PDRobstacle::CUBOID_1 :
+        case PDRobstacle::CUBOID :
+        {
+            boundBox box(obs.pt, obs.pt + obs.span);
+
+            pointField pts(box.points());
+            faceList fcs(boundBox::faces);
+
+            surf.transfer(pts, fcs);
+
+            break;
+        }
+
+        case PDRobstacle::DIAG_BEAM :
+        {
+            boundBox box(Zero);
+
+            switch (orient)
+            {
+                case vector::X:
+                {
+                    box.min() = vector(0, -0.5*obs.wa, -0.5*obs.wb);
+                    box.max() = vector(obs.len(), 0.5*obs.wa,  0.5*obs.wb);
+                    break;
+                }
+
+                case vector::Y:
+                {
+                    box.min() = vector(-0.5*obs.wb, 0, -0.5*obs.wa);
+                    box.max() = vector(0.5*obs.wb, obs.len(), 0.5*obs.wa);
+                    break;
+                }
+
+                case vector::Z:
+                {
+                    box.min() = vector(-0.5*obs.wa, -0.5*obs.wb, 0);
+                    box.max() = vector(0.5*obs.wa,  0.5*obs.wb, obs.len());
+                    break;
+                }
+            }
+
+            coordinateSystem cs
+            (
+                obs.pt,
+                coordinateRotations::axisAngle
+                (
+                    vector::components(obs.orient),
+                    obs.theta(),
+                    false
+                )
+            );
+
+            pointField pts0(box.points());
+            faceList fcs(boundBox::faces);
+
+            pointField pts(cs.globalPosition(pts0));
+
+            surf.transfer(pts, fcs);
+
+            break;
+        }
+
+        case PDRobstacle::CYLINDER :
+        {
+            // Tessellation 12 looks fairly reasonable
+
+            constexpr int nDiv = 12;
+
+            point org(obs.pt);
+
+            direction e1 = vector::X;
+            direction e2 = vector::Y;
+            direction e3 = vector::Z;
+
+            if (obs.orient == vector::X)
+            {
+                e1 = vector::Y;
+                e2 = vector::Z;
+                e3 = vector::X;
+            }
+            else if (obs.orient == vector::Y)
+            {
+                e1 = vector::Z;
+                e2 = vector::X;
+                e3 = vector::Y;
+            }
+
+            pointField pts(2*nDiv, org);
+            faceList   fcs(2 + nDiv);
+
+            // Origin for back
+            org[e3] += obs.len();
+            SubList<point>(pts, nDiv, nDiv) = org;
+
+            const scalar radius = 0.5*obs.dia();
+
+            for (label i=0; i < nDiv; ++i)
+            {
+                const scalar angle = (i * mathematical::twoPi) / nDiv;
+                const scalar s = radius * sin(angle);
+                const scalar c = radius * cos(angle);
+
+                pts[i][e1] += s;
+                pts[i][e2] += c;
+
+                pts[nDiv+i][e1] += s;
+                pts[nDiv+i][e2] += c;
+            }
+
+            // Side-faces
+            for (label facei=0; facei < nDiv; ++facei)
+            {
+                face& f = fcs[facei];
+                f.resize(4);
+
+                f[0] = facei;
+                f[3] = (facei + 1) % nDiv;
+                f[1] = f[0] + nDiv;
+                f[2] = f[3] + nDiv;
+            }
+
+            {
+                // Front face
+                face& f1 = fcs[nDiv];
+                f1.resize(nDiv);
+
+                f1[0] = 0;
+                for (label pti=1; pti < nDiv; ++pti)
+                {
+                    f1[pti] = nDiv-pti;
+                }
+
+                // Back face
+                labelList& f2 = fcs[nDiv+1];
+                f2 = identity(nDiv, nDiv);
+            }
+
+            surf.transfer(pts, fcs);
+
+            break;
+        }
+
+        case PDRobstacle::RECT_PATCH :
+        {
+            pointField pts(4, obs.span);
+            pts[0] = Zero;
+
+            switch (obs.inlet_dirn)
+            {
+                case -1:
+                case 1:
+                {
+                    for (auto& p : pts)
+                    {
+                        p.x() = 0;
+                    }
+
+                    pts[1].z() = 0;
+                    pts[3].y() = 0;
+                    break;
+                }
+                case -2:
+                case 2:
+                {
+                    for (auto& p : pts)
+                    {
+                        p.y() = 0;
+                    }
+
+                    pts[1].x() = 0;
+                    pts[3].z() = 0;
+                    break;
+                }
+                default:
+                {
+                    for (auto& p : pts)
+                    {
+                        p.z() = 0;
+                    }
+
+                    pts[1].y() = 0;
+                    pts[3].x() = 0;
+                    break;
+                }
+            }
+
+            // pts += obs.pt;
+
+            faceList fcs(one(), face(identity(4)));
+
+            surf.transfer(pts, fcs);
+
+            break;
+        }
+
+        default:
+            break;
+
+//         LOUVRE_BLOWOFF = 5,
+//         WALL_BEAM =  7,
+//         GRATING   =  8,
+//         CIRC_PATCH  = 12,
+//         MESH_PLANE  = 46,
+    }
+
+    return surf;
+}
+
+
+Foam::label Foam::PDRobstacle::addPieces
+(
+    vtk::surfaceWriter& surfWriter,
+    const UList<PDRobstacle>& list,
+    label pieceId
+)
+{
+    for (const PDRobstacle& obs : list)
+    {
+        meshedSurface surf(obs.surface());
+
+        if (!surf.empty())
+        {
+            surfWriter.piece(surf.points(), surf.surfFaces());
+
+            surfWriter.writeGeometry();
+            surfWriter.beginCellData(2);
+            surfWriter.writeUniform("group", label(obs.groupId));
+            surfWriter.writeUniform("type", label(obs.typeId));
+            surfWriter.writeUniform("obstacle", pieceId);
+            ++pieceId;
+        }
+    }
+
+    return pieceId;
+}
+
+
+void Foam::PDRobstacle::generateVtk
+(
+    const fileName& outputDir,
+    const UList<PDRobstacle>& obslist,
+    const UList<PDRobstacle>& cyllist
+)
+{
+    label pieceId = 0;
+
+    meshedSurf::emptySurface dummy;
+
+    vtk::surfaceWriter surfWriter
+    (
+        dummy.points(),
+        dummy.faces(),
+        // vtk::formatType::INLINE_ASCII,
+        (outputDir / "Obstacles"),
+        false  // serial only
+    );
+
+    pieceId = addPieces(surfWriter, obslist, pieceId);
+    pieceId = addPieces(surfWriter, cyllist, pieceId);
+
+    Info<< "Wrote " << pieceId << " obstacles (VTK) to "
+        << outputDir/"Obstacles" << nl;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const InfoProxy<PDRobstacle>& iproxy
+)
+{
+    const PDRobstacle& obs = iproxy.t_;
+
+    switch (obs.typeId)
+    {
+        case PDRobstacle::CUBOID_1 :
+        case PDRobstacle::CUBOID :
+            os  << "box  { point " << obs.pt
+                << "; size " << obs.span
+                << "; }";
+            break;
+
+        case PDRobstacle::CYLINDER :
+            os  << "cyl { point " << obs.pt
+                << "; length " << obs.len() << "; diameter " << obs.dia()
+                << "; direction " << vector::componentNames[obs.orient]
+                << "; }";
+            break;
+
+        case PDRobstacle::DIAG_BEAM :
+            os  << "diag { point " << obs.pt
+                << "; length " << obs.len()
+                << "; width (" << obs.wa << ' ' << obs.wb << ')'
+                << "; angle " << radToDeg(obs.theta())
+                << "; direction " << vector::componentNames[obs.orient]
+                << "; }";
+            break;
+
+        case PDRobstacle::WALL_BEAM :
+            os  << "wallbeam { point " << obs.pt
+                << " size " << obs.span
+                << "; }";
+            break;
+
+        case PDRobstacle::GRATING :
+            os  << "grate { point " << obs.pt
+                << "; size " << obs.span
+                << "; slats " << obs.slat_width
+                << "; }";
+            break;
+
+        case PDRobstacle::LOUVER_BLOWOFF :
+            os  << "louver { point " << obs.pt
+                << "; size " << obs.span
+                << "; pressure " << paToBar(obs.blowoff_press)
+                << "; }";
+            break;
+
+        case PDRobstacle::RECT_PATCH :
+            os  << "patch { " << obs.pt
+                << "; size " << obs.span
+                << "; name " << obs.identifier
+                << "; }";
+            break;
+
+        case PDRobstacle::OLD_INLET :
+        case PDRobstacle::OLD_BLOWOFF :
+        case PDRobstacle::IGNITION :
+            os  << "/* ignored: " << obs.typeId << " */";
+            break;
+
+        default:
+            os  << "/* unknown: " << obs.typeId << " */";
+            break;
+
+    }
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+bool Foam::operator<(const PDRobstacle& a, const PDRobstacle& b)
+{
+    return (a.pt.x() + a.sortBias) < (b.pt.x() + b.sortBias);
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacle.H

@@ -0,0 +1,477 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PDRobstacle
+
+Description
+    Obstacle definitions for PDR
+
+SourceFiles
+    PDRobstacle.C
+    PDRobstacleIO.C
+    PDRobstacleRead.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRobstacle_H
+#define PDRobstacle_H
+
+#include "InfoProxy.H"
+#include "labelPair.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
+#include "boundBox.H"
+#include "DynamicList.H"
+#include "pointField.H"
+#include "volumeType.H"
+#include "memberFunctionSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class boundBox;
+class PDRobstacle;
+
+Istream& operator>>(Istream& is, PDRobstacle& obs);
+Ostream& operator<<(Ostream& os, const InfoProxy<PDRobstacle>& info);
+
+namespace vtk
+{
+    class surfaceWriter;
+}
+
+
+/*---------------------------------------------------------------------------*\
+                         Class PDRobstacle Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRobstacle
+{
+public:
+
+    //- Obstacle types (legacy numbering)
+    enum legacyTypes
+    {
+        NONE =  0,          //!< Placeholder
+        CUBOID_1  =  1,
+        CYLINDER  =  2,
+        LOUVER_BLOWOFF = 5,
+        LOUVRE_BLOWOFF = 5,
+        CUBOID    =  6,
+        WALL_BEAM =  7,
+        GRATING   =  8,
+        OLD_INLET =  9,     //!< ignored (old)
+        OLD_BLOWOFF = 10,   //!< ignored (old)
+        CIRC_PATCH  = 12,
+        RECT_PATCH  = 16,
+        DIAG_BEAM   = 22,
+        IGNITION    = 41,   //!< ignored (old)
+        MESH_PLANE  = 46,
+        IGNORE = 200
+    };
+
+
+    // Static Data Members
+
+        //- The max blowoff pressure [bar]
+        //  Primarily to catch accidental input in Pa or mbar
+        static constexpr int maxBlowoffPressure = 10;
+
+
+    // Data Members
+
+        //- The group-id
+        label groupId;
+
+        //- The obstacle type-id
+        int typeId;
+
+        //- The x/y/z orientation (0,1,2)
+        direction orient;
+
+        //- Bias for position sorting
+        scalar sortBias;
+
+        //- The obstacle location.
+        //  Lower corner for boxes, end-centre for cylinders
+        point pt;
+
+        //- The obstacle dimensions (for boxes)
+        vector span;
+
+        // Accessors for cylinders and diagonal blocks
+
+        inline scalar dia()   const { return span[vector::X]; }
+        inline scalar theta() const { return span[vector::Y]; }
+        inline scalar len()   const { return span[vector::Z]; }
+
+        inline scalar& dia()   { return span[vector::X]; }
+        inline scalar& theta() { return span[vector::Y]; }
+        inline scalar& len()   { return span[vector::Z]; }
+
+        union
+        {
+            scalar wa;
+            scalar slat_width;
+            scalar blowoff_press;
+        };
+        union
+        {
+            scalar wb;
+            scalar blowoff_time;
+        };
+        scalar vbkge;
+        scalar xbkge;
+        scalar ybkge;
+        scalar zbkge;
+
+        union
+        {
+            int blowoff_type;
+            int inlet_dirn;
+        };
+
+        string identifier;
+
+public:
+
+    // Constructors
+
+        //- Construct zero-initialized
+        PDRobstacle();
+
+        //- Read construct as named dictionary
+        explicit PDRobstacle(Istream& is);
+
+
+    // Member Function Selectors
+
+        declareMemberFunctionSelectionTable
+        (
+            void,
+            PDRobstacle,
+            read,
+            dictRead,
+            (
+                PDRobstacle& obs,
+                const dictionary& dict
+            ),
+            (obs, dict)
+        );
+
+
+    // Static Member Functions
+
+        //- Read obstacle files and add to the lists
+        //  \return the total volume
+        static scalar legacyReadFiles
+        (
+            const fileName& obsFileDir,
+            const wordList& obsFileNames,
+            const boundBox& meshBb,
+            // output
+            DynamicList<PDRobstacle>& blocks,
+            DynamicList<PDRobstacle>& cylinders
+        );
+
+        //- Read obstacle files and set the lists
+        //  \return the total volume
+        static scalar readFiles
+        (
+            const fileName& obsFileDir,
+            const wordList& obsFileNames,
+            const boundBox& meshBb,
+            // output
+            DynamicList<PDRobstacle>& blocks,
+            DynamicList<PDRobstacle>& cylinders
+        );
+
+
+    // Member Functions
+
+        //- Read name / dictionary
+        bool read(Istream& is);
+
+        //- Read the 'name' identifier if present
+        void readProperties(const dictionary& dict);
+
+        //- Obstacle position accessors
+        inline scalar x() const { return pt.x(); }
+        inline scalar y() const { return pt.y(); }
+        inline scalar z() const { return pt.z(); }
+        inline scalar& x() { return pt.x(); }
+        inline scalar& y() { return pt.y(); }
+        inline scalar& z() { return pt.z(); }
+
+
+        //- Is obstacle type id cylinder-like?
+        inline static bool isCylinder(const label id);
+
+        //- Is obstacle cylinder-like?
+        inline bool isCylinder() const;
+
+        //- Reset to a zero obstacle
+        void clear();
+
+        //- Scale obstacle dimensions by specified scaling factor
+        //  Zero and negative factors are ignored
+        void scale(const scalar factor);
+
+        //- Volume of the obstacle
+        scalar volume() const;
+
+        //- True if the obstacle is considered to be too small
+        bool tooSmall(const scalar minWidth) const;
+
+        //- Set values from single-line, multi-column format.
+        //  The only input format, but termed 'legacy' since it may
+        //  be replaced in the near future.
+        //  \return false if the scanning failed or if the obstacle type
+        //      is not supported (or no longer supported)
+        bool setFromLegacy
+        (
+            const int groupTypeId,
+            const string& buffer,
+            const label lineNo = -1,
+            const word& inputFile = word::null
+        );
+
+        //- Trim obstacle to ensure it is within the specified bounding box
+        volumeType trim(const boundBox& bb);
+
+        //- Surface (points, faces) representation
+        meshedSurface surface() const;
+
+        //- Add pieces to vtp output
+        static label addPieces
+        (
+            vtk::surfaceWriter& surfWriter,
+            const UList<PDRobstacle>& list,
+            label pieceId = 0
+        );
+
+        //- Generate multi-piece VTK (vtp) file of obstacles
+        static void generateVtk
+        (
+            const fileName& outputDir,
+            const UList<PDRobstacle>& obslist,
+            const UList<PDRobstacle>& cyllist
+        );
+
+
+    // IOstream Operators
+
+        //- Return info proxy.
+        InfoProxy<PDRobstacle> info() const
+        {
+            return *this;
+        }
+
+        friend Istream& operator>>(Istream& is, PDRobstacle& obs);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Global Operators
+
+//- Compare according to x0 position
+bool operator<(const PDRobstacle& a, const PDRobstacle& b);
+
+//- For list output, assert that no obstacles are identical
+inline bool operator!=(const PDRobstacle& a, const PDRobstacle& b)
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace PDRlegacy
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class obstacleGrouping Declaration
+\*---------------------------------------------------------------------------*/
+
+//- Locations for each instance of an obstacle group.
+class obstacleGrouping
+:
+    public DynamicList<point>
+{
+    //- Number of obstacles counted
+    label nObstacle_;
+
+    //- Number of cylinder-like obstacles counted
+    label nCylinder_;
+
+
+public:
+
+    //- Construct null
+    obstacleGrouping()
+    :
+        nObstacle_(0),
+        nCylinder_(0)
+    {}
+
+    //- Construct with one location (instance)
+    explicit obstacleGrouping(const vector& origin)
+    :
+        obstacleGrouping()
+    {
+        append(origin);
+    }
+
+    //- Clear obstacle count and locations
+    void clear()
+    {
+        nObstacle_ = 0;
+        nCylinder_ = 0;
+        DynamicList<point>::clear();
+    }
+
+    //- Increment the number of obstacles
+    void addObstacle()
+    {
+        ++nObstacle_;
+    }
+
+    //- Increment the number of cylinder-like obstacles
+    void addCylinder()
+    {
+        ++nCylinder_;
+    }
+
+    //- The number of obstacles
+    label nObstacle() const
+    {
+        return nObstacle_;
+    }
+
+    //- The number of cylinder-like obstacles
+    label nCylinder() const
+    {
+        return nCylinder_;
+    }
+
+    //- The number of locations x number of obstacles
+    label nTotalObstacle() const
+    {
+        return size() * nObstacle_;
+    }
+
+    //- The number of locations x number of cylinder-like obstacles
+    label nTotalCylinder() const
+    {
+        return size() * nCylinder_;
+    }
+
+    //- The number of locations x number of obstacles
+    label nTotal() const
+    {
+        return size() * (nObstacle_ + nCylinder_);
+    }
+
+    //- Add a location
+    using DynamicList<point>::append;
+
+    //- Add a location
+    void append(const scalar x, const scalar y, const scalar z)
+    {
+        append(point(x, y, z));
+    }
+};
+
+
+// Service Functions
+
+//- Read obstacle files, do counting only.
+//  \return nObstacle, nCylinder read
+labelPair readObstacleFiles
+(
+    const fileName& obsFileDir,
+    const wordList& obsFileNames,
+    Map<obstacleGrouping>& groups
+);
+
+
+//- Read obstacle files and add to the lists
+//  \return the total volume
+scalar readObstacleFiles
+(
+    const fileName& obsFileDir,
+    const wordList& obsFileNames,
+    const Map<obstacleGrouping>& groups,
+    const boundBox& meshBb,
+    // output
+    DynamicList<PDRobstacle>& blocks,
+    DynamicList<PDRobstacle>& cylinders
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace PDRlegacy
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Global Operators
+
+//- Locations for each instance of an obstacle group.
+inline Ostream& operator<<
+(
+    Ostream& os,
+    const PDRlegacy::obstacleGrouping& group
+)
+{
+    os  << token::BEGIN_LIST
+        << group.size() << token::SPACE
+        << group.nObstacle() << token::SPACE
+        << group.nCylinder() << token::END_LIST;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "PDRobstacleI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleI.H

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::PDRobstacle::isCylinder(const label id)
+{
+    return
+    (
+        id == PDRobstacle::CYLINDER
+     || id == PDRobstacle::DIAG_BEAM
+    );
+}
+
+
+inline bool Foam::PDRobstacle::isCylinder() const
+{
+    return isCylinder(typeId);
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleIO.C

@@ -0,0 +1,353 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRsetFields.H"
+#include "PDRobstacle.H"
+#include "volumeType.H"
+
+using namespace Foam;
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::PDRobstacle::read(Istream& is)
+{
+    this->clear();
+
+    const word obsType(is);
+    const dictionary dict(is);
+
+    const auto mfIter = readdictReadMemberFunctionTablePtr_->cfind(obsType);
+
+    if (!mfIter.good())
+    {
+        FatalIOErrorInFunction(is)
+            << "Unknown obstacle type: " << obsType << nl
+            << "Valid types:" << nl
+            << readdictReadMemberFunctionTablePtr_->sortedToc() << nl
+            << exit(FatalIOError);
+    }
+
+    mfIter()(*this, dict);
+
+    return true;
+}
+
+
+Foam::scalar Foam::PDRobstacle::readFiles
+(
+    const fileName& obsFileDir,
+    const wordList& obsFileNames,
+    const boundBox& meshBb,
+
+    DynamicList<PDRobstacle>& blocks,
+    DynamicList<PDRobstacle>& cylinders
+)
+{
+    blocks.clear();
+    cylinders.clear();
+
+    scalar totVolume = 0;
+    label nOutside = 0;
+    label nProtruding = 0;
+
+    scalar shift = pars.obs_expand;
+
+    if (!obsFileNames.empty())
+    {
+        Info<< "Reading obstacle files" << nl;
+    }
+
+    label maxGroup = -1;
+
+    for (const word& inputFile : obsFileNames)
+    {
+        Info<< "    file: " << inputFile << nl;
+
+        fileName path = (obsFileDir / inputFile);
+
+        IFstream is(path);
+        dictionary inputDict(is);
+
+        const scalar scaleFactor = inputDict.getOrDefault<scalar>("scale", 0);
+
+        const label verbose = inputDict.getOrDefault<label>("verbose", 0);
+
+        for (const entry& dEntry : inputDict)
+        {
+            if (!dEntry.isDict())
+            {
+                // ignore non-dictionary entry
+                continue;
+            }
+
+            const dictionary& dict = dEntry.dict();
+
+            if (!dict.getOrDefault("enabled", true))
+            {
+                continue;
+            }
+
+            label obsGroupId = 0;
+            if (dict.readIfPresent("groupId", obsGroupId))
+            {
+                maxGroup = max(maxGroup, obsGroupId);
+            }
+            else
+            {
+                obsGroupId = ++maxGroup;
+            }
+
+
+            pointField pts;
+            dict.readIfPresent("locations", pts);
+            if (pts.empty())
+            {
+                pts.resize(1, Zero);
+            }
+
+            List<PDRobstacle> obsInput;
+            dict.readEntry("obstacles", obsInput);
+
+            label nCyl = 0; // The number of cylinders vs blocks
+
+            for (PDRobstacle& obs : obsInput)
+            {
+                obs.groupId = obsGroupId;
+                obs.scale(scaleFactor);
+
+                if (obs.isCylinder())
+                {
+                    ++nCyl;
+                }
+            }
+
+            const label nBlock = (obsInput.size() - nCyl);
+
+            blocks.reserve(blocks.size() + nBlock*pts.size());
+            cylinders.reserve(cylinders.size() + nCyl*pts.size());
+
+            if (verbose)
+            {
+                Info<< "Read " << obsInput.size() << " obstacles ("
+                    << nCyl << " cylinders) with "
+                    << pts.size() << " locations" << nl;
+
+                if (verbose > 1)
+                {
+                    Info<< "locations " << pts << nl
+                        << "obstacles " << obsInput << nl;
+                }
+            }
+
+            for (const PDRobstacle& scanObs : obsInput)
+            {
+                // Reject anything below minimum width
+                if (scanObs.tooSmall(pars.min_width))
+                {
+                    continue;
+                }
+
+                for (const point& origin : pts)
+                {
+                    // A different (very small) shift for each obstacle
+                    // so that faces cannot be coincident
+
+                    shift += floatSMALL;
+                    const scalar shift2 = shift * 2.0;
+
+
+                switch (scanObs.typeId)
+                {
+                    case PDRobstacle::CYLINDER:
+                    {
+                        // Make a copy
+                        PDRobstacle obs(scanObs);
+
+                        // Offset for the group position
+                        obs.pt += origin;
+
+                        // Shift the end outwards so, if exactly on
+                        // cell boundary, now overlap cell.
+                        // So included in Aw.
+                        obs.pt -= point::uniform(shift);
+                        obs.len() += shift2;
+                        obs.dia() -= floatSMALL;
+
+
+                        // Trim against the mesh bounds
+                        // - ignore if it doesn't overlap
+                        const volumeType vt = obs.trim(meshBb);
+
+                        switch (vt)
+                        {
+                            case volumeType::OUTSIDE:
+                                ++nOutside;
+                                continue; // Can ignore the rest
+                                break;
+
+                            case volumeType::MIXED:
+                                ++nProtruding;
+                                break;
+
+                            default:
+                                break;
+                        }
+
+                        // Later for position sorting
+                        switch (obs.orient)
+                        {
+                            case vector::X:
+                                obs.sortBias = obs.len();
+                                break;
+                            case vector::Y:
+                                obs.sortBias = 0.5*obs.dia();
+                                break;
+                            case vector::Z:
+                                obs.sortBias = 0.5*obs.dia();
+                                break;
+                        }
+
+                        totVolume += obs.volume();
+                        cylinders.append(obs);
+
+                        break;
+                    }
+
+                    case PDRobstacle::DIAG_BEAM:
+                    {
+                        // Make a copy
+                        PDRobstacle obs(scanObs);
+
+                        // Offset for the group position
+                        obs.pt += origin;
+
+                        // Shift the end outwards so, if exactly on
+                        // cell boundary, now overlap cell.
+                        // So included in Aw.
+                        obs.pt -= point::uniform(shift);
+                        obs.len() += shift2;
+                        obs.wa += shift2;
+                        obs.wb += shift2;
+
+                        totVolume += obs.volume();
+                        cylinders.append(obs);
+
+                        break;
+                    }
+
+                    case PDRobstacle::CUBOID_1:
+                    case PDRobstacle::LOUVRE_BLOWOFF:
+                    case PDRobstacle::CUBOID:
+                    case PDRobstacle::WALL_BEAM:
+                    case PDRobstacle::GRATING:
+                    case PDRobstacle::RECT_PATCH:
+                    {
+                        // Make a copy
+                        PDRobstacle obs(scanObs);
+
+                        // Offset for the position of the group
+                        obs.pt += origin;
+
+                        if (obs.typeId == PDRobstacle::GRATING)
+                        {
+                            if (obs.slat_width <= 0)
+                            {
+                                obs.slat_width = pars.def_grating_slat_w;
+                            }
+                        }
+
+                        // Shift the end outwards so, if exactly on
+                        // cell boundary, now overlap cell.
+                        // So included in Aw.
+                        obs.pt -= point::uniform(shift);
+                        obs.span += point::uniform(shift2);
+
+
+                        // Trim against the mesh bounds
+                        // - ignore if it doesn't overlap
+                        const volumeType vt = obs.trim(meshBb);
+
+                        switch (vt)
+                        {
+                            case volumeType::OUTSIDE:
+                                ++nOutside;
+                                continue; // Can ignore the rest
+                                break;
+
+                            case volumeType::MIXED:
+                                ++nProtruding;
+                                break;
+
+                            default:
+                                break;
+                        }
+
+                        totVolume += obs.volume();
+
+                        blocks.append(obs);
+
+                        break;
+                    }
+                }
+            }
+            }
+
+            // Info<< "Cylinders: " << cylinders << nl;
+        }
+
+        if (nOutside || nProtruding)
+        {
+            Info<< "Warning: " << nOutside << " obstacles outside domain, "
+                << nProtruding << " obstacles partly outside domain" << nl;
+        }
+    }
+
+    // #ifdef FULLDEBUG
+    // Info<< blocks << nl << cylinders << nl;
+    // #endif
+
+
+    Info<< "Number of obstacles: "
+        << (blocks.size() + cylinders.size()) << " ("
+        << cylinders.size() << " cylinders)" << nl;
+
+    return totVolume;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, PDRobstacle& obs)
+{
+    obs.read(is);
+
+    return is;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleLegacyIO.C

@@ -0,0 +1,475 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRobstacle.H"
+#include "vector.H"
+#include "doubleVector.H"
+#include "stringOps.H"
+#include "unitConversion.H"
+#include <cmath>
+
+#define ReportLineInfo(line, file)                                            \
+    if (line >= 0 && !file.empty())                                           \
+    {                                                                         \
+        Info<< " Line " << line << " of file '" << file << '\'';              \
+    }                                                                         \
+    Info<< nl;
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::PDRobstacle::setFromLegacy
+(
+    const int groupTypeId,
+    const string& buffer,
+    const label lineNo,
+    const word& inputFile
+)
+{
+    // Handling the optional identifier string can be a pain.
+    // Generally can only say that it exists if there are 15 or
+    // more columns.
+    //
+    // Cylinder has 8 normal entries
+    // Cuboid, diagonal beam etc have 14 normal entries
+    // However, reject anything that looks like a slipped numeric
+
+    double dummy1;
+
+    string in_ident;
+
+    const auto columns = stringOps::splitSpace(buffer);
+
+    for (std::size_t coli = 14; coli < columns.size(); ++coli)
+    {
+        // See if it can parse into a numerical value
+        if (!readDouble(columns[coli].str(), dummy1))
+        {
+            // Not a numeric value. This must be our identifier
+            in_ident = buffer.substr(columns[coli].first - buffer.begin());
+
+            #ifdef FULLDEBUG
+            Info<< "Identifier: " << in_ident << nl;
+            #endif
+            break;
+        }
+    }
+
+    // Strip off group number
+    groupId = groupTypeId / 100;
+    typeId  = groupTypeId % 100;
+
+    // This is a safe value
+    orient = vector::X;
+
+    switch (typeId)
+    {
+        case PDRobstacle::CYLINDER:
+        {
+            // 8 Tokens
+            // "%d %lf %lf %lf %lf %lf %d %lf"
+            // USP 13/8/14  Read vbkge in case a negative cyl to punch a circular hole
+
+            int in_typeId;
+            double in_x, in_y, in_z;
+            double in_len, in_dia;
+            int in_orient;
+            double in_poro;
+
+            int nread =
+                sscanf
+                (
+                    buffer.c_str(),
+                    "%d %lf %lf %lf %lf %lf %d %lf",
+                    &in_typeId, &in_x, &in_y, &in_z,
+                    &in_len, &in_dia, &in_orient,
+                    &in_poro
+                );
+
+            if (nread < 8)
+            {
+                Info<< "Expected 8 items, but read in " << nread;
+                ReportLineInfo(lineNo, inputFile);
+            }
+
+            identifier = in_ident;
+            pt = point(in_x, in_y, in_z);
+
+            len() = in_len;
+            dia() = in_dia;
+
+            orient = vector::X;  // Check again later
+
+            // Read porosity. Convert to blockage.
+            vbkge = 1.0 - in_poro;
+
+            // Orientation (1,2,3) on input -> (0,1,2)
+            // - sortBias for later position sorting
+            switch (in_orient)
+            {
+                case 1:
+                    orient = vector::X;
+                    sortBias = len();
+                    break;
+                case 2:
+                    orient = vector::Y;
+                    sortBias = 0.5*dia();
+                    break;
+                case 3:
+                    orient = vector::Z;
+                    sortBias = 0.5*dia();
+                    break;
+                default:
+                    sortBias = len();
+                    Info<< "Unexpected orientation " << in_orient;
+                    ReportLineInfo(lineNo, inputFile);
+                    break;
+            }
+        }
+        break;
+
+        case PDRobstacle::DIAG_BEAM:
+        {
+            // A diagonal block
+
+            // 14 columns + identifier
+            // "%d %lf %lf %lf %lf %lf %d %lf %lf %lf %lf %lf %d %lf %s"
+            // vbkge (porosity at this stage) should be 0. Not used (yet)
+
+            int in_typeId;
+            double in_x, in_y, in_z;
+            double in_len, in_theta;
+            int in_orient;
+            double in_wa, in_wb, in_poro;
+            double col_11, col_12, col_14;
+            int col_13;
+
+            int nread =
+                sscanf
+                (
+                    buffer.c_str(),
+                    "%d %lf %lf %lf %lf %lf %d %lf %lf %lf %lf %lf %d %lf",
+                    &in_typeId, &in_x, &in_y, &in_z,
+                    &in_len, &in_theta, &in_orient,
+                    &in_wa, &in_wb, &in_poro,
+                    &col_11, &col_12, &col_13, &col_14
+                );
+
+            if (nread < 14)
+            {
+                Info<< "Expected min 10 items, but read in " << nread;
+                ReportLineInfo(lineNo, inputFile);
+            }
+
+            identifier = in_ident;
+            pt = point(in_x, in_y, in_z);
+
+            len() = in_len;
+            dia() = 0;
+            theta() = 0; // Fix later on
+
+            orient = vector::X;  // Check again later
+
+            wa  = in_wa;
+            wb  = in_wb;
+
+            // Degrees on input, limit to range [0, PI]
+            while (in_theta > 180) in_theta -= 180;
+            while (in_theta < 0) in_theta += 180;
+
+            // Swap axes when theta > PI/2
+            // For 89-90 degrees it becomes -ve, which is picked up
+            // in next section
+            if (in_theta > 89)
+            {
+                in_theta -= 90;
+                // Swap wa <-> wb
+                wa = in_wb;
+                wb = in_wa;
+            }
+
+            theta() = degToRad(in_theta);
+
+            // Orientation (1,2,3) on input -> (0,1,2)
+            // - sortBias for later position sorting
+            switch (in_orient)
+            {
+                case 1:
+                    orient = vector::X;
+                    sortBias = len();
+                    break;
+
+                case 2:
+                    orient = vector::Y;
+                    sortBias = 0.5*(wa * sin(theta()) + wb * cos(theta()));
+                    break;
+
+                case 3:
+                    orient = vector::Z;
+                    sortBias = 0.5*(wa * cos(theta()) + wb * sin(theta()));
+                    break;
+
+                default:
+                    sortBias = len();
+                    Info<< "Unexpected orientation " << in_orient;
+                    ReportLineInfo(lineNo, inputFile);
+                    break;
+            }
+
+
+            // If very nearly aligned with axis, turn it into normal block,
+            // to avoid 1/tan(theta) blowing up
+            if (in_theta < 1)
+            {
+                switch (orient)
+                {
+                    case vector::X:
+                        span = vector(len(), wa, wb);
+                        // Was end center, now lower corner
+                        pt.y() = pt.y() - 0.5 * span.y();
+                        pt.z() = pt.z() - 0.5 * span.z();
+                        break;
+
+                    case vector::Y:
+                        span = vector(wb, len(), wa);
+                        // Was end center, now lower corner
+                        pt.z() = pt.z() - 0.5 * span.z();
+                        pt.x() = pt.x() - 0.5 * span.x();
+                        break;
+
+                    case vector::Z:
+                        span = vector(wa, wb, len());
+                        // Was end center, now lower corner
+                        pt.x() = pt.x() - 0.5 * span.x();
+                        pt.y() = pt.y() - 0.5 * span.y();
+                        break;
+                }
+
+                typeId = PDRobstacle::CUBOID;
+                sortBias = 0;
+                xbkge = ybkge = zbkge = vbkge = 1.0;
+                blowoff_type = 0;
+
+                Info<< "... changed to type cuboid" << nl;
+                break;
+            }
+        }
+        break;
+
+        case PDRobstacle::CUBOID_1:       // Cuboid "Type 1"
+        case PDRobstacle::LOUVRE_BLOWOFF: // Louvred wall or blow-off panel
+        case PDRobstacle::CUBOID:         // Cuboid
+        case PDRobstacle::WALL_BEAM:      // Beam against wall (treated here as normal cuboid)
+        case PDRobstacle::GRATING:        // Grating
+        case PDRobstacle::RECT_PATCH:     // Inlet, outlet, other b.c. (rectangular)
+        {
+            // 14 columns + identifier
+            // "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %lf"
+
+            int in_typeId;
+            double in_x, in_y, in_z;
+            double in_delx, in_dely, in_delz;
+            double in_poro, in_porox, in_poroy, in_poroz;
+            double col_12;
+            int col_13;
+            double in_blowoff_time = 0;
+
+            int nread =
+                sscanf
+                (
+                    buffer.c_str(),
+                    "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %lf",
+                    &in_typeId, &in_x, &in_y, &in_z,
+                    &in_delx, &in_dely, &in_delz,
+                    &in_poro, &in_porox, &in_poroy, &in_poroz,
+                    &col_12, &col_13, &in_blowoff_time
+                );
+
+            blowoff_time = scalar(in_blowoff_time);
+
+            if (nread < 14)
+            {
+                Info<< "Expected 14 items, but read in " << nread;
+                ReportLineInfo(lineNo, inputFile);
+            }
+
+            identifier = in_ident;
+            pt = point(in_x, in_y, in_z);
+
+            span = vector(in_delx, in_dely, in_delz);
+
+            // Read porosity. Convert to blockage.
+            vbkge = 1.0 - in_poro;
+            xbkge = 1.0 - in_porox;
+            ybkge = 1.0 - in_poroy;
+            zbkge = 1.0 - in_poroz;
+
+            if
+            (
+                typeId == PDRobstacle::CUBOID_1
+             || typeId == PDRobstacle::WALL_BEAM
+             || typeId == PDRobstacle::RECT_PATCH
+            )
+            {
+                // Check for invalid input
+
+                if (vbkge != 1.0 || xbkge != 1.0 || ybkge != 1.0 || zbkge != 1.0)
+                {
+                    Info<< "Type " << typeId << " is porous (setting to blockage).";
+                    ReportLineInfo(lineNo, inputFile);
+
+                    vbkge = 1;
+                    xbkge = 1;
+                    ybkge = 1;
+                    zbkge = 1;
+                }
+                if (typeId == PDRobstacle::RECT_PATCH)
+                {
+                    // Correct interpretation of column 13
+                    inlet_dirn = col_13;
+
+                    if (identifier.empty())
+                    {
+                        FatalErrorInFunction
+                            << "RECT_PATCH without a patch name"
+                            << exit(FatalError);
+                    }
+                }
+            }
+            else if (typeId == PDRobstacle::CUBOID)
+            {
+            }
+            else
+            {
+                if (!equal(cmptProduct(span), 0))
+                {
+                    Info<< "Type " << typeId << " has non-zero thickness.";
+                    ReportLineInfo(lineNo, inputFile);
+                }
+            }
+
+            if (typeId == PDRobstacle::LOUVRE_BLOWOFF)
+            {
+                // Blowoff panel
+                blowoff_press = barToPa(col_12);
+                blowoff_type = col_13;
+
+                if (blowoff_type == 1)
+                {
+                    Info<< "Type " << typeId
+                        << ": blowoff-type 1 not yet implemented.";
+                    ReportLineInfo(lineNo, inputFile);
+
+                    if (blowoff_time != 0)
+                    {
+                        Info<< "Type " << typeId << " has blowoff time set,"
+                            << " not set to blow off cell-by-cell";
+                        ReportLineInfo(lineNo, inputFile);
+                    }
+                }
+                else
+                {
+                    if
+                    (
+                        (blowoff_type == 1 || blowoff_type == 2)
+                     && (col_12 > 0)
+                    )
+                    {
+                        if (col_12 > maxBlowoffPressure)
+                        {
+                            Info<< "Blowoff pressure (" << col_12
+                                << ") too high for blowoff type "
+                                << blowoff_type;
+                            ReportLineInfo(lineNo, inputFile);
+                        }
+                    }
+                    else
+                    {
+                        Info<< "Problem with blowoff parameters";
+                        ReportLineInfo(lineNo, inputFile);
+                        Info<< "Col12 " << col_12
+                            << " Blowoff type " << blowoff_type
+                            << ", blowoff pressure " << blowoff_press << nl;
+                    }
+                }
+            }
+            else if (typeId == PDRobstacle::WALL_BEAM)
+            {
+                // WALL_BEAM against walls only contribute half to drag
+                // if ((col_12 == 1) || (col_12 == -1)) { against_wall_fac = 0.5; }
+            }
+            else if (typeId == PDRobstacle::GRATING)
+            {
+                if (col_12 > 0)
+                {
+                    slat_width = col_12;
+                }
+                else
+                {
+                    slat_width = 0;
+                }
+
+                // Set orientation
+                if (equal(span.x(), 0))
+                {
+                    orient = vector::X;
+                }
+                else if (equal(span.y(), 0))
+                {
+                    orient = vector::Y;
+                }
+                else
+                {
+                    orient = vector::Z;
+                }
+            }
+        }
+        break;
+
+        case 0:  // Group location
+        case PDRobstacle::OLD_INLET: // Ventilation source only
+            return false;
+            break;
+
+        case PDRobstacle::IGNITION:  // Ignition (now ignored. 2019-04)
+            Info<< "Ignition cell type ignored";
+            ReportLineInfo(lineNo, inputFile);
+            return false;
+            break;
+
+        default:
+            Info<< "Unexpected type " << typeId;
+            ReportLineInfo(lineNo, inputFile);
+            return false;
+            break;
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleLegacyRead.C

@@ -0,0 +1,567 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 Shell Research Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRsetFields.H"
+#include "PDRobstacle.H"
+#include "volumeType.H"
+
+using namespace Foam;
+using namespace Foam::constant;
+
+#undef USE_ZERO_INSTANCE_GROUPS
+// #define USE_ZERO_INSTANCE_GROUPS
+
+
+// Counting
+Foam::labelPair Foam::PDRlegacy::readObstacleFiles
+(
+    const fileName& obsFileDir,
+    const wordList& obsFileNames,
+    Map<obstacleGrouping>& groups
+)
+{
+    // Default group with single instance and position (0,0,0)
+    groups(0).clear();
+    groups(0).append(point::zero);
+
+    string buffer;
+
+    if (!obsFileNames.empty())
+    {
+        Info<< "Counting groups in obstacle files" << nl;
+    }
+    for (const word& inputFile : obsFileNames)
+    {
+        Info<< "    file: " << inputFile << nl;
+
+        fileName path = (obsFileDir / inputFile);
+
+        IFstream is(path);
+
+        while (is.good())
+        {
+            // Process each line of obstacle files
+            is.getLine(buffer);
+
+            const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");
+
+            if
+            (
+                firstCh == std::string::npos
+             || buffer[firstCh] == '#'
+            )
+            {
+                // Empty line or comment line
+                continue;
+            }
+
+            int typeId;
+            double x, y, z; // Double (not scalar) to match sscanf spec
+
+            if
+            (
+                sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z)<4
+             || typeId == 0
+             || typeId == PDRobstacle::MESH_PLANE
+            )
+            {
+                continue;
+            }
+
+            x *= pars.scale;
+            y *= pars.scale;
+            z *= pars.scale;
+
+            const label groupId = typeId / 100;
+            typeId %= 100;
+
+            if (typeId == PDRobstacle::OLD_INLET)
+            {
+                Info<< "Ignored old-inlet type" << nl;
+                continue;
+            }
+            else if (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
+            {
+                Info<< "Ignored grating" << nl;
+                continue;
+            }
+
+            if (typeId == 0)
+            {
+                // Defining a group location
+                groups(groupId).append(x, y, z);
+            }
+            else if (PDRobstacle::isCylinder(typeId))
+            {
+                // Increment cylinder count for the group
+                groups(groupId).addCylinder();
+            }
+            else
+            {
+                // Increment obstacle count for the group
+                groups(groupId).addObstacle();
+            }
+        }
+    }
+
+
+    label nTotalObs = 0;
+    label nTotalCyl = 0;
+
+    label nMissedObs = 0;
+    label nMissedCyl = 0;
+
+    forAllConstIters(groups, iter)
+    {
+        const auto& group = iter.val();
+
+        nTotalObs += group.nTotalObstacle();
+        nTotalCyl += group.nTotalCylinder();
+
+        if (group.empty())
+        {
+            nMissedObs += group.nObstacle();
+            nMissedCyl += group.nCylinder();
+        }
+    }
+
+    for (const label groupId : groups.sortedToc())
+    {
+        const auto& group = groups[groupId];
+
+        if (groupId)
+        {
+            if (group.size())
+            {
+                Info<< "Found " << group.size()
+                    << " instances of group " << groupId << " ("
+                    << group.nObstacle() << " obstacles "
+                    << group.nCylinder() << " cylinders)"
+                    << nl;
+            }
+        }
+        else
+        {
+            // The group 0 is for ungrouped obstacles
+            Info<< "Found "
+                << group.nObstacle() << " obstacles "
+                << group.nCylinder() << " cylinders not in groups" << nl;
+        }
+    }
+
+    Info<< "Number of obstacles: "
+        << (nTotalObs + nTotalCyl) << " ("
+        << nTotalCyl << " cylinders)" << nl;
+
+    if (nMissedObs + nMissedCyl)
+    {
+        #ifdef USE_ZERO_INSTANCE_GROUPS
+
+        nTotalObs += nMissedObs;
+        nTotalCyl += nMissedCyl;
+        Info<< "Adding " << (nMissedObs + nMissedCyl)
+            << " obstacles in groups without instances to default group" << nl;
+
+        #else
+
+        Warning
+            << nl << "Found " << (nMissedObs + nMissedCyl)
+            << " obstacles in groups without instances" << nl << nl;
+
+        if (pars.debugLevel > 1)
+        {
+            for (const label groupId : groups.sortedToc())
+            {
+                const auto& group = groups[groupId];
+
+                if
+                (
+                    groupId && group.empty()
+                 && (group.nObstacle() || group.nCylinder())
+                )
+                {
+                    Info<< "    Group " << groupId << " ("
+                        << group.nObstacle() << " obstacles "
+                        << group.nCylinder() << " cylinders)"
+                        << nl;
+                }
+            }
+        }
+        #endif
+    }
+
+    return labelPair(nTotalObs, nTotalCyl);
+}
+
+
+Foam::scalar Foam::PDRlegacy::readObstacleFiles
+(
+    const fileName& obsFileDir,
+    const wordList& obsFileNames,
+    const Map<obstacleGrouping>& groups,
+    const boundBox& meshBb,
+
+    DynamicList<PDRobstacle>& blocks,
+    DynamicList<PDRobstacle>& cylinders
+)
+{
+    // Catch programming errors
+    if (!groups.found(0))
+    {
+        FatalErrorInFunction
+            << "No default group 0 defined!" << nl
+            << exit(FatalError);
+    }
+
+    scalar totVolume = 0;
+    label nOutside = 0;
+    label nProtruding = 0;
+
+    scalar shift = pars.obs_expand;
+
+    string buffer;
+
+    if (!obsFileNames.empty())
+    {
+        Info<< "Reading obstacle files" << nl;
+    }
+
+    for (const word& inputFile : obsFileNames)
+    {
+        Info<< "    file: " << inputFile << nl;
+
+        fileName path = (obsFileDir / inputFile);
+
+        IFstream is(path);
+
+        label lineNo = 0;
+        while (is.good())
+        {
+            // Process each line of obstacle files
+            ++lineNo;
+            is.getLine(buffer);
+
+            const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");
+
+            if
+            (
+                firstCh == std::string::npos
+             || buffer[firstCh] == '#'
+            )
+            {
+                // Empty line or comment line
+                continue;
+            }
+
+            // Quick reject
+
+            int typeId;     // Int (not label) to match sscanf spec
+            double x, y, z; // Double (not scalar) to match sscanf spec
+
+            if
+            (
+                sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z) < 4
+             || typeId == 0
+             || typeId == PDRobstacle::MESH_PLANE
+            )
+            {
+                continue;
+            }
+
+            int groupId = typeId / 100;
+            typeId %= 100;
+
+            if
+            (
+                typeId == PDRobstacle::OLD_INLET
+             || (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
+            )
+            {
+                // Silent - already warned during counting
+                continue;
+            }
+
+            if (typeId == 0)
+            {
+                // Group location - not an obstacle
+                continue;
+            }
+
+            if (!groups.found(groupId))
+            {
+                // Catch programming errors.
+                // - group should be there after the previous read
+                Warning
+                    << "Encountered undefined group: " << groupId << nl;
+                continue;
+            }
+
+            #ifdef USE_ZERO_INSTANCE_GROUPS
+            const obstacleGrouping& group =
+            (
+                groups[groups[groupId].size() ? groupId : 0]
+            );
+            #else
+            const obstacleGrouping& group = groups[groupId];
+            #endif
+
+            // Add the obstacle to the list with different position
+            // offsets according to its group. Non-group obstacles
+            // are treated as group 0, which has a single instance
+            // with position (0,0,0) and are added only once.
+
+            PDRobstacle scanObs;
+
+            if
+            (
+                !scanObs.setFromLegacy
+                (
+                    (groupId * 100) + typeId,
+                    buffer,
+                    lineNo,
+                    inputFile
+                )
+            )
+            {
+                continue;
+            }
+
+            scanObs.scale(pars.scale);
+
+            // Ignore anything below minimum width
+            if (scanObs.tooSmall(pars.min_width))
+            {
+                continue;
+            }
+
+
+            for (const point& origin : group)
+            {
+                // A different (very small) shift for each obstacle
+                // so that faces cannot be coincident
+                shift += floatSMALL;
+                const scalar shift2 = shift * 2.0;
+
+                switch (typeId)
+                {
+                    case PDRobstacle::CYLINDER:
+                    {
+                        // Make a copy
+                        PDRobstacle obs(scanObs);
+
+                        // Offset for the position of the group
+                        obs.pt += origin;
+
+                        // Shift the end outwards so, if exactly on
+                        // cell boundary, now overlap cell.
+                        // So included in Aw.
+                        obs.pt -= point::uniform(shift);
+                        obs.len() += shift2;
+                        obs.dia() -= floatSMALL;
+
+
+                        // Trim against the mesh bounds
+                        // - ignore if it doesn't overlap
+                        const volumeType vt = obs.trim(meshBb);
+
+                        switch (vt)
+                        {
+                            case volumeType::OUTSIDE:
+                                ++nOutside;
+                                continue; // Can ignore the rest
+                                break;
+
+                            case volumeType::MIXED:
+                                ++nProtruding;
+                                break;
+
+                            default:
+                                break;
+                        }
+
+                        // Later for position sorting
+                        switch (obs.orient)
+                        {
+                            case vector::X:
+                                obs.sortBias = obs.len();
+                                break;
+                            case vector::Y:
+                                obs.sortBias = 0.5*obs.dia();
+                                break;
+                            case vector::Z:
+                                obs.sortBias = 0.5*obs.dia();
+                                break;
+                        }
+
+                        totVolume += obs.volume();
+                        cylinders.append(obs);
+
+                        break;
+                    }
+
+                    case PDRobstacle::DIAG_BEAM:
+                    {
+                        // Make a copy
+                        PDRobstacle obs(scanObs);
+
+                        // Offset for the position of the group
+                        obs.pt += origin;
+
+                        // Shift the end outwards so, if exactly on
+                        // cell boundary, now overlap cell.
+                        // So included in Aw.
+                        obs.pt -= point::uniform(shift);
+                        obs.len() += shift2;
+                        obs.wa += shift2;
+                        obs.wb += shift2;
+
+                        totVolume += obs.volume();
+                        cylinders.append(obs);
+
+                        break;
+                    }
+
+                    case PDRobstacle::CUBOID_1:       // Cuboid "Type 1"
+                    case PDRobstacle::LOUVRE_BLOWOFF: // Louvred wall or blow-off panel
+                    case PDRobstacle::CUBOID:         // Cuboid
+                    case PDRobstacle::WALL_BEAM:      // Beam against wall (treated here as normal cuboid)
+                    case PDRobstacle::GRATING:        // Grating
+                    case PDRobstacle::RECT_PATCH:     // Inlet, outlet or ather b.c. (rectangular)
+                    {
+                        // Make a copy
+                        PDRobstacle obs(scanObs);
+
+                        // Offset for the position of the group
+                        obs.pt += origin;
+
+                        if (typeId == PDRobstacle::GRATING)
+                        {
+                            if (obs.slat_width <= 0)
+                            {
+                                obs.slat_width = pars.def_grating_slat_w;
+                            }
+                        }
+
+                        // Shift the end outwards so, if exactly on
+                        // cell boundary, now overlap cell.
+                        // So included in Aw.
+                        obs.pt -= point::uniform(shift);
+                        obs.span += point::uniform(shift2);
+
+
+                        // Trim against the mesh bounds
+                        // - ignore if it doesn't overlap
+                        const volumeType vt = obs.trim(meshBb);
+
+                        switch (vt)
+                        {
+                            case volumeType::OUTSIDE:
+                                ++nOutside;
+                                continue; // Can ignore the rest
+                                break;
+
+                            case volumeType::MIXED:
+                                ++nProtruding;
+                                break;
+
+                            default:
+                                break;
+                        }
+
+                        totVolume += obs.volume();
+
+                        blocks.append(obs);
+
+                        break;
+                    }
+                }
+            }
+        }
+
+        if (nOutside || nProtruding)
+        {
+            Info<< "Warning: " << nOutside << " obstacles outside domain, "
+                << nProtruding << " obstacles partly outside domain" << nl;
+        }
+    }
+
+
+    // #ifdef FULLDEBUG
+    // Info<< blocks << nl << cylinders << nl;
+    // #endif
+
+    return totVolume;
+}
+
+
+Foam::scalar Foam::PDRobstacle::legacyReadFiles
+(
+    const fileName& obsFileDir,
+    const wordList& obsFileNames,
+    const boundBox& meshBb,
+    DynamicList<PDRobstacle>& blocks,
+    DynamicList<PDRobstacle>& cylinders
+)
+{
+    // Still just with legacy reading
+
+    // Count the obstacles and get the group locations
+    Map<PDRlegacy::obstacleGrouping> groups;
+
+    const labelPair obsCounts =
+        PDRlegacy::readObstacleFiles(obsFileDir, obsFileNames, groups);
+
+    const label nObstacle = obsCounts.first();
+    const label nCylinder = obsCounts.second();
+
+    // Info<< "grouping: " << groups << endl;
+
+    if (!nObstacle && !nCylinder)
+    {
+        FatalErrorInFunction
+            << "No obstacles in domain" << nl
+            << exit(FatalError);
+    }
+
+    blocks.clear();
+    blocks.reserve(4 * max(4, nObstacle));
+
+    cylinders.clear();
+    cylinders.reserve(4 * max(4, nCylinder));
+
+    return PDRlegacy::readObstacleFiles
+    (
+        obsFileDir, obsFileNames, groups,
+        meshBb,
+        blocks,
+        cylinders
+    );
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleTypes.C

@@ -0,0 +1,512 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRobstacleTypes.H"
+#include "PDRobstacleTypes.H"
+#include "Enum.H"
+#include "unitConversion.H"
+#include "addToMemberFunctionSelectionTable.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define addObstacleReader(obsType, obsName)                                   \
+    namespace Foam                                                            \
+    {                                                                         \
+    namespace PDRobstacles                                                    \
+    {                                                                         \
+        addNamedToMemberFunctionSelectionTable                                \
+        (                                                                     \
+            PDRobstacle,                                                      \
+            obsType,                                                          \
+            read,                                                             \
+            dictRead,                                                         \
+            obsName                                                           \
+        );                                                                    \
+    }                                                                         \
+    }
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Read porosity, change to blockage. Clamp values [0-1] silently
+static const scalarMinMax limits01(scalarMinMax::zero_one());
+
+// Volume porosity -> blockage
+inline scalar getPorosity(const dictionary& dict)
+{
+    return 1 - limits01.clip(dict.getOrDefault<scalar>("porosity", 0));
+}
+
+// Direction porosities -> blockage
+inline vector getPorosities(const dictionary& dict)
+{
+    vector blockage(vector::one);
+
+    if (dict.readIfPresent("porosities", blockage))
+    {
+        for (scalar& val : blockage)
+        {
+            val = 1 - limits01.clip(val);
+        }
+    }
+
+    return blockage;
+}
+
+
+// Check for "porosity", or "porosities"
+// inline static bool hasPorosity(const dictionary& dict)
+// {
+//     return dict.found("porosity") || dict.found("porosities");
+// }
+
+
+static const Foam::Enum<Foam::vector::components>
+vectorComponentsNames
+({
+    { vector::components::X, "x" },
+    { vector::components::Y, "y" },
+    { vector::components::Z, "z" },
+});
+
+
+enum inletDirnType
+{
+    _X = -1,    // -ve x
+    _Y = -2,    // -ve y
+    _Z = -3,    // -ve z
+    X = 1,      // +ve x
+    Y = 2,      // +ve y
+    Z = 3,      // +ve z
+};
+
+static const Foam::Enum<inletDirnType>
+inletDirnNames
+({
+    { inletDirnType::_X, "-x" },
+    { inletDirnType::_Y, "-y" },
+    { inletDirnType::_Z, "-z" },
+    { inletDirnType::_X, "_x" },
+    { inletDirnType::_Y, "_y" },
+    { inletDirnType::_Z, "_z" },
+    { inletDirnType::X, "+x" },
+    { inletDirnType::Y, "+y" },
+    { inletDirnType::Z, "+z" },
+    { inletDirnType::X, "x" },
+    { inletDirnType::Y, "y" },
+    { inletDirnType::Z, "z" },
+});
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(cylinder, cyl);
+addObstacleReader(cylinder, cylinder);
+
+void Foam::PDRobstacles::cylinder::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    obs.PDRobstacle::readProperties(dict);
+    obs.typeId = enumTypeId;
+
+    // Enforce complete blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+    // if (hasPorosity(dict)) ... warn?
+
+
+    dict.readEntry("point", obs.pt);
+    dict.readEntry("length", obs.len());
+    dict.readEntry("diameter", obs.dia());
+
+    obs.orient = vectorComponentsNames.get("direction", dict);
+
+    // The sortBias for later position sorting
+    switch (obs.orient)
+    {
+        case vector::X:
+            obs.sortBias = obs.len();
+            break;
+
+        default:
+            obs.sortBias = 0.5*obs.dia();
+            break;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(diagbeam, diag);
+addObstacleReader(diagbeam, diagbeam);
+
+void Foam::PDRobstacles::diagbeam::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    obs.PDRobstacle::readProperties(dict);
+    obs.typeId = enumTypeId;
+
+    // Enforce complete blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+    // if (hasPorosity(dict)) ... warn?
+
+
+    dict.readEntry("point", obs.pt);
+    dict.readEntry("length", obs.len());
+    obs.dia() = Zero;
+    obs.theta() = Zero;   // Fix later on
+
+    obs.orient = vectorComponentsNames.get("direction", dict);
+
+    // Angle (degrees) on input, limit to range [0, PI]
+    scalar angle = dict.get<scalar>("angle");
+
+    while (angle > 180) angle -= 180;
+    while (angle < 0) angle += 180;
+
+    labelPair dims;
+    dict.readEntry("width", dims);
+
+    // Swap axes when theta > PI/2
+    // For 89-90 degrees it becomes -ve, which is picked up in following section
+    if (angle > 89)
+    {
+        angle -= 90;
+        dims.flip();  // Swap dimensions
+    }
+
+    obs.theta() = degToRad(angle);
+
+    obs.wa = dims.first();
+    obs.wb = dims.second();
+
+    const scalar ctheta = cos(obs.theta());
+    const scalar stheta = sin(obs.theta());
+
+    // The sortBias for later position sorting
+    switch (obs.orient)
+    {
+        case vector::X:
+            obs.sortBias = obs.len();
+            break;
+
+        case vector::Y:
+            obs.sortBias = 0.5*(obs.wa * stheta + obs.wb * ctheta);
+            break;
+
+        case vector::Z:
+            obs.sortBias = 0.5*(obs.wa * ctheta + obs.wb * stheta);
+            break;
+    }
+
+    // If very nearly aligned with axis, turn it into normal block,
+    // to avoid 1/tan(theta) blowing up
+    if (angle < 1)
+    {
+        Info<< "... changed diag-beam to box" << nl;
+
+        switch (obs.orient)
+        {
+            case vector::X:
+                obs.span = vector(obs.len(), obs.wa, obs.wb);
+                break;
+
+            case vector::Y:
+                obs.span = vector(obs.wb, obs.len(), obs.wa);
+                break;
+
+            case vector::Z:
+                obs.span = vector(obs.wa, obs.wb, obs.len());
+                break;
+        }
+
+        // The pt was end centre (cylinder), now lower corner
+        vector adjustPt = -0.5*obs.span;
+        adjustPt[obs.orient] = 0;
+
+        obs.pt -= adjustPt;
+
+        obs.typeId = PDRobstacles::cuboid::enumTypeId;
+        obs.sortBias = 0;
+        obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1.0;
+        obs.blowoff_type = 0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(cuboid, box);
+
+void Foam::PDRobstacles::cuboid::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    obs.PDRobstacle::readProperties(dict);
+    obs.typeId = enumTypeId;
+
+    // Default - full blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+
+
+    dict.readEntry("point", obs.pt);
+    dict.readEntry("size", obs.span);
+
+    // Optional
+    obs.vbkge = getPorosity(dict);
+
+    // Optional
+    const vector blockages = getPorosities(dict);
+    obs.xbkge = blockages.x();
+    obs.ybkge = blockages.y();
+    obs.zbkge = blockages.z();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(wallbeam, wallbeam);
+
+void Foam::PDRobstacles::wallbeam::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    PDRobstacles::cuboid::read(obs, dict);
+    obs.typeId = enumTypeId;
+
+    // Enforce complete blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+
+    // if (hasPorosity(dict)) ... warn?
+
+    // Additional adjustment for drag etc.
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(grating, grating);
+addObstacleReader(grating, grate);
+
+void Foam::PDRobstacles::grating::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    obs.PDRobstacle::readProperties(dict);
+    obs.typeId = enumTypeId;
+
+    // Initialize to full blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+
+    dict.readEntry("point", obs.pt);
+    dict.readEntry("size", obs.span);
+
+    // TODO: better error handling
+    // if (!equal(cmptProduct(obs.span), 0))
+    // {
+    //     Info<< "Type " << typeId << " has non-zero thickness.";
+    //     ReportLineInfo(lineNo, inputFile);
+    // }
+
+    obs.vbkge = getPorosity(dict);
+
+    const vector blockages = getPorosities(dict);
+    obs.xbkge = blockages.x();
+    obs.ybkge = blockages.y();
+    obs.zbkge = blockages.z();
+
+    // TODO: Warning if porosity was not specified?
+
+
+    // TBD: Default slat width from PDR.params
+    obs.slat_width = dict.getOrDefault<scalar>("slats", Zero);
+
+    // Determine the orientation
+    if (equal(obs.span.x(), 0))
+    {
+        obs.orient = vector::X;
+    }
+    else if (equal(obs.span.y(), 0))
+    {
+        obs.orient = vector::Y;
+    }
+    else
+    {
+        obs.orient = vector::Z;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(louver, louver);
+addObstacleReader(louver, louvre);
+
+void Foam::PDRobstacles::louver::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    obs.PDRobstacle::readProperties(dict);
+    obs.typeId = enumTypeId;
+
+    // Initialize to full blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+
+    dict.readEntry("point", obs.pt);
+    dict.readEntry("size", obs.span);
+
+    // TODO: better error handling
+    // if (!equal(cmptProduct(obs.span), 0))
+    // {
+    //     Info<< "Type " << typeId << " has non-zero thickness.";
+    //     ReportLineInfo(lineNo, inputFile);
+    // }
+
+    obs.vbkge = getPorosity(dict);
+
+    const vector blockages = getPorosities(dict);
+    obs.xbkge = blockages.x();
+    obs.ybkge = blockages.y();
+    obs.zbkge = blockages.z();
+
+    // TODO: Warning if porosity was not specified?
+
+
+    // Blowoff pressure [bar]
+    const scalar blowoffPress = dict.get<scalar>("pressure");
+
+    obs.blowoff_press = barToPa(blowoffPress);
+    obs.blowoff_time = dict.getOrDefault<scalar>("time", 0);
+    obs.blowoff_type = dict.getOrDefault<label>("type", 2);
+
+    if (obs.blowoff_type == 1)
+    {
+        Info<< "Louver : blowoff-type 1 not yet implemented." << nl;
+        // ReportLineInfo(lineNo, inputFile);
+
+        if (obs.blowoff_time != 0)
+        {
+            Info<< "Louver : has blowoff time set,"
+                << " not set to blow off cell-by-cell" << nl;
+            // ReportLineInfo(lineNo, inputFile);
+        }
+    }
+    else
+    {
+        if
+        (
+            (obs.blowoff_type == 1 || obs.blowoff_type == 2)
+         && (blowoffPress > 0)
+        )
+        {
+            if (blowoffPress > maxBlowoffPressure)
+            {
+                Info<< "Blowoff pressure (" << blowoffPress
+                    << ") too high for blowoff type "
+                    << obs.blowoff_type << nl;
+                // ReportLineInfo(lineNo, inputFile);
+            }
+        }
+        else
+        {
+            Info<< "Problem with blowoff parameters" << nl;
+            // ReportLineInfo(lineNo, inputFile);
+            Info<< "Pressure[bar] " << blowoffPress
+                << " Blowoff type " << obs.blowoff_type
+                << ", blowoff pressure " << obs.blowoff_press << nl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+addObstacleReader(patch, patch);
+
+void Foam::PDRobstacles::patch::read
+(
+    PDRobstacle& obs,
+    const dictionary& dict
+)
+{
+    obs.PDRobstacle::readProperties(dict);
+    obs.typeId = enumTypeId;
+
+    const auto nameLen = obs.identifier.length();
+
+    word patchName = word::validate(obs.identifier);
+
+    if (patchName.empty())
+    {
+        FatalErrorInFunction
+            << "RECT_PATCH without a patch name"
+            << exit(FatalError);
+    }
+    else if (patchName.length() != nameLen)
+    {
+        WarningInFunction
+            << "RECT_PATCH stripped invalid characters from patch name: "
+            << obs.identifier
+            << exit(FatalError);
+
+        obs.identifier = std::move(patchName);
+    }
+
+    // Enforce complete blockage
+    obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
+
+    dict.readEntry("point", obs.pt);
+    dict.readEntry("size", obs.span);
+    obs.inlet_dirn = inletDirnNames.get("direction", dict);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#undef addObstacleReader
+
+// ************************************************************************* //

applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleTypes.H

@@ -0,0 +1,300 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::PDRobstacles
+
+Description
+    Reader classes for concrete PDRsetFields obstacle types.
+
+SourceFiles
+    PDRobstacleTypes.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRobstacleTypes_H
+#define PDRobstacleTypes_H
+
+#include "PDRobstacle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PDRobstacles
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class cylinder Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief A cylinder, selectable as \c cyl or \c cylinder
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    point       | The end centre-point                      | yes |
+    length      | The cylinder length                       | yes |
+    diameter    | The cylinder diameter                     | yes |
+    direction   | The x/y/z direction                       | yes |
+\endtable
+
+Example,
+\verbatim
+cyl
+{
+     point       (0 0 0);
+     length      0.95;
+     diameter    0.05;
+     direction   x;
+}
+\endverbatim
+
+\note Does not have porosity.
+*/
+struct cylinder : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::CYLINDER;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+
+/*---------------------------------------------------------------------------*\
+                          Class diagbeam Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief A diagonal beam, which is cylinder-like,
+ selectable as \c diag or \c diagbeam
+
+ If the beam angle is close to zero, it will be changed to a box
+ (PDRobstacles::cuboid)
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    point       | The end centre-point                      | yes |
+    length      | The beam length                           | yes |
+    width       | The beam cross-dimensions                 | yes |
+    angle       | The beam angle (degrees)                  | yes |
+    direction   | The x/y/z direction                       | yes |
+\endtable
+
+Example,
+\verbatim
+diag
+{
+     point       (0 0 0);
+     length      0.95;
+     width       (0.05 0.01);
+     angle       25;
+     direction   x;
+}
+\endverbatim
+
+\note Does not have porosity.
+*/
+struct diagbeam : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::DIAG_BEAM;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+/*---------------------------------------------------------------------------*\
+                           Class cuboid Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief A cuboid, selectable as \c box
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    point       | The lower left corner                     | yes |
+    size        | The (x y z) dimensions                    | yes |
+    porosity    | The volumetric porosity                   | no  | 0
+    porosities  | The directional porosities                | no  | (0 0 0)
+\endtable
+
+Example,
+\verbatim
+box
+{
+     point       (0 0 0);
+     size        (0.05 0.05 2);
+}
+\endverbatim
+*/
+struct cuboid : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::CUBOID;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+
+/*---------------------------------------------------------------------------*\
+                           Class wallbeam Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief A wallbeam, selectable as \c wallbeam which is currently identical
+ to a box (PDRobstacles::cuboid)
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    point       | The lower left corner                     | yes |
+    size        | The (x y z) dimensions                    | yes |
+\endtable
+
+Example,
+\verbatim
+wallbeam
+{
+     point       (0 0 0);
+     size        (0.05 0.05 2);
+}
+\endverbatim
+
+\note Does not have porosity.
+*/
+struct wallbeam : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::WALL_BEAM;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+/*---------------------------------------------------------------------------*\
+                           Class grating Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief A grating, selectable as \c grate or \c grating
+
+The dimensions must include one component that is zero,
+which establishes the orientation.
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    point       | The lower left corner                     | yes |
+    size        | The (x y z) dimensions                    | yes |
+    slats       | The slat width                            | no  | 0
+\endtable
+
+Example,
+\verbatim
+grating
+{
+     point       (0 0 0);
+     size        (0.1 0.1 0);
+     slats       0.005;
+}
+\endverbatim
+*/
+struct grating : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::GRATING;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+/*---------------------------------------------------------------------------*\
+                           Class louver Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief Louver blowoff, selectable as \c louver or \c louvre
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    point       | The lower left corner                     | yes |
+    size        | The (x y z) dimensions                    | yes |
+    pressure    | The blowoff pressure (bar)                | yes |
+    time        | The blowoff time                          | no  | 0
+    type        | The blowoff type (1, 2)                   | no  | 1
+\endtable
+
+Example,
+\verbatim
+louver
+{
+     point       (0 0 0);
+     size        (0.1 0.1);
+     pressure    3;
+     type        1;
+}
+\endverbatim
+*/
+struct louver : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::LOUVER_BLOWOFF;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+/*---------------------------------------------------------------------------*\
+                           Class patch Declaration
+\*---------------------------------------------------------------------------*/
+
+/*! \brief Rectangular patch, selectable as \c patch
+
+Dictionary controls
+\table
+    Property    | Description                               | Required | Default
+    name        | The patch name corner                     | yes |
+    point       | The lower left corner                     | yes |
+    size        | The (x y z) dimensions                    | yes |
+    direction   | The direction (x/y/z, _x/_y/_z or "-x"/"-y"/"-z" | yes |
+\endtable
+
+Example,
+\verbatim
+patch
+{
+     name        inlet;
+     direction   _x;
+     point       (0 0 0);
+     size        (0 0.05 1.0);
+}
+\endverbatim
+*/
+struct patch : public PDRobstacle
+{
+    static constexpr int enumTypeId = legacyTypes::RECT_PATCH;
+    static void read(PDRobstacle& obs, const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace PDRobstacles
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //