From 6e509c10fca93443fd2be05f4d3076d968521177 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 15 Feb 2022 15:36:50 +0100
Subject: [PATCH] ENH: handle manifold cells in VTK, Ensight output (#2396)

- also disables PointData if manifold cells are detected.
  This is a partial workaround for volPointInterpolation problems
  with handling manifold cells.
---
 .../foamToVTK/convertVolumeFields.H           | 11 ++++++++++
 .../ensight/output/ensightOutput.C            | 13 ++++++++----
 .../ensight/part/cells/ensightCells.C         |  5 +++++
 .../ensight/part/cells/ensightCells.H         |  6 ++++++
 .../ensight/part/cells/ensightCellsAddr.C     | 18 +++++++++++------
 .../ensight/part/cells/ensightCellsI.H        |  8 +++++++-
 src/fileFormats/vtk/part/foamVtuSizing.C      | 20 ++++++++++++++-----
 src/fileFormats/vtk/part/foamVtuSizing.H      |  6 ++++++
 src/fileFormats/vtk/part/foamVtuSizingI.H     |  8 +++++++-
 src/fileFormats/vtk/part/foamVtuSizingImpl.C  | 13 ++++++++----
 .../utilities/ensightWrite/ensightWrite.H     |  1 -
 .../utilities/vtkWrite/vtkWrite.C             | 20 +++++++++++++------
 12 files changed, 101 insertions(+), 28 deletions(-)

diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
index db25c78cdc..d224b5d01c 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
@@ -86,6 +86,17 @@ Description
         {
             // Use the appropriate mesh (baseMesh or subMesh)
             vtuMeshCells.reset(meshProxy.mesh());
+
+            if (doPointValues && vtuMeshCells.manifold())
+            {
+                doPointValues = false;
+                nPointFields = 0;
+                Warning
+                    << nl
+                    << "Manifold cells detected (eg, AMI) - disabling PointData"
+                    << nl
+                    << endl;
+            }
         }
 
         internalWriter = autoPtr<vtk::internalWriter>::New
diff --git a/src/fileFormats/ensight/output/ensightOutput.C b/src/fileFormats/ensight/output/ensightOutput.C
index 27e12cba9f..df07a35655 100644
--- a/src/fileFormats/ensight/output/ensightOutput.C
+++ b/src/fileFormats/ensight/output/ensightOutput.C
@@ -31,6 +31,7 @@ License
 #include "face.H"
 #include "polyMesh.H"
 #include "ListOps.H"
+#include "manifoldCellsMeshObject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,7 +81,8 @@ Foam::labelList Foam::ensightOutput::Detail::getPolysNFaces
     const labelUList& addr
 )
 {
-    const cellList& meshCells = mesh.cells();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
 
     labelList list(addr.size());
 
@@ -103,7 +105,8 @@ Foam::labelList Foam::ensightOutput::Detail::getPolysNPointsPerFace
     const labelUList& addr
 )
 {
-    const cellList& meshCells = mesh.cells();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
     const faceList& meshFaces = mesh.faces();
 
     // Count the number of faces per element
@@ -211,7 +214,8 @@ Foam::ensightOutput::Detail::getPolysFacePoints
     const labelList& pointMap
 )
 {
-    const cellList& meshCells = mesh.cells();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
     const faceList& meshFaces = mesh.faces();
     const labelList& owner = mesh.faceOwner();
 
@@ -286,7 +290,8 @@ void Foam::ensightOutput::writePolysPoints
     const labelList& pointMap
 )
 {
-    const cellList& meshCells = mesh.cells();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
     const faceList& meshFaces = mesh.faces();
     const labelList& owner = mesh.faceOwner();
 
diff --git a/src/fileFormats/ensight/part/cells/ensightCells.C b/src/fileFormats/ensight/part/cells/ensightCells.C
index 92164841d6..edbcde79c2 100644
--- a/src/fileFormats/ensight/part/cells/ensightCells.C
+++ b/src/fileFormats/ensight/part/cells/ensightCells.C
@@ -30,6 +30,7 @@ License
 #include "bitSet.H"
 #include "polyMesh.H"
 #include "cellModel.H"
+#include "manifoldCellsMeshObject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -78,6 +79,7 @@ void Foam::ensightCells::resizeAll()
 Foam::ensightCells::ensightCells()
 :
     ensightPart(),
+    manifold_(false),
     offsets_(Zero),
     sizes_(Zero)
 {}
@@ -123,6 +125,7 @@ void Foam::ensightCells::clear()
 
     ensightPart::clear();
 
+    manifold_ = false;
     sizes_ = Zero;
     offsets_ = Zero;
 }
