ENH: improvements for makeFaMesh, checkFaMesh

- added -dry-run, -write-vtk options.
  Additional mesh information after creation.

- add parallel reductions and more information for checkFaMesh

ENH: minor cleanup of faPatch internals

- align pointLabels and pointEdges creation more closely with coding
  patterns used in PrimitivePatch

- use fileHandler when loading "S0" field.

applications/utilities/finiteArea/checkFaMesh/Make/options

@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
-    -lfiniteVolume \
+    -lfileFormats \
+    -lsurfMesh \
     -lmeshTools

applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,14 +30,22 @@ Application
 Description
     Check a finiteArea mesh
 
-Author
+Original Authors
     Zeljko Tukovic, FAMENA
     Hrvoje Jasak, Wikki Ltd.
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
-#include "faCFD.H"
+#include "Time.H"
+#include "argList.H"
+#include "faMesh.H"
+#include "polyMesh.H"
+#include "areaFaMesh.H"
+#include "edgeFaMesh.H"
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "processorFaPatch.H"
+#include "foamVtkUIndPatchWriter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,37 +60,32 @@ int main(int argc, char *argv[])
         "Check a finiteArea mesh"
     );
 
-    #include "addRegionOption.H"
+    argList::addBoolOption
+    (
+        "write-vtk",
+        "Write mesh as a vtp (vtk) file for display or debugging"
+    );
 
+    #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createNamedMesh.H"
-    #include "createFaMesh.H"
+    #include "createNamedPolyMesh.H"
+
+
+    // Create
+    faMesh aMesh(mesh);
 
     Info<< "Time = " << runTime.timeName() << nl << endl;
 
-    // General mesh statistics
-    Info<< "Number of points: " << aMesh.nPoints() << nl
-        << "Number of internal edges: " << aMesh.nInternalEdges() << nl
-        << "Number of edges: " << aMesh.nEdges() << nl
-        << "Number of faces: " << aMesh.nFaces() << nl
-        << endl;
+    #include "printMeshSummary.H"
 
-    // Check geometry
-    Info<< "Face area:  min = " << min(aMesh.S().field())
-        << " max = "  << max(aMesh.S().field()) << nl
-        << "Internal edge length: min = "
-        << min(aMesh.magLe().internalField()) << nl
-        << " max = "  << max(aMesh.magLe().internalField()) << nl
-        << "Edge length: min = "
-        << min(aMesh.magLe()).value() << nl
-        << " max = "  << max(aMesh.magLe()).value() << nl
-        << "Face area normals:  min = " << min(aMesh.faceAreaNormals().field())
-        << " max = "  << max(aMesh.faceAreaNormals().field()) << nl
-        << endl;
+    if (args.found("write-vtk"))
+    {
+        #include "faMeshWriteVTK.H"
+    }
 
+    Info<< "\nEnd\n" << endl;
 
-    Info << "\nEnd" << endl;
     return 0;
 }
 

applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    VTK output of faMesh with some geometric or debug fields
+
+\*---------------------------------------------------------------------------*/
+
+{
+    // finiteArea
+    vtk::uindirectPatchWriter writer
+    (
+        aMesh.patch(),
+        // vtk::formatType::INLINE_ASCII,
+        fileName
+        (
+            aMesh.mesh().time().globalPath() / "finiteArea"
+        )
+    );
+
+    writer.writeGeometry();
+
+    // CellData
+    writer.beginCellData(3);
+    writer.writeProcIDs();
+    writer.write("area", aMesh.S().field());
+    writer.write("normal", aMesh.faceAreaNormals());
+
+    // PointData
+    writer.beginPointData(1);
+    writer.write("normal", aMesh.pointAreaNormals());
+
+    Info<< nl
+        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+}
+
+
+// ************************************************************************* //

applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Summary of faMesh information
+
+\*---------------------------------------------------------------------------*/
+
+{
+    const faBoundaryMesh& patches = aMesh.boundary();
+    const label nNonProcessor = patches.nNonProcessor();
+    const label nPatches = patches.size();
+
+    Info<< "----------------" << nl
+        << "Mesh Information" << nl
+        << "----------------" << nl
+        << "  " << "boundingBox: " << boundBox(aMesh.points()) << nl;
+
+    Info<< "  Number of points: "
+        << returnReduce(aMesh.nPoints(), sumOp<label>()) << nl
+        << "  Number of faces: "
+        << returnReduce(aMesh.nFaces(), sumOp<label>()) << nl;
+
+    Info<< "  Number of edges: "
+        << returnReduce(aMesh.nEdges(), sumOp<label>()) << nl
+        << "  Number of internal edges: "
+        << returnReduce(aMesh.nInternalEdges(), sumOp<label>()) << nl;
+
+
+    label nProcEdges = 0;
+
+    if (Pstream::parRun())
+    {
+        for (const faPatch& fap : patches)
+        {
+            const auto* cpp = isA<processorFaPatch>(fap);
+
+            if (cpp)
+            {
+                nProcEdges += fap.nEdges();
+            }
+        }
+    }
+
+    const label nBoundEdges = aMesh.nBoundaryEdges() - nProcEdges;
+
+    Info<< "  Number of boundary edges: "
+        << returnReduce(nBoundEdges - nProcEdges, sumOp<label>()) << nl;
+
+    if (Pstream::parRun())
+    {
+        Info<< "  Number of processor edges: "
+            << returnReduce(nProcEdges, sumOp<label>()) << nl;
+    }
+
+
+    Info<< "----------------" << nl
+        << "Patches" << nl
+        << "----------------" << nl;
+
+    for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+    {
+        const faPatch& p = patches[patchi];
+
+         Info<< "  " << "patch " << p.index()
+             << " (size: " << returnReduce(p.size(), sumOp<label>())
+             << ") name: " << p.name()
+             << nl;
+    }
+
+
+    // Geometry information
+    Info<< nl;
+    {
+        scalarMinMax limit(gMinMax(aMesh.S().field()));
+        Info<< "Face area:" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+    }
+
+    {
+        scalarMinMax limit(minMax(aMesh.magLe().primitiveField()));
+
+        // Include processor boundaries into 'internal' edges
+        if (Pstream::parRun())
+        {
+            for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
+            {
+                limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+            }
+
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length (internal):" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+
+
+        // Include (non-processor) boundaries
+        for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+        {
+            limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length:" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+    }
+
+    // Not particularly meaningful
+    #if 0
+    {
+        MinMax<vector> limit(gMinMax(aMesh.faceAreaNormals().field()));
+
+        Info<< "Face area normals:" << nl
+            << "    min = " << limit.min() << " max = " << limit.max() << nl;
+    }
+    #endif
+}
+
+
+// ************************************************************************* //

applications/utilities/finiteArea/makeFaMesh/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/faDecompose/lnInclude \
     -I$(LIB_SRC)/parallel/reconstruct/faReconstruct/lnInclude
@@ -8,6 +9,7 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteArea \
     -lfileFormats \
+    -lsurfMesh \
     -lmeshTools \
     -lfaDecompose \
     -lfaReconstruct

applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H

@@ -17,9 +17,11 @@ Description
 
 if (Pstream::parRun())
 {
-    faMeshReconstructor reconstructor(areaMesh);
+    faMeshReconstructor reconstructor(aMesh);
     reconstructor.writeAddressing();
 
+    Info<< "Wrote proc-addressing" << nl << endl;
+
     // Handle area fields
     // ------------------
 
@@ -72,22 +74,53 @@ if (Pstream::parRun())
 
     if (nAreaFields)
     {
-        Info<< "Decomposing " << nAreaFields << " area fields" << nl << endl;
+        Info<< "Decomposing " << nAreaFields << " area fields" << nl;
 
         faFieldDecomposer fieldDecomposer
         (
             fullMesh,
-            areaMesh,
+            aMesh,
             reconstructor.edgeProcAddressing(),
             reconstructor.faceProcAddressing(),
             reconstructor.boundaryProcAddressing()
         );
 
-        fieldDecomposer.decomposeFields(areaScalarFields);
-        fieldDecomposer.decomposeFields(areaVectorFields);
-        fieldDecomposer.decomposeFields(areaSphTensorFields);
-        fieldDecomposer.decomposeFields(areaSymmTensorFields);
-        fieldDecomposer.decomposeFields(areaTensorFields);
+        if (areaScalarFields.size())
+        {
+            Info<< "  scalars: "
+                << flatOutput(PtrListOps::names(areaScalarFields)) << nl;
+            fieldDecomposer.decomposeFields(areaScalarFields);
+        }
+
+        if (areaVectorFields.size())
+        {
+            Info<< "  vectors: "
+                << flatOutput(PtrListOps::names(areaVectorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaVectorFields);
+        }
+
+        if (areaSphTensorFields.size())
+        {
+            Info<< "  sphTensors: "
+                << flatOutput(PtrListOps::names(areaSphTensorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaSphTensorFields);
+        }
+
+        if (areaSymmTensorFields.size())
+        {
+            Info<< "  symmTensors: "
+                << flatOutput(PtrListOps::names(areaSymmTensorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaSymmTensorFields);
+        }
+
+        if (areaTensorFields.size())
+        {
+            Info<< "  tensors: "
+                << flatOutput(PtrListOps::names(areaTensorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaTensorFields);
+        }
+
+        Info<< endl;
     }
 }
 

applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H

@@ -16,15 +16,18 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    Info<< "Writing edges in obj format" << endl;
+    Info<< nl
+        << "Writing edges in obj format" << endl;
 
-    word outputName("faMesh-edges.obj");
+    word outputName("finiteArea-edges.obj");
 
     if (Pstream::parRun())
     {
         outputName = word
         (
-            "faMesh-edges-" + Foam::name(Pstream::myProcNo()) + ".obj"
+            "finiteArea-edges-proc"
+          + Foam::name(Pstream::myProcNo())
+          + ".obj"
         );
     }
 
@@ -36,7 +39,7 @@ Description
         false
     );
 
-    os.write(areaMesh.patch().edges(), areaMesh.patch().localPoints());
+    os.write(aMesh.patch().edges(), aMesh.patch().localPoints());
 }
 
 

applications/utilities/finiteArea/makeFaMesh/faMeshWriteVTK.H

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    VTK output of faMesh with some geometric or debug fields
+
+\*---------------------------------------------------------------------------*/
+
+{
+    // finiteArea
+    vtk::uindirectPatchWriter writer
+    (
+        aMesh.patch(),
+        // vtk::formatType::INLINE_ASCII,
+        fileName
+        (
+            aMesh.mesh().time().globalPath() / "finiteArea"
+        )
+    );
+
+    writer.writeGeometry();
+
+    // CellData
+    writer.beginCellData(3);
+    writer.writeProcIDs();
+    writer.write("area", aMesh.S().field());
+    writer.write("normal", aMesh.faceAreaNormals());
+
+    // PointData
+    writer.beginPointData(1);
+    writer.write("normal", aMesh.pointAreaNormals());
+
+    Info<< nl
+        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+}
+
+
+// ************************************************************************* //

applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C

@@ -32,7 +32,7 @@ Description
     When called in parallel, it will also try to act like decomposePar,
     create procAddressing and decompose serial finite-area fields.
 
-Author
+Original Authors
     Zeljko Tukovic, FAMENA
     Hrvoje Jasak, Wikki Ltd.
 
@@ -48,6 +48,8 @@ Author
 #include "areaFields.H"
 #include "faFieldDecomposer.H"
 #include "faMeshReconstructor.H"
+#include "PtrListOps.H"
+#include "foamVtkUIndPatchWriter.H"
 #include "OBJstream.H"
 
 using namespace Foam;
@@ -69,11 +71,21 @@ int main(int argc, char *argv[])
     );
     argList::addOption("dict", "file", "Alternative faMeshDefinition");
 
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Create but do not write"
+    );
+    argList::addBoolOption
+    (
+        "write-vtk",
+        "Write mesh as a vtp (vtk) file for display or debugging"
+    );
     argList::addBoolOption
     (
         "write-edges-obj",
-        "Write mesh edges as obj files and exit",
-        false  // could make an advanced option
+        "Write mesh edges as obj files (one per processor)",
+        true  // advanced option (debugging only)
     );
 
     #include "addRegionOption.H"
@@ -81,6 +93,8 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createNamedPolyMesh.H"
 
+    const bool dryrun = args.found("dry-run");
+
     // Reading faMeshDefinition dictionary
     #include "findMeshDefinitionDict.H"
 
@@ -91,33 +105,40 @@ int main(int argc, char *argv[])
         meshDefDict.add("emptyPatch", patchName, true);
     }
 
-    // Create
-    faMesh areaMesh(mesh, meshDefDict);
 
-    bool quickExit = false;
+    // Create
+    faMesh aMesh(mesh, meshDefDict);
+
+    // Mesh information
+    #include "printMeshSummary.H"
 
     if (args.found("write-edges-obj"))
     {
-        quickExit = true;
         #include "faMeshWriteEdgesOBJ.H"
     }
 
-    if (quickExit)
+    if (args.found("write-vtk"))
     {
-        Info<< "\nEnd\n" << endl;
-        return 0;
+        #include "faMeshWriteVTK.H"
     }
 
-    // Set the precision of the points data to 10
-    IOstream::defaultPrecision(10);
+    if (dryrun)
+    {
+        Info<< "\ndry-run: not writing mesh or decomposing fields\n" << nl;
+    }
+    else
+    {
+        // Set the precision of the points data to 10
+        IOstream::defaultPrecision(10);
 
-    Info<< nl << "Write finite area mesh." << nl;
-    areaMesh.write();
-    Info<< endl;
+        Info<< nl << "Write finite area mesh." << nl;
+        aMesh.write();
 
-    #include "decomposeFaFields.H"
+        Info<< endl;
+        #include "decomposeFaFields.H"
+    }
 
-    Info << "\nEnd\n" << endl;
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }

applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Summary of faMesh information
+
+\*---------------------------------------------------------------------------*/
+
+{
+    const faBoundaryMesh& patches = aMesh.boundary();
+    const label nNonProcessor = patches.nNonProcessor();
+    const label nPatches = patches.size();
+
+    Info<< "----------------" << nl
+        << "Mesh Information" << nl
+        << "----------------" << nl
+        << "  " << "boundingBox: " << boundBox(aMesh.points()) << nl
+        << "  " << "nFaces: " << returnReduce(aMesh.nFaces(), sumOp<label>())
+        << nl;
+
+
+    Info<< "----------------" << nl
+        << "Patches" << nl
+        << "----------------" << nl;
+
+    for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+    {
+        const faPatch& p = patches[patchi];
+
+         Info<< "  " << "patch " << p.index()
+             << " (size: " << returnReduce(p.size(), sumOp<label>())
+             << ") name: " << p.name()
+             << nl;
+    }
+
+    // Geometry information
+    Info<< nl;
+    {
+        scalarMinMax limit(gMinMax(aMesh.S().field()));
+        Info<< "Face area:" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+    }
+
+    {
+        scalarMinMax limit(minMax(aMesh.magLe().primitiveField()));
+
+        // Include processor boundaries into 'internal' edges
+        if (Pstream::parRun())
+        {
+            for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
+            {
+                limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+            }
+
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length (internal):" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+
+
+        // Include (non-processor) boundaries
+        for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+        {
+            limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length:" << nl
+            << "    min = " << limit.min()
+            << " max = "  << limit.max() << nl;
+    }
+
+    // Not particularly meaningful
+    #if 0
+    {
+        MinMax<vector> limit(gMinMax(aMesh.faceAreaNormals().field()));
+
+        Info<< "Face area normals:" << nl
+            << "    min = " << limit.min() << " max = " << limit.max() << nl;
+    }
+    #endif
+}
+
+
+// ************************************************************************* //