ENH: improve time-handling for collated ensight surface writer (#1418)

- now maintains a correct list of geometry instances, without assuming
  contiguous file numbering. However, if the numbering is contiguous,
  the more compact case representation will be used.

src/surfMesh/writers/ensight/ensightSurfaceWriter.C

@@ -62,7 +62,7 @@ void Foam::surfaceWriters::ensightWriter::printTimeset
         << "time set:               " << ts << nl
         << "number of steps:        " << 1 << nl;
 
-    // Assume to be contiguous numbering
+    // Single value - starts at index 0
     os  << "filename start number:  0" << nl
         << "filename increment:     1" << nl
         << "time values:" << nl;
@@ -79,25 +79,91 @@ void Foam::surfaceWriters::ensightWriter::printTimeset
     const UList<scalar>& values
 )
 {
-    label count = values.size();
+    label pos_;
+
     os
         << "time set:               " << ts << nl
-        << "number of steps:        " << count << nl;
+        << "number of steps:        " << values.size() << nl;
 
-    // Assume to be contiguous numbering
+    // Assume contiguous numbering - starts at index 0
     os  << "filename start number:  0" << nl
-        << "filename increment:     1" << nl
-        << "time values:" << nl;
+        << "filename increment:     1" << nl;
+
+
+    os  << "time values:" << nl;
+    pos_ = 0;
+    for (const scalar& val : values)
+    {
+        if (pos_ == 6)
+        {
+            os  << nl;
+            pos_ = 0;
+        }
+        ++pos_;
+
+        os  << ' ' << setf(ios_base::right) << setw(12) << val;
+    }
+    os  << nl << nl;
+}
+
+
+void Foam::surfaceWriters::ensightWriter::printTimeset
+(
+    OSstream& os,
+    const label ts,
+    const UList<scalar>& values,
+    const bitSet& indices
+)
+{
+    label pos_;
 
-    count = 0;
-    for (const scalar& t : values)
+    // Check if continuous numbering can be used
+    if
+    (
+        values.empty()
+     || (indices.size() == values.size() && indices.all())
+    )
+    {
+        // Can simply emit as 0-based with increment
+        printTimeset(os, ts, values);
+        return;
+    }
+
+
+    // Generate time set
+    os
+        << "time set:               " << ts << nl
+        << "number of steps:        " << indices.count() << nl;
+
+
+    os  << "filename numbers:" << nl;
+    pos_ = 0;
+    for (const label& idx : indices)
     {
-        os << ' ' << setw(12) << t;
+        if (pos_ == 6)
+        {
+            os  << nl;
+            pos_ = 0;
+        }
+        ++pos_;
+
+        os  << ' ' << setf(ios_base::right) << setw(8) << idx;
+    }
+    os  << nl;
 
-        if (++count % 6 == 0)
+
+    os  << "time values:" << nl;
+    pos_ = 0;
+    for (const label& idx : indices)
+    {
+        if (pos_ == 6)
         {
-            os << nl;
+            os  << nl;
+            pos_ = 0;
         }
+        ++pos_;
+
+        os  << ' ' << setf(ios_base::right) << setw(12) << values[idx];
     }
     os  << nl << nl;
 }
@@ -121,7 +187,7 @@ Foam::surfaceWriters::ensightWriter::ensightWriter
     surfaceWriter(options),
     writeFormat_
     (
-        IOstreamOption::formatNames.lookupOrDefault
+        IOstreamOption::formatNames.getOrDefault
         (
             "format",
             options,
@@ -129,7 +195,7 @@ Foam::surfaceWriters::ensightWriter::ensightWriter
             true  // Failsafe behaviour
         )
     ),
-    collateTimes_(options.lookupOrDefault("collateTimes", true))
+    collateTimes_(options.getOrDefault("collateTimes", true))
 {}
 
 
@@ -164,6 +230,14 @@ Foam::surfaceWriters::ensightWriter::ensightWriter
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::surfaceWriters::ensightWriter::close()
+{
+    times_.clear();
+    meshes_.clear();
+    surfaceWriter::close();
+}
+
+
 // Note that ensight does supports geometry in a separate file,
 // but setting this true leaves mesh files in the wrong places
 // (when there are fields).

src/surfMesh/writers/ensight/ensightSurfaceWriter.H