@@ -166,6 +169,8 @@ void Foam::ensightCells::classifyImpl
     const Addressing& cellIds
 )
 {
+    manifold_ = manifoldCellsMeshObject::New(mesh).manifold();
+
     // References to cell shape models
     const cellModel& tet   = cellModel::ref(cellModel::TET);
     const cellModel& pyr   = cellModel::ref(cellModel::PYR);
diff --git a/src/fileFormats/ensight/part/cells/ensightCells.H b/src/fileFormats/ensight/part/cells/ensightCells.H
index c36a6222e6..f2f7edce17 100644
--- a/src/fileFormats/ensight/part/cells/ensightCells.H
+++ b/src/fileFormats/ensight/part/cells/ensightCells.H
@@ -89,6 +89,9 @@ private:
 
     // Private Data
 
+        //- Manifold cells detected
+        bool manifold_;
+
         //- Begin/end offsets for address of each element type
         FixedList<label, nTypes+1> offsets_;
 
@@ -166,6 +169,9 @@ public:
 
     // Access
 
+        //- Manifold mesh cells detected? Globally consistent quantity.
+        inline bool manifold() const noexcept;
+
         //- Processor-local size of all elements.
         using ensightPart::size;
 
diff --git a/src/fileFormats/ensight/part/cells/ensightCellsAddr.C b/src/fileFormats/ensight/part/cells/ensightCellsAddr.C
index 7633ba2679..ff8a7019f8 100644
--- a/src/fileFormats/ensight/part/cells/ensightCellsAddr.C
+++ b/src/fileFormats/ensight/part/cells/ensightCellsAddr.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,23 +30,26 @@ License
 #include "globalIndex.H"
 #include "globalMeshData.H"
 #include "ListOps.H"
+#include "manifoldCellsMeshObject.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::Map<Foam::label>
 Foam::ensightCells::meshPointMap(const polyMesh& mesh) const
 {
-    const label nEstimate = 8*this->size();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
+    const faceList& meshFaces = mesh.faces();
 
-    Map<label> pointMap(nEstimate);
+    Map<label> pointMap(8*this->size());
 
     // Pass 1: markup used points from cells
 
     for (const label celli : this->cellIds())
     {
-        for (const label facei : mesh.cells()[celli])
+        for (const label facei : meshCells[celli])
         {
-            for (const label pointi : mesh.faces()[facei])
+            for (const label pointi : meshFaces[facei])
             {
                 pointMap.insert(pointi, 0);
             }
@@ -72,6 +75,9 @@ Foam::label Foam::ensightCells::meshPointMapppings
     bool parallel
 ) const
 {
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
+
     labelList pointToGlobal;
 
     const bool rewritePointMap = notNull(pointToGlobalRequest);
@@ -177,7 +183,7 @@ Foam::label Foam::ensightCells::meshPointMapppings
 
             for (const label celli : this->cellIds())
             {
-                for (const label facei : mesh.cells()[celli])
+                for (const label facei : meshCells[celli])
                 {
                     for (const label pointi : mesh.faces()[facei])
                     {
diff --git a/src/fileFormats/ensight/part/cells/ensightCellsI.H b/src/fileFormats/ensight/part/cells/ensightCellsI.H
index 18ece2b549..52515d0c14 100644
--- a/src/fileFormats/ensight/part/cells/ensightCellsI.H
+++ b/src/fileFormats/ensight/part/cells/ensightCellsI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,12 @@ inline Foam::label Foam::ensightCells::add(const elemType etype, label id)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline bool Foam::ensightCells::manifold() const noexcept
+{
+    return manifold_;
+}
+
+
 inline const char* Foam::ensightCells::key(const elemType etype)
 {
     return elemNames[etype];
diff --git a/src/fileFormats/vtk/part/foamVtuSizing.C b/src/fileFormats/vtk/part/foamVtuSizing.C
index c1a1688e7a..87a99cac4b 100644
--- a/src/fileFormats/vtk/part/foamVtuSizing.C
+++ b/src/fileFormats/vtk/part/foamVtuSizing.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
 #include "foamVtkCore.H"
 #include "polyMesh.H"
 #include "cellShape.H"
+#include "manifoldCellsMeshObject.H"
 
 // Only used in this file
 #include "foamVtuSizingImpl.C"
@@ -222,6 +223,7 @@ Foam::vtk::vtuSizing::vtuSizing
 void Foam::vtk::vtuSizing::clear() noexcept
 {
     decompose_   = false;
+    manifold_    = false;
     selectionMode_ = FULL_MESH;
     nCells_      = 0;
     nPoints_     = 0;
@@ -263,6 +265,9 @@ void Foam::vtk::vtuSizing::reset
     const cellModel& tetWedge = cellModel::ref(cellModel::TETWEDGE);
 
     const cellShapeList& shapes = mesh.cellShapes();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
+    const faceList& meshFaces = mesh.faces();
 
     // Unique vertex labels per polyhedral
     labelHashSet hashUniqId(2*256);
@@ -285,6 +290,9 @@ void Foam::vtk::vtuSizing::reset
         selectionMode_ = selectionModeType::FULL_MESH;
     }
 
+    // Manifold cells detected?
+    manifold_ = manifoldCellsMeshObject::New(mesh).manifold();
+
     const label nInputCells =
     (
         isSubsetMesh
@@ -337,10 +345,10 @@ void Foam::vtk::vtuSizing::reset
             // Count vertices into first decomposed cell
             bool first = true;
 
-            const cell& cFaces = mesh.cells()[celli];
+            const cell& cFaces = meshCells[celli];
             for (const label facei : cFaces)
             {
-                const face& f = mesh.faces()[facei];
+                const face& f = meshFaces[facei];
 
                 // Face decomposed into triangles and quads
                 // Tri -> Tet, Quad -> Pyr
@@ -365,7 +373,7 @@ void Foam::vtk::vtuSizing::reset
         {
             // Polyhedral: Not decomposed
 
-            const labelList& cFaces = mesh.cells()[celli];
+            const labelList& cFaces = meshCells[celli];
 
             // Unique node ids used (XML/INTERNAL, not needed for LEGACY)
             hashUniqId.clear();
@@ -376,7 +384,7 @@ void Foam::vtk::vtuSizing::reset
 
             for (const label facei : cFaces)
             {
-                const face& f = mesh.faces()[facei];
+                const face& f = meshFaces[facei];
                 nFaceLabels_ += f.size();
 
                 hashUniqId.insert(f);
@@ -411,6 +419,8 @@ void Foam::vtk::vtuSizing::resetShapes
     const cellModel& hex      = cellModel::ref(cellModel::HEX);
 
     decompose_ = false;  // Disallow decomposition
+    manifold_ = false;   // Assume no manifold cells possible
+
     selectionMode_ = SHAPE_MESH;
 
     const label nInputCells = shapes.size();
diff --git a/src/fileFormats/vtk/part/foamVtuSizing.H b/src/fileFormats/vtk/part/foamVtuSizing.H
index 7b2c5a6081..dbedb1f615 100644
--- a/src/fileFormats/vtk/part/foamVtuSizing.H
+++ b/src/fileFormats/vtk/part/foamVtuSizing.H
@@ -169,6 +169,9 @@ private:
         //- Polyhedral decomposition requested
         bool decompose_;
 
+        //- Manifold cells detected
+        bool manifold_;
+
         //- How the mesh cells have been selected or defined
         selectionModeType selectionMode_;
 
@@ -344,6 +347,9 @@ public:
         //- Query the decompose flag (normally off)
         inline bool decompose() const noexcept;
 
+        //- Manifold mesh cells detected? Globally consistent quantity.
+        inline bool manifold() const noexcept;
+
         //- Query how the mesh cells have been selected or defined
         inline selectionModeType selectionMode() const noexcept;
 
diff --git a/src/fileFormats/vtk/part/foamVtuSizingI.H b/src/fileFormats/vtk/part/foamVtuSizingI.H
index 6c20571487..67550e3de1 100644
--- a/src/fileFormats/vtk/part/foamVtuSizingI.H
+++ b/src/fileFormats/vtk/part/foamVtuSizingI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,12 @@ inline bool Foam::vtk::vtuSizing::decompose() const noexcept
 }
 
 
+inline bool Foam::vtk::vtuSizing::manifold() const noexcept
+{
+    return manifold_;
+}
+
+
 inline Foam::vtk::vtuSizing::selectionModeType
 Foam::vtk::vtuSizing::selectionMode() const noexcept
 {
diff --git a/src/fileFormats/vtk/part/foamVtuSizingImpl.C b/src/fileFormats/vtk/part/foamVtuSizingImpl.C
index 69abf42df5..ea83469251 100644
--- a/src/fileFormats/vtk/part/foamVtuSizingImpl.C
+++ b/src/fileFormats/vtk/part/foamVtuSizingImpl.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
 #include "foamVtkCore.H"
 #include "polyMesh.H"
 #include "cellShape.H"
+#include "manifoldCellsMeshObject.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -229,6 +230,9 @@ void Foam::vtk::vtuSizing::populateArrays
     const cellModel& tetWedge = cellModel::ref(cellModel::TETWEDGE);
 
     const cellShapeList& shapes = mesh.cellShapes();
+    ///const cellList& meshCells = mesh.cells();
+    const cellList& meshCells = manifoldCellsMeshObject::New(mesh).cells();
+    const faceList& meshFaces = mesh.faces();
 
     // The face owner is needed to determine the face orientation
     const labelList& owner = mesh.faceOwner();
@@ -451,11 +455,11 @@ void Foam::vtk::vtuSizing::populateArrays
             // Whether to insert cell in place of original or not.
             bool firstCell = true;
 
-            const labelList& cFaces = mesh.cells()[celli];
+            const labelList& cFaces = meshCells[celli];
 
             for (const label facei : cFaces)
             {
-                const face& f = mesh.faces()[facei];
+                const face& f = meshFaces[facei];
                 const bool isOwner = (owner[facei] == celli);
 
                 // Count triangles/quads in decomposition
@@ -579,7 +583,8 @@ void Foam::vtk::vtuSizing::populateArrays
             // face-stream
             //   [nFaces, nFace0Pts, id1, id2, ..., nFace1Pts, id1, id2, ...]
             cellTypes[cellIndex] = vtk::cellType::VTK_POLYHEDRON;
-            const labelList& cFaces = mesh.cells()[celli];
+
+            const labelList& cFaces = meshCells[celli];
 
             const label startLabel = faceIndexer;
 
diff --git a/src/functionObjects/utilities/ensightWrite/ensightWrite.H b/src/functionObjects/utilities/ensightWrite/ensightWrite.H
index f91aadc040..1c4675d48e 100644
--- a/src/functionObjects/utilities/ensightWrite/ensightWrite.H
+++ b/src/functionObjects/utilities/ensightWrite/ensightWrite.H
@@ -98,7 +98,6 @@ Description
         directory   | The output directory name     | no | postProcessing/NAME
         overwrite   | Remove existing directory             | no  | false
         consecutive | Consecutive output numbering          | no  | false
-        nodeValues  | Write values at nodes                 | no  | false
     \endtable
 
     \heading Output Selection
diff --git a/src/functionObjects/utilities/vtkWrite/vtkWrite.C b/src/functionObjects/utilities/vtkWrite/vtkWrite.C
index 66939d61cc..5f92c2a7cd 100644
--- a/src/functionObjects/utilities/vtkWrite/vtkWrite.C
+++ b/src/functionObjects/utilities/vtkWrite/vtkWrite.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -347,15 +347,18 @@ bool Foam::functionObjects::vtkWrite::write()
 
         if (doInternal_)
         {
-            if (interpolate_)
-            {
-                pInterp.reset(new volPointInterpolation(meshProxy.mesh()));
-            }
-
             if (vtuMeshCells.empty())
             {
                 // Use the appropriate mesh (baseMesh or subMesh)
                 vtuMeshCells.reset(meshProxy.mesh());
+
+                if (interpolate_ && vtuMeshCells.manifold())
+                {
+                    interpolate_ = false;
+                    WarningInFunction
+                        << "Manifold cells detected - disabling PointData"
+                        << endl;
+                }
             }
 
             internalWriter = autoPtr<vtk::internalWriter>::New
@@ -385,6 +388,11 @@ bool Foam::functionObjects::vtkWrite::write()
 
             internalWriter->writeTimeValue(timeValue);
             internalWriter->writeGeometry();
+
+            if (interpolate_)
+            {
+                pInterp.reset(new volPointInterpolation(meshProxy.mesh()));
+            }
         }