From 9d8af49c1661af00c280a1a9b7f036ccefed1c4d Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 16 Oct 2018 13:50:13 +0200
Subject: [PATCH] ENH: parallel and xml output for surfaces (issue #926)

---
 src/meshTools/Make/files                      |   1 +
 src/meshTools/output/foamVtkSurfaceWriter.C   | 477 ++++++++++++++++++
 src/meshTools/output/foamVtkSurfaceWriter.H   | 239 +++++++++
 .../output/foamVtkSurfaceWriterTemplates.C    | 130 +++++
 4 files changed, 847 insertions(+)
 create mode 100644 src/meshTools/output/foamVtkSurfaceWriter.C
 create mode 100644 src/meshTools/output/foamVtkSurfaceWriter.H
 create mode 100644 src/meshTools/output/foamVtkSurfaceWriterTemplates.C

diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index e89cf987fc..4c1116e0b6 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -298,6 +298,7 @@ meshStructure/topoDistanceData.C
 meshStructure/pointTopoDistanceData.C
 
 output/foamVtkIndPatchWriter.C
+output/foamVtkSurfaceWriter.C
 output/foamVtkWriteTopoSet.C
 output/foamVtkWriteFaceSet.C
 output/foamVtkWritePointSet.C
diff --git a/src/meshTools/output/foamVtkSurfaceWriter.C b/src/meshTools/output/foamVtkSurfaceWriter.C
new file mode 100644
index 0000000000..7bf887d232
--- /dev/null
+++ b/src/meshTools/output/foamVtkSurfaceWriter.C
@@ -0,0 +1,477 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "foamVtkSurfaceWriter.H"
+#include "foamVtkOutput.H"
+#include "globalIndex.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::vtk::surfaceWriter::beginPiece()
+{
+    // Basic sizes
+    nLocalPoints_ = points_.size();
+    nLocalFaces_  = faces_.size();
+    nLocalVerts_  = 0;
+
+    for (const face& f : faces_)
+    {
+        nLocalVerts_ += f.size();
+    }
+
+    numberOfPoints_ = nLocalPoints_;
+    numberOfCells_  = nLocalFaces_;
+
+    if (parallel_)
+    {
+        reduce(numberOfPoints_, sumOp<label>());
+        reduce(numberOfCells_,  sumOp<label>());
+    }
+
+
+    // Nothing else to do for legacy
+    if (legacy()) return;
+
+    if (format_)
+    {
+        format().tag
+        (
+            vtk::fileTag::PIECE,
+            vtk::fileAttr::NUMBER_OF_POINTS, numberOfPoints_,
+            vtk::fileAttr::NUMBER_OF_POLYS,  numberOfCells_
+        );
+    }
+}
+
+
+void Foam::vtk::surfaceWriter::writePoints()
+{
+    if (format_)
+    {
+        if (legacy())
+        {
+            legacy::beginPoints(os_, numberOfPoints_);
+        }
+        else
+        {
+            const uint64_t payLoad = vtk::sizeofData<float, 3>(numberOfPoints_);
+
+            format()
+                .tag(vtk::fileTag::POINTS)
+                .beginDataArray<float,3>(vtk::dataArrayAttr::POINTS);
+
+            format().writeSize(payLoad);
+        }
+    }
+
+
+    if (parallel_ ? Pstream::master() : true)
+    {
+        {
+            vtk::writeList(format(), points_);
+        }
+    }
+
+
+    if (parallel_)
+    {
+        if (Pstream::master())
+        {
+            pointField recv;
+
+            // Receive each point field and write
+            for
+            (
+                int slave=Pstream::firstSlave();
+                slave<=Pstream::lastSlave();
+                ++slave
+            )
+            {
+                IPstream fromSlave(Pstream::commsTypes::blocking, slave);
+
+                {
+                    fromSlave >> recv;
+
+                    vtk::writeList(format(), recv);
+                }
+            }
+        }
+        else
+        {
+            // Send to master
+            OPstream toMaster
+            (
+                Pstream::commsTypes::blocking,
+                Pstream::masterNo()
+            );
+
+            {
+                toMaster << points_;
+            }
+        }
+    }
+
+
+    if (format_)
+    {
+        format().flush();
+        format().endDataArray();
+
+        if (!legacy())
+        {
+            format()
+                .endTag(vtk::fileTag::POINTS);
+        }
+    }
+}
+
+
+void Foam::vtk::surfaceWriter::writePolysLegacy
+(
+    const globalIndex& pointOffsets
+)
+{
+    // Connectivity count without additional storage (done internally)
+
+    label nFaces = nLocalFaces_;
+    label nVerts = nLocalVerts_;
+
+    if (parallel_)
+    {
+        reduce(nFaces, sumOp<label>());
+        reduce(nVerts, sumOp<label>());
+    }
+
+    if (nFaces != numberOfCells_)
+    {
+        FatalErrorInFunction
+            << "Expecting " << numberOfCells_
+            << " faces, but found " << nFaces
+            << exit(FatalError);
+    }
+
+    legacy::beginPolys(os_, nFaces, nVerts);
+
+    labelList vertLabels(nLocalFaces_ + nLocalVerts_);
+
+    {
+        // Legacy: size + connectivity together
+        // [nPts, id1, id2, ..., nPts, id1, id2, ...]
+
+        auto iter = vertLabels.begin();
+
+        label off = pointOffsets.localStart();
+
+        {
+            for (const face& f : faces_)
+            {
+                *iter = f.size();       // The size prefix
+                ++iter;
+
+                for (const label pfi : f)
+                {
+                    *iter = pfi + off;  // Face vertex label
+                    ++iter;
+                }
+            }
+            off += points_.size();
+        }
+    }
+
+
+    if (parallel_)
+    {
+        vtk::writeListParallel(format_.ref(), vertLabels);
+    }
+    else
+    {
+        vtk::writeList(format(), vertLabels);
+    }
+
+    if (format_)
+    {
+        format().flush();
+    }
+}
+
+
+void Foam::vtk::surfaceWriter::writePolys
+(
+    const globalIndex& pointOffsets
+)
+{
+    if (format_)
+    {
+        format().tag(vtk::fileTag::POLYS);
+    }
+
+    //
+    // 'connectivity'
+    //
+    {
+        labelList vertLabels(nLocalVerts_);
+
+        label nVerts = nLocalVerts_;
+
+        if (parallel_)
+        {
+            reduce(nVerts, sumOp<label>());
+        }
+
+        if (format_)
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nVerts);
+
+            format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
+            format().writeSize(payLoad * sizeof(label));
+        }
+
+        {
+            // XML: connectivity only
+            // [id1, id2, ..., id1, id2, ...]
+
+            auto iter = vertLabels.begin();
+
+            label off = pointOffsets.localStart();
+
+            {
+                for (const face& f : faces_)
+                {
+                    for (const label pfi : f)
+                    {
+                        *iter = pfi + off;  // Face vertex label
+                        ++iter;
+                    }
+                }
+                off += points_.size();
+            }
+        }
+
+
+        if (parallel_)
+        {
+            vtk::writeListParallel(format_.ref(), vertLabels);
+        }
+        else
+        {
+            vtk::writeList(format(), vertLabels);
+        }
+
+        if (format_)
+        {
+            format().flush();
+            format().endDataArray();
+        }
+    }
+
+
+    //
+    // 'offsets'  (connectivity offsets)
+    //
+    {
+        labelList vertOffsets(nLocalFaces_);
+        label nOffs = vertOffsets.size();
+
+        // global connectivity offsets
+        const globalIndex procOffset(nLocalVerts_);
+
+        if (parallel_)
+        {
+            reduce(nOffs, sumOp<label>());
+        }
+
+        if (format_)
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nOffs);
+
+            format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
+            format().writeSize(payLoad);
+        }
+
+
+        label off = procOffset.localStart();
+
+        auto iter = vertOffsets.begin();
+
+        {
+            for (const face& f : faces_)
+            {
+                off += f.size();   // End offset
+                *iter = off;
+                ++iter;
+            }
+        }
+
+
+        if (parallel_)
+        {
+            vtk::writeListParallel(format_.ref(), vertOffsets);
+        }
+        else
+        {
+            vtk::writeList(format_.ref(), vertOffsets);
+        }
+
+
+        if (format_)
+        {
+            format().flush();
+            format().endDataArray();
+        }
+    }
+
+    if (format_)
+    {
+        format().endTag(vtk::fileTag::POLYS);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::vtk::surfaceWriter::surfaceWriter
+(
+    const pointField& points,
+    const faceList& faces,
+    const vtk::outputOptions opts
+)
+:
+    vtk::fileWriter(vtk::fileTag::POLY_DATA, opts),
+    points_(points),
+    faces_(faces),
+    numberOfPoints_(0),
+    numberOfCells_(0),
+    nLocalPoints_(0),
+    nLocalFaces_(0),
+    nLocalVerts_(0),
+    instant_()
+{
+    // We do not currently support append mode
+    opts_.append(false);
+}
+
+
+Foam::vtk::surfaceWriter::surfaceWriter
+(
+    const pointField& points,
+    const faceList& faces,
+    const fileName& file,
+    bool parallel
+)
+:
+    surfaceWriter(points, faces)
+{
+    open(file, parallel);
+}
+
+
+Foam::vtk::surfaceWriter::surfaceWriter
+(
+    const pointField& points,
+    const faceList& faces,
+    const vtk::outputOptions opts,
+    const fileName& file,
+    bool parallel
+)
+:
+    surfaceWriter(points, faces, opts)
+{
+    open(file, parallel);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::vtk::surfaceWriter::setTime(const instant& inst)
+{
+    instant_ = inst;
+}
+
+
+bool Foam::vtk::surfaceWriter::beginFile(std::string title)
+{
+    if (title.size())
+    {
+        return vtk::fileWriter::beginFile(title);
+    }
+
+    if (instant_.name().size())
+    {
+        return vtk::fileWriter::beginFile
+        (
+            "time='" + instant_.name() + "'"
+        );
+    }
+
+    // Provide default title
+    return vtk::fileWriter::beginFile("surface");
+}
+
+
+bool Foam::vtk::surfaceWriter::writeGeometry()
+{
+    enter_Piece();
+
+    beginPiece();
+
+    writePoints();
+
+    const globalIndex globalPointOffset(nLocalPoints_);
+
+    if (legacy())
+    {
+        writePolysLegacy(globalPointOffset);
+    }
+    else
+    {
+        writePolys(globalPointOffset);
+    }
+
+    return true;
+}
+
+
+bool Foam::vtk::surfaceWriter::beginCellData(label nFields)
+{
+    return enter_CellData(numberOfCells_, nFields);
+}
+
+
+bool Foam::vtk::surfaceWriter::beginPointData(label nFields)
+{
+    return enter_PointData(numberOfPoints_, nFields);
+}
+
+
+void Foam::vtk::surfaceWriter::writeTimeValue()
+{
+    if (instant_.name().size())
+    {
+        vtk::fileWriter::writeTimeValue(instant_.value());
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/output/foamVtkSurfaceWriter.H b/src/meshTools/output/foamVtkSurfaceWriter.H
new file mode 100644
index 0000000000..f6c373559c
--- /dev/null
+++ b/src/meshTools/output/foamVtkSurfaceWriter.H
@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::vtk::surfaceWriter
+
+Description
+    Write faces/points (optionally with fields)
+    as a vtp file or a legacy vtk file.
+
+    The file output states are managed by the Foam::vtk::fileWriter class.
+    FieldData (eg, TimeValue) must appear before any geometry pieces.
+
+Note
+    Parallel output is combined into a single Piece without point merging,
+    which is similar to using multi-piece data sets, but allows more
+    convenient creation as a streaming process.
+    In the future, the duplicate points at processor connections
+    may be addressed using ghost points.
+
+SourceFiles
+    foamVtkSurfaceWriter.C
+    foamVtkSurfaceWriterTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef foamVtkSurfaceWriter_H
+#define foamVtkSurfaceWriter_H
+
+#include "foamVtkFileWriter.H"
+#include "pointField.H"
+#include "faceList.H"
+#include "instant.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class globalIndex;
+
+namespace vtk
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class vtk::surfaceWriter Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceWriter
+:
+    public vtk::fileWriter
+{
+    // Private Member Data
+
+        //- Reference to the points
+        const pointField& points_;
+
+        //- Reference to the faces
+        const faceList& faces_;
+
+        //- The numer of field points for the current Piece
+        label numberOfPoints_;
+
+        //- The numer of field cells (faces) for the current Piece
+        label numberOfCells_;
+
+        //- Local number of points
+        label nLocalPoints_;
+
+        //- Local number of faces
+        label nLocalFaces_;
+
+        //- Local face vertices (connectivity) count. Sum of face sizes.
+        label nLocalVerts_;
+
+        //- Time name/value
+        instant instant_;
+
+
+    // Private Member Functions
+
+        //- Determing sizes (nLocalPoints_, nLocalFaces_, nLocalVerts_),
+        //- and begin piece
+        void beginPiece();
+
+        //- Write patch points
+        void writePoints();
+
+        //- Write patch faces, legacy format
+        void writePolysLegacy(const globalIndex& pointOffsets);
+
+        //- Write patch faces
+        void writePolys(const globalIndex& pointOffsets);
+
+
+        //- No copy construct
+        surfaceWriter(const surfaceWriter&) = delete;
+
+        //- No copy assignment
+        void operator=(const surfaceWriter&) = delete;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components (default format INLINE_BASE64)
+        surfaceWriter
+        (
+            const pointField& pts,
+            const faceList& faces,
+            const vtk::outputOptions opts = vtk::formatType::INLINE_BASE64
+        );
+
+        //- Construct from components (default format INLINE_BASE64),
+        //- and open the file for writing.
+        //  The file name is with/without an extension.
+        surfaceWriter
+        (
+            const pointField& pts,
+            const faceList& faces,
+            const fileName& file,
+            bool parallel = Pstream::parRun()
+        );
+
+        //- Construct from components and open the file for writing.
+        //  The file name is with/without an extension.
+        surfaceWriter
+        (
+            const pointField& pts,
+            const faceList& faces,
+            const vtk::outputOptions opts,
+            const fileName& file,
+            bool parallel = Pstream::parRun()
+        );
+
+
+    //- Destructor
+    virtual ~surfaceWriter() = default;
+
+
+    // Member Functions
+
+        //- File extension for current format type.
+        using vtk::fileWriter::ext;
+
+        //- File extension for given output type
+        inline static word ext(vtk::outputOptions opts)
+        {
+            return opts.ext(vtk::fileTag::POLY_DATA);
+        }
+
+        //- Define a time name/value for the output
+        virtual void setTime(const instant& inst);
+
+        //- Write file header (non-collective)
+        //  \note Expected calling states: (OPENED).
+        virtual bool beginFile(std::string title = "");
+
+        //- Write patch topology
+        //  Also writes the file header if not previously written.
+        //  \note Must be called prior to writing CellData or PointData
+        virtual bool writeGeometry();
+
+        //- Begin CellData output section for specified number of fields.
+        //  Must be called prior to writing any cell data fields.
+        //  \param nFields is for legacy format only.
+        //      When nFields=0, this a no-op for legacy format.
+        //  \note Expected calling states: (PIECE | POINT_DATA).
+        //
+        //  \return True if the state changed
+        virtual bool beginCellData(label nFields = 0);
+
+        //- Begin PointData for specified number of fields.
+        //  Must be called prior to writing any point data fields.
+        //  \param nFields is for legacy format only.
+        //      When nFields=0, this a no-op for legacy format.
+        //  \note Expected calling states: (PIECE | CELL_DATA).
+        //
+        //  \return True if the state changed
+        virtual bool beginPointData(label nFields = 0);
+
+        //- Write "TimeValue" FieldData (name as per Catalyst output)
+        //  Must be called within the FIELD_DATA state.
+        //  \note As a convenience this can also be called from
+        //      (OPENED | DECLARED) states, in which case it invokes
+        //      beginFieldData(1) internally.
+        using vtk::fileWriter::writeTimeValue;
+
+        //- Write the currently set time as "TimeValue" FieldData
+        void writeTimeValue();
+
+
+    // Write
+
+        //- Write a list of Cell (Face) or Point values
+        template<class Type>
+        void write(const word& fieldName, const UList<Type>& field);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace vtk
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "foamVtkSurfaceWriterTemplates.C"
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/output/foamVtkSurfaceWriterTemplates.C b/src/meshTools/output/foamVtkSurfaceWriterTemplates.C
new file mode 100644
index 0000000000..ace78f733f
--- /dev/null
+++ b/src/meshTools/output/foamVtkSurfaceWriterTemplates.C
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::vtk::surfaceWriter::write
+(
+    const word& fieldName,
+    const UList<Type>& field
+)
+{
+    if (isState(outputState::CELL_DATA))
+    {
+        ++nCellData_;
+    }
+    else if (isState(outputState::POINT_DATA))
+    {
+        ++nPointData_;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Bad writer state (" << stateNames[state_]
+            << ") - should be (" << stateNames[outputState::CELL_DATA]
+            << ") or (" << stateNames[outputState::POINT_DATA]
+            << ") for field " << fieldName << nl << endl
+            << exit(FatalError);
+    }
+
+    static_assert
+    (
+        (
+            std::is_same<typename pTraits<Type>::cmptType,label>::value
+         || std::is_floating_point<typename pTraits<Type>::cmptType>::value
+        ),
+        "Label and Floating-point vector space only"
+    );
+
+    const bool isLabel =
+        std::is_same<typename pTraits<Type>::cmptType(), label>::value;
+
+
+    const direction nCmpt(pTraits<Type>::nComponents);
+
+    label nValues = field.size();
+
+    // Could check sizes:
+    //     nValues == nLocalFaces (CELL_DATA)
+    //     nValues == nLocalPoints (POINT_DATA)
+
+    if (parallel_)
+    {
+        reduce(nValues, sumOp<label>());
+    }
+
+    if (format_)
+    {
+        if (isLabel)
+        {
+            if (legacy())
+            {
+                legacy::intField<nCmpt>(format(), fieldName, nValues);
+            }
+            else
+            {
+                const uint64_t payLoad = vtk::sizeofData<label, nCmpt>(nValues);
+
+                format().beginDataArray<label, nCmpt>(fieldName);
+                format().writeSize(payLoad);
+            }
+        }
+        else
+        {
+            if (legacy())
+            {
+                legacy::floatField<nCmpt>(format(), fieldName, nValues);
+            }
+            else
+            {
+                const uint64_t payLoad = vtk::sizeofData<float, nCmpt>(nValues);
+
+                format().beginDataArray<float, nCmpt>(fieldName);
+                format().writeSize(payLoad);
+            }
+        }
+    }
+
+
+    if (parallel_)
+    {
+        vtk::writeListParallel(format_.ref(), field);
+    }
+    else
+    {
+        vtk::writeList(format(), field);
+    }
+
+
+    if (format_)
+    {
+        format().flush();
+        format().endDataArray();
+    }
+}
+
+
+// ************************************************************************* //