@@ -47,6 +47,10 @@ Description
         collateTimes | use common geometry for times        | no  | true
     \endtable
 
+    The collated format maintains an internal list of the known times
+    as well as a file-cached version with the field information.
+    The information is used for restarts.
+
 SourceFiles
     ensightSurfaceWriter.C
 
@@ -56,6 +60,8 @@ SourceFiles
 #define ensightSurfaceWriter_H
 
 #include "surfaceWriter.H"
+#include "bitSet.H"
+#include "DynamicList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,10 +86,35 @@ class ensightWriter
         //- Collate times (default: true)
         bool collateTimes_;
 
+        //- The collated output times
+        DynamicList<scalar> times_;
+
+        //- Indices in times_ when a geometry (mesh) has been written
+        bitSet meshes_;
+
+        // NOTE: Could also maintain fieldsDict here,
+        // or as a HashTable of info.
+
 
     // Private Member Functions
 
-        //- Print time-set for ensight case file
+        //- Read time information from baseDir / dictName.
+        //  Returns timeIndex corresponding to timeValue
+        label readPreviousTimes
+        (
+            const fileName& baseDir,
+            const word& dictName,
+            const scalar& timeValue
+        );
+
+        //- The geometry can be any of the following:
+        //
+        // 0: constant/static
+        // 1: moving, with the same frequency as the data
+        // 2: moving, with different frequency as the data
+        int geometryTimeset() const;
+
+        //- Print time-set for ensight case file with a single time
         static void printTimeset
         (
             OSstream& os,
@@ -91,7 +122,17 @@ class ensightWriter
             const scalar& timeValue
         );
 
-        //- Print time-set for ensight case file
+        //- Print time-set for ensight case file, with N times and 0-based
+        //- file numbering
+        //
+        //  \verbatim
+        //  TIME
+        //  time set: ts
+        //  number of steps: ns
+        //  filename start number: 0
+        //  filename increment: 1
+        //  time values: time_1 time_2 ... time_ns
+        //  \endverbatim
         static void printTimeset
         (
             OSstream& os,
@@ -100,6 +141,25 @@ class ensightWriter
         );
 
 
+        //- Print time-set for ensight case file, with N times, 0-based
+        //- file numbering but perhaps non-contiguous
+        //
+        //  \verbatim
+        //  TIME
+        //  time set: ts
+        //  number of steps: ns
+        //  filename numbers: idx_1 idx_2 ... idx_ns
+        //  time values: time_1 time_2 ... time_ns
+        //  \endverbatim
+        static void printTimeset
+        (
+            OSstream& os,
+            const label ts,
+            const UList<scalar>& times,
+            const bitSet& indices
+        );
+
+
         //- Write geometry
         fileName writeCollated();
 
@@ -170,6 +230,10 @@ public:
 
     // Member Functions
 
+        //- Finish output, clears output times.
+        //  Later reuse will rebuild times from fieldsDict file cache.
+        virtual void close(); // override
+
         //- True if the surface format supports geometry in a separate file.
         //  False if geometry and field must be in a single file
         virtual bool separateGeometry() const;

src/surfMesh/writers/ensight/ensightSurfaceWriterCollated.C

@@ -35,6 +35,7 @@ static const equalOp<scalar> equalTimes(ROOTSMALL);
 
 // Use ListOps findLower (with tolerance), to find the location of the next
 // time-related index.
