From 2ca47927346de8736a6e8d8eed5fddaa78c5e746 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 21 Aug 2019 12:13:51 +0200
Subject: [PATCH] TUT: add coded example for sampling cellZoneId (#1383)

---
 .../angledDuct/implicit/system/controlDict    |   5 +
 .../angledDuct/implicit/system/sampling       | 123 ++++++++++++++++++
 2 files changed, 128 insertions(+)
 create mode 100644 tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/sampling

diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/controlDict b/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/controlDict
index 247b89daa2..8f5cee03fe 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/controlDict
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/controlDict
@@ -47,5 +47,10 @@ graphFormat     raw;
 
 runTimeModifiable true;
 
+functions
+{
+    #include "sampling"
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/sampling b/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/sampling
new file mode 100644
index 0000000000..6bc32218c5
--- /dev/null
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuct/implicit/system/sampling
@@ -0,0 +1,123 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+cellZoneID
+{
+    type    coded;
+    libs    ("utilityFunctionObjects");
+    name    cellZoneID;
+
+    executeControl  writeTime;
+
+    localCode
+    #{
+        static const word fieldName("cellZoneID");
+
+        // Create and populate a "cellZoneID" volScalarField for cellZones.
+        // Return nullptr if there are not any such zones.
+        //
+        // Could be improved for mesh motion etc
+        // but is fairly low overhead anyhow
+        volScalarField* getZoneField(const fvMesh& mesh)
+        {
+            volScalarField* volZonePtr =
+                mesh.getObjectPtr<volScalarField>(fieldName);
+
+            const cellZoneMesh& zones = mesh.cellZones();
+
+            if (!zones.empty())
+            {
+                if (!volZonePtr)
+                {
+                    volZonePtr = new volScalarField
+                    (
+                        IOobject
+                        (
+                            fieldName,
+                            mesh.time().timeName(),
+                            mesh,
+                            IOobject::NO_READ,
+                            IOobject::NO_WRITE,
+                            true  // Register
+                        ),
+                        mesh,
+                        dimless,
+                        zeroGradientFvPatchField<scalar>::typeName
+                    );
+
+                    volZonePtr->store();
+
+                    Info<< "Creating " << fieldName
+                        << " field for postProcessing" << nl;
+                }
+            }
+            else
+            {
+                Info<< "No cellZones - not creating " << fieldName
+                    << " field for postProcessing" << nl;
+            }
+
+            if (volZonePtr)
+            {
+                // Copy zone id as scalar.
+                // For consistency, do this whenever the volField exists,
+                // even if there are no zones.
+
+                auto& field = *volZonePtr;
+
+                auto& internal = field.primitiveFieldRef();
+                internal = scalar(-1);
+
+                for (const cellZone& zn : zones)
+                {
+                    UIndirectList<scalar>(internal, zn) = scalar(zn.index());
+                }
+
+                field.correctBoundaryConditions();
+            }
+
+            return volZonePtr;
+        }
+    #};
+
+    codeExecute
+    #{
+        // Don't need the return value, just the contents on the registry.
+        (void) getZoneField(mesh());
+    #};
+}
+
+
+plane
+{
+    type    surfaces;
+    libs    ("sampling");
+
+    writeControl    writeTime;
+
+    surfaceFormat   vtk;
+    fields          ( cellZoneID );
+
+    surfaces
+    (
+        zNormal
+        {
+            type            cuttingPlane;
+            planeType       pointAndNormal;
+            pointAndNormalDict
+            {
+                point   (0 0 0);
+                normal  (0 0 1);
+            }
+            interpolate     false;
+        }
+    );
+}
+
+
+// ************************************************************************* //