+// The returned index is always 0 or larger (no negative values).
 static label findTimeIndex(const UList<scalar>& list, const scalar val)
 {
     label idx =
@@ -60,6 +61,96 @@ static label findTimeIndex(const UList<scalar>& list, const scalar val)
 } // End namespace Foam
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::label Foam::surfaceWriters::ensightWriter::readPreviousTimes
+(
+    const fileName& baseDir,
+    const word& dictName,
+    const scalar& timeValue
+)
+{
+    // In 1906 and earlier, the fieldsDict contained "meshes" and "times"
+    // entries, each with their own time values.
+    // This makes it more difficult to define the exact correspondence
+    // between geometry intervals and times.
+    //
+    // We now instead track used geometry intervals as a bitSet.
+
+
+    // Only called from master
+
+    label timeIndex = 0;
+
+    labelList geomIndices;
+    scalarList meshTimes;
+
+    dictionary dict;
+
+    const fileName dictFile(baseDir/dictName);
+
+    if (isFile(dictFile))
+    {
+        IFstream is(dictFile);
+
+        if (is.good() && dict.read(is))
+        {
+            meshes_.clear();
+
+            dict.readIfPresent("times", times_);
+            timeIndex = findTimeIndex(times_, timeValue);
+
+            if (dict.readIfPresent("geometry", geomIndices))
+            {
+                // Convert indices to bitSet entries
+                meshes_.set(geomIndices);
+            }
+            else if (dict.readIfPresent("meshes", meshTimes))
+            {
+                WarningInFunction
+                    << nl
+                    << "Setting geometry timeset information from time values"
+                    << " (fieldsDict from an older OpenFOAM version)." << nl
+                    << "This may not be fully reliable." << nl
+                    << nl;
+
+                for (const scalar& meshTime : meshTimes)
+                {
+                    const label meshIndex = findTimeIndex(times_, meshTime);
+                    meshes_.set(meshIndex);
+                }
+            }
+
+            // Make length consistent with time information.
+            // We read/write the indices instead of simply dumping the bitSet.
+            // This makes the contents more human readable.
+            meshes_.resize(times_.size());
+        }
+    }
+
+    return timeIndex;
+}
+
+
+int Foam::surfaceWriters::ensightWriter::geometryTimeset() const
+{
+    if (meshes_.count() <= 1)
+    {
+        // Static
+        return 0;
+    }
+
+    if (meshes_.size() == times_.size() && meshes_.all())
+    {
+        // Geometry changing is the same as fields changing
+        return 1;
+    }
+
+    // Geometry changing differently from fields
+    return 2;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::fileName Foam::surfaceWriters::ensightWriter::writeCollated()
@@ -121,6 +212,7 @@ Foam::fileName Foam::surfaceWriters::ensightWriter::writeCollated
     // Mesh changed since last output? Do before any merging.
     const bool meshChanged = (!upToDate_);
 
+
     // geometry merge() implicit
     tmp<Field<Type>> tfield = mergeField(localValues);
 
@@ -133,96 +225,95 @@ Foam::fileName Foam::surfaceWriters::ensightWriter::writeCollated
             mkDir(outputFile.path());
         }
 
-        scalar meshValue = timeValue;
-        label meshIndex = 0;
-        label timeIndex = 0;
+        bool stateChanged = meshChanged;
+
+        const label timeIndex =
+        (
+            times_.empty()
+          ? readPreviousTimes(baseDir, "fieldsDict", timeValue)
+          : findTimeIndex(times_, timeValue)
+        );
+
+
+        // Update stored times list and mesh index
+
+        if (timeIndex < meshes_.size()-1)
+        {
+            // Clear old content when shrinking
+            meshes_.unset(timeIndex);
+        }
+
+        // Extend or truncate list
+        meshes_.resize(timeIndex+1);
+        times_.resize(timeIndex+1, VGREAT);
 
-        fileName geometryName;
+        if (meshChanged)
+        {
+            meshes_.set(timeIndex);
+        }
 
-        // Do case file
+        if (!equalTimes(times_[timeIndex], timeValue))
         {
-            dictionary dict;
-            scalarList meshes;
-            scalarList times;
-            bool stateChanged = meshChanged;
+            stateChanged = true;
+            times_[timeIndex] = timeValue;
+        }
 
-            if (isFile(baseDir/"fieldsDict"))
-            {
-                IFstream is(baseDir/"fieldsDict");
-                if (is.good() && dict.read(is))
-                {
-                    dict.readIfPresent("meshes", meshes);
-                    dict.readIfPresent("times", times);
-
-                    timeIndex = findTimeIndex(times, timeValue);
-
-                    if (meshChanged)
-                    {
-                        meshValue = timeValue;
-                        meshIndex = findTimeIndex(meshes, meshValue);
-                    }
-                    else if (meshes.size())
-                    {
-                        meshIndex = meshes.size()-1;
-                        meshValue = meshes.last();
-                    }
-                }
-            }
 
-            // Update stored times list
-            meshes.resize(meshIndex+1, -1);
-            times.resize(timeIndex+1, -1);
+        // The most current geometry index
+        const label geomIndex(max(0, meshes_.find_last()));
 
-            if
-            (
-                !equalTimes(meshes[meshIndex], meshValue)
-             || !equalTimes(times[timeIndex], timeValue)
-            )
-            {
-                stateChanged = true;
-            }
 
-            meshes[meshIndex] = meshValue;
-            times[timeIndex] = timeValue;
+        // This will be used for the name of a static geometry,
+        // or just the masking part for moving geometries.
+        const fileName geometryName
+        (
+            "data"/word::printf(fmt, geomIndex)/ensightCase::geometryName
+        );
 
-            geometryName =
-                "data"/word::printf(fmt, meshIndex)/ensightCase::geometryName;
 
+        // Do case file
+        {
+            dictionary dict;
 
-            // Add my information to dictionary
-            {
-                dict.set("meshes", meshes);
-                dict.set("times", times);
-                if (dict.found("fields"))
-                {
-                    dictionary& fieldsDict = dict.subDict("fields");
-                    if (!fieldsDict.found(fieldName))
-                    {
-                        dictionary fieldDict;
-                        fieldDict.set("type", ensightPTraits<Type>::typeName);
-                        fieldDict.set("name", varName); // ensight variable name
+            // Add time information to dictionary
+            dict.set("geometry", meshes_.sortedToc());
+            dict.set("times", times_);
 
-                        fieldsDict.set(fieldName, fieldDict);
+            // Debugging, or if needed for older versions:
+            //// dict.set
+            //// (
+            ////     "meshes",
+            ////     IndirectList<scalar>(times_, meshes_.sortedToc())
+            //// );
 
-                        stateChanged = true;
-                    }
-                }
-                else
-                {
-                    dictionary fieldDict;
-                    fieldDict.set("type", ensightPTraits<Type>::typeName);
-                    fieldDict.set("name", varName); // ensight variable name
 
-                    dictionary fieldsDict;
-                    fieldsDict.set(fieldName, fieldDict);
+            // Add field information to dictionary?
+            bool hasField = false;
 
-                    dict.set("fields", fieldsDict);
+            dictionary* fieldsDictPtr = dict.findDict("fields");
+            if (fieldsDictPtr)
+            {
+                hasField = fieldsDictPtr->found(fieldName);
+            }
+            else
+            {
+                // No "fields" dictionary - create first
+                fieldsDictPtr = dict.set("fields", dictionary())->dictPtr();
+            }
 
-                    stateChanged = true;
-                }
+            if (!hasField)
+            {
+                dictionary fieldDict;
+                fieldDict.set("type", ensightPTraits<Type>::typeName);
+                fieldDict.set("name", varName); // ensight variable name
+
+                fieldsDictPtr->set(fieldName, fieldDict);
+                stateChanged = true;
             }
 
 
+            const dictionary& fieldsDict = *fieldsDictPtr;
+
             if (stateChanged)
             {
                 if (verbose_)
@@ -252,8 +343,7 @@ Foam::fileName Foam::surfaceWriters::ensightWriter::writeCollated
                 // 1: moving, with the same frequency as the data
                 // 2: moving, with different frequency as the data
 
-                const label tsGeom =
-                    (meshes.size() == 1 ? 0 : meshes == times ? 1 : 2);
+                const label tsGeom = geometryTimeset();
 
                 osCase
                     << "FORMAT" << nl
@@ -281,16 +371,19 @@ Foam::fileName Foam::surfaceWriters::ensightWriter::writeCollated
                     << nl
                     << "VARIABLE" << nl;
 
-                const dictionary& fieldsDict = dict.subDict("fields");
+
                 for (const entry& dEntry : fieldsDict)
                 {
                     const dictionary& subDict = dEntry.dict();
 
                     const word fieldType(subDict.get<word>("type"));
-                    const word varName = subDict.lookupOrDefault
+                    const word varName
                     (
-                        "name",
-                        dEntry.keyword() // fieldName as fallback
+                        subDict.getOrDefault<word>
+                        (
+                            "name",
+                            dEntry.keyword() // fieldName as fallback
+                        )
                     );
 
                     osCase
@@ -309,10 +402,10 @@ Foam::fileName Foam::surfaceWriters::ensightWriter::writeCollated
                     << nl
                     << "TIME" << nl;
 
-                printTimeset(osCase, 1, times);
+                printTimeset(osCase, 1, times_);
                 if (tsGeom == 2)
                 {
-                    printTimeset(osCase, 2, meshes);
+                    printTimeset(osCase, 2, times_, meshes_);
                 }
 
                 osCase << "# end" << nl;