From 3944c7a628e3e55668d44f7071ea9f60bc26b3f1 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Fri, 24 Nov 2023 16:01:56 +0000
Subject: [PATCH 1/3] ENH: coordSetWriter - brought across field scaling and
 offsets from surface writers

Example:

    formatOptions
    {
        <writer>
        {
            // Apply offsets to field values
            fieldLevel
            {
                T   273.15; // Convert from K to C by subtracting 273.15
            }

            // Note: scale applied after application of field level
            fieldScale
            {
                p   0.001;  // Convert pressure from Pa to kPa by scaling by 0.001
            }
        }
    }
---
 .../coordSet/writers/common/coordSetWriter.C  |  28 ++++-
 .../coordSet/writers/common/coordSetWriter.H  |  23 ++++
 .../writers/common/coordSetWriterTemplates.C  | 102 ++++++++++++++++++
 3 files changed, 152 insertions(+), 1 deletion(-)

diff --git a/src/meshTools/coordSet/writers/common/coordSetWriter.C b/src/meshTools/coordSet/writers/common/coordSetWriter.C
index d3e7c5cae4..9ea4cd5775 100644
--- a/src/meshTools/coordSet/writers/common/coordSetWriter.C
+++ b/src/meshTools/coordSet/writers/common/coordSetWriter.C
@@ -117,7 +117,12 @@ Foam::coordSetWriter::coordSetWriter()
     verbose_(false),
     nFields_(0),
     currTime_(),
-    outputPath_()
+    outputPath_(),
+    geometryScale_(1),
+    geometryCentre_(Zero),
+    geometryTransform_(),
+    fieldLevel_(),
+    fieldScale_()
 {}
 
 
@@ -126,6 +131,27 @@ Foam::coordSetWriter::coordSetWriter(const dictionary& options)
     coordSetWriter()
 {
     options.readIfPresent("verbose", verbose_);
+
+    geometryScale_ = 1;
+    geometryCentre_ = Zero;
+    geometryTransform_.clear();
+
+    options.readIfPresent("scale", geometryScale_);
+
+    // Optional cartesian coordinate system transform
+    const auto* dictptr = options.findDict("transform", keyType::LITERAL);
+
+    if (dictptr)
+    {
+        dictptr->readIfPresent("rotationCentre", geometryCentre_);
+
+        // 'origin' is optional within sub-dictionary
+        geometryTransform_ =
+            coordSystem::cartesian(*dictptr, IOobjectOption::LAZY_READ);
+    }
+
+    fieldLevel_ = options.subOrEmptyDict("fieldLevel");
+    fieldScale_ = options.subOrEmptyDict("fieldScale");
 }
 
 
diff --git a/src/meshTools/coordSet/writers/common/coordSetWriter.H b/src/meshTools/coordSet/writers/common/coordSetWriter.H
index f98cf887b1..8069fadb7f 100644
--- a/src/meshTools/coordSet/writers/common/coordSetWriter.H
+++ b/src/meshTools/coordSet/writers/common/coordSetWriter.H
@@ -66,6 +66,7 @@ SourceFiles
 #include "UPtrList.H"
 #include "instant.H"
 #include "InfoProxy.H"
+#include "cartesianCS.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -125,6 +126,21 @@ protected:
         //- The full output directory and file (coords) name
         fileName outputPath_;
 
+        //- Output geometry scaling after rotate/translate
+        scalar geometryScale_;
+
+        //- The centre of rotation (untranslate, translate)
+        point geometryCentre_;
+
+        //- Local coordinate system transformation
+        coordSystem::cartesian geometryTransform_;
+
+        //- Field level to remove (on output)
+        dictionary fieldLevel_;
+
+        //- Field scaling (on output)
+        dictionary fieldScale_;
+
 
     // Buffering
 
@@ -212,6 +228,13 @@ protected:
 
     // Helpers
 
+        template<class Type>
+        tmp<Field<Type>> adjustFieldTemplate
+        (
+            const word& fieldName,
+            const tmp<Field<Type>>& tfield
+        ) const;
+
         //- Repackage field into a UPtrList
         template<class Type>
         static UPtrList<const Field<Type>>
diff --git a/src/meshTools/coordSet/writers/common/coordSetWriterTemplates.C b/src/meshTools/coordSet/writers/common/coordSetWriterTemplates.C
index 84c5b0c336..28e8f08de5 100644
--- a/src/meshTools/coordSet/writers/common/coordSetWriterTemplates.C
+++ b/src/meshTools/coordSet/writers/common/coordSetWriterTemplates.C
@@ -25,8 +25,110 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "transformField.H"
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::coordSetWriter::adjustFieldTemplate
+(
+    const word& fieldName,
+    const tmp<Field<Type>>& tfield
+) const
+{
+    if (verbose_)
+    {
+        Info<< "Writing field " << fieldName;
+    }
+
+    tmp<Field<Type>> tadjusted;
+
+    // Output scaling for the variable, but not for integer types
+    // which are typically ids etc.
+    if (!std::is_integral<Type>::value)
+    {
+        scalar value;
+
+        // Remove *uniform* reference level
+        if
+        (
+            fieldLevel_.readIfPresent(fieldName, value, keyType::REGEX)
+         && !equal(value, 0)
+        )
+        {
+            // Could also detect brackets (...) and read accordingly
+            // or automatically scale by 1/sqrt(nComponents) instead ...
+
+            Type refLevel;
+            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+            {
+                setComponent(refLevel, cmpt) = value;
+            }
+
+            if (verbose_)
+            {
+                Info<< " [level " << refLevel << ']';
+            }
+
+            if (!tadjusted)
+            {
+                // Steal or clone
+                tadjusted.reset(tfield.ptr());
+            }
+
+            // Remove offset level
+            tadjusted.ref() -= refLevel;
+        }
+
+        // Apply scaling
+        if
+        (
+            fieldScale_.readIfPresent(fieldName, value, keyType::REGEX)
+         && !equal(value, 1)
+        )
+        {
+            if (verbose_)
+            {
+                Info<< " [scaling " << value << ']';
+            }
+
+            if (!tadjusted)
+            {
+                // Steal or clone
+                tadjusted.reset(tfield.ptr());
+            }
+
+            // Apply scaling
+            tadjusted.ref() *= value;
+        }
+
+        // Rotate fields (vector and non-spherical tensors)
+        if
+        (
+            (pTraits<Type>::rank != 0 && pTraits<Type>::nComponents > 1)
+         && geometryTransform_.valid()
+         && !geometryTransform_.R().is_identity()
+        )
+        {
+            if (!tadjusted)
+            {
+                // Steal or clone
+                tadjusted.reset(tfield.ptr());
+            }
+
+           Foam::transform
+           (
+               tadjusted.ref(),
+               geometryTransform_.R(),
+               tadjusted()
+           );
+        }
+    }
+
+    return (tadjusted ? tadjusted : tfield);
+}
+
+
 template<class Type>
 Foam::UPtrList<const Foam::Field<Type>>
 Foam::coordSetWriter::repackageFields(const Field<Type>& field)

From d3b3a5a41e9b59b0b5d2c7c8276106ff1fadf0ef Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Fri, 24 Nov 2023 16:02:40 +0000
Subject: [PATCH 2/3] ENH: Added new abaqus coordSet writer

Write coordSet(s) as Abaqus point fields

Example usage

    T
    {
        type        sets;
        setFormat   abaqus;
        fields      (T);
        sets
        {
            ...
        }
    }
    \endverbatim

    Optional format options
    \verbatim
    formatOptions
    {
        abaqus
        {
            format          ascii;

            // Optional entries

            // Custom header: $ entries are substituions
            header
            (
                "** OpenFOAM abaqus output"
                "** Project $FOAM_CASE"
                "** File $FILE_NAME"
                "** $FIELD_NAME Time t=$TIME"
            );

            // Write geometry in addition to field data
            writeGeometry   yes;

            // Null value when sample value is not found
            // Default is scalar::min
            nullValue       0;

            // Insert additional time sub-directory in the output path
            // - yes : postProcessing/<fo-name>/<time>/<file>
            // - no  : postProcessing/<fo-name>/<file>
            useTimeDir      no;

            // Available when 'useTimeDir' is 'no' to disambiguate file names

            // Time base for output file names:
            // - 'time'      : <base>.inp_<field>.<time>
            // - 'iteration' : <base>.inp_<field>.<iteration>
            timeBase        iteration;

            // Optional start counters when using timeBase iteration
            writeIndex
            (
                T 1
            );

            ...
        }
    }
---
 src/meshTools/Make/files                      |   1 +
 .../writers/abaqus/abaqusCoordSetWriter.C     | 409 ++++++++++++++++++
 .../writers/abaqus/abaqusCoordSetWriter.H     | 251 +++++++++++
 3 files changed, 661 insertions(+)
 create mode 100644 src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.C
 create mode 100644 src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.H

diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index bd0112da39..6218645011 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -116,6 +116,7 @@ coordSet/coordSet.C
 
 setWriters = coordSet/writers
 
+$(setWriters)/abaqus/abaqusCoordSetWriter.C
 $(setWriters)/common/coordSetWriter.C
 $(setWriters)/common/coordSetWriterBuffers.C
 $(setWriters)/common/coordSetWriterNew.C
diff --git a/src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.C b/src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.C
new file mode 100644
index 0000000000..4fc418c718
--- /dev/null
+++ b/src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.C
@@ -0,0 +1,409 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "abaqusCoordSetWriter.H"
+#include "coordSet.H"
+#include "IOmanip.H"
+#include "OFstream.H"
+#include "OSspecific.H"
+#include "stringOps.H"
+#include "coordSetWriterMethods.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace coordSetWriters
+{
+    defineTypeName(abaqusWriter);
+    addToRunTimeSelectionTable(coordSetWriter, abaqusWriter, word);
+    addToRunTimeSelectionTable(coordSetWriter, abaqusWriter, wordDict);
+}
+}
+
+const Foam::Enum<Foam::coordSetWriters::abaqusWriter::timeBase>
+Foam::coordSetWriters::abaqusWriter::timeBaseNames_
+({
+    { timeBase::time, "time" },
+    { timeBase::iter, "iteration" },
+});
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+template<class Type>
+static inline void putValue(Ostream& os, const Type& value, const int width)
+{
+    if (width) os << setw(width);
+    os << value;
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::string Foam::coordSetWriters::abaqusWriter::replaceUserEntries
+(
+    const string& str,
+    const dictionary& vars
+) const
+{
+    string result = str;
+
+    const bool allowEnv = true;
+    const bool allowEmpty = false;
+
+    stringOps::inplaceExpand(result, vars, allowEnv, allowEmpty);
+
+    return result;
+}
+
+
+void Foam::coordSetWriters::abaqusWriter::appendTimeName
+(
+    const word& fieldName,
+    fileName& fName
+) const
+{
+    if (useTimeDir())
+    {
+        return;
+    }
+
+    switch (timeBase_)
+    {
+        case timeBase::time:
+        {
+            fName.ext(timeName());
+            break;
+        }
+        case timeBase::iter:
+        {
+            fName.ext(Foam::name(writeIndex_[fieldName]));
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unhandled enumeration " << timeBaseNames_[timeBase_]
+                << ". Available options: " << timeBaseNames_.sortedToc()
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::coordSetWriters::abaqusWriter::abaqusWriter()
+:
+    coordSetWriter(),
+    outputHeader_(),
+    writeGeometry_(false),
+    nullValue_(pTraits<scalar>::min),
+    useLocalTimeDir_(coordSetWriter::useTimeDir()),
+    timeBase_(timeBase::time),
+    writeIndex_(0)
+{}
+
+
+Foam::coordSetWriters::abaqusWriter::abaqusWriter(const dictionary& options)
+:
+    coordSetWriter(options),
+    outputHeader_(),
+    writeGeometry_(false),
+    nullValue_(pTraits<scalar>::min),
+    useLocalTimeDir_(coordSetWriter::useTimeDir()),
+    timeBase_(timeBase::time),
+    writeIndex_(0)
+{
+    options.readIfPresent("header", outputHeader_);
+
+    options.readIfPresent("useTimeDir", useLocalTimeDir_);
+
+    if (!useLocalTimeDir_)
+    {
+        timeBaseNames_.readIfPresent("timeBase", options, timeBase_);
+        options.readIfPresent("writeIndex", writeIndex_);
+    }
+
+    options.readIfPresent("writeGeometry", writeGeometry_);
+
+    options.readIfPresent("nullValue", nullValue_);
+}
+
+
+Foam::coordSetWriters::abaqusWriter::abaqusWriter
+(
+    const coordSet& coords,
+    const fileName& outputPath,
+    const dictionary& options
+)
+:
+    abaqusWriter(options)
+{
+    open(coords, outputPath);
+}
+
+
+Foam::coordSetWriters::abaqusWriter::abaqusWriter
+(
+    const UPtrList<coordSet>& tracks,
+    const fileName& outputPath,
+    const dictionary& options
+)
+:
+    abaqusWriter(options)
+{
+    open(tracks, outputPath);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::coordSetWriters::abaqusWriter::~abaqusWriter()
+{
+    close();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::fileName Foam::coordSetWriters::abaqusWriter::path() const
+{
+    // 1) rootdir/<TIME>/setName.{inp}
+    // 2) rootdir/setName.{inp}
+
+    return getExpectedPath("inp");
+}
+
+
+void Foam::coordSetWriters::abaqusWriter::writeGeometry
+(
+    Ostream& os,
+    label nTracks
+) const
+{
+    if (!writeGeometry_ || coords_.empty())
+    {
+        return;
+    }
+
+    os  << "** Geometry" << nl
+        << "**" << nl
+        << "** Points" << nl
+        << "**" << nl;
+
+    // Write points
+    label globalPointi = 1;
+    for (const coordSet& coords : coords_)
+    {
+        for (const point& p : coords)
+        {
+            const point tp = p*geometryScale_;
+
+            os  << globalPointi << ", "
+                << tp[0] << ", " << tp[1] << ", " << tp[2]  << nl;
+
+            ++globalPointi;
+        }
+    }
+
+    if (nTracks)
+    {
+        WarningInFunction
+            << "Tracks not implemented for " << typeName << endl;
+    }
+
+    wroteGeom_ = true;
+}
+
+
+template<class Type>
+Foam::fileName Foam::coordSetWriters::abaqusWriter::writeTemplate
+(
+    const word& fieldName,
+    const Field<Type>& values
+)
+{
+    // useTimeDir(options.getOrDefault("useTimeDir", true));
+    useTimeDir(useLocalTimeDir_);
+
+    // Note - invalid samples are set to pTraits<Type>::max in sampledSetsImpl.C
+
+    checkOpen();
+    if (coords_.empty())
+    {
+        return fileName::null;
+    }
+
+    fileName outputFile = path();
+
+    if (!wroteGeom_)
+    {
+        if (!writeIndex_.insert(fieldName, 0))
+        {
+            ++writeIndex_[fieldName];
+        }
+
+        if (!isDir(outputFile.path()))
+        {
+            mkDir(outputFile.path());
+        }
+
+        if (writeGeometry_)
+        {
+            const word geomName("geometry");
+            if (!writeIndex_.insert(geomName, 0))
+            {
+                ++writeIndex_[geomName];
+            }
+
+            fileName geomFileName(outputFile.lessExt() + ".inp");
+
+            appendTimeName("geometry", geomFileName);
+
+            if (verbose_)
+            {
+                Info<< "Writing abaqus geometry to " << geomFileName << endl;
+            }
+
+            OFstream osGeom(geomFileName);
+
+            writeGeometry(osGeom, (useTracks_ ? coords_.size() : 0));
+        }
+
+        fileName fieldFileName(outputFile.lessExt() + ".inp_" + fieldName);
+        appendTimeName(fieldName, fieldFileName);
+
+        if (verbose_)
+        {
+            Info<< "Writing field data to " << fieldFileName << endl;
+        }
+
+        OFstream os(fieldFileName);
+
+        if (!outputHeader_.empty())
+        {
+            dictionary vars;
+            vars.add("TIME", timeName());
+            vars.add("FIELD_NAME", fieldName);
+            vars.add("FILE_NAME", fieldFileName);
+
+            for (const auto& s : outputHeader_)
+            {
+                os  << replaceUserEntries(s, vars).c_str() << nl;
+            }
+        }
+        else
+        {
+            os  << "** OpenFOAM " << fieldFileName.nameLessExt() << nl
+                << "** Project " << outputFile << nl
+                << "** Field=" << fieldName << " Time=" << timeName() << nl;
+        }
+
+        tmp<Field<Type>> tfield(values);
+        tfield = adjustFieldTemplate(fieldName, tfield);
+        const auto& field = tfield();
+
+        forAll(field, samplei)
+        {
+            os  << (samplei+1);
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
+            {
+                // Work-around to set null values - set to pTraits<Type>::max
+                // by default
+                scalar s = component(field[samplei], cmpt);
+                if (s > 0.5*pTraits<scalar>::max)
+                {
+                    s = nullValue_;
+                }
+
+                os  << ", " << s;
+            }
+            os  << nl;
+        }
+    }
+
+    return outputFile;
+}
+
+
+template<class Type>
+Foam::fileName Foam::coordSetWriters::abaqusWriter::writeTemplate
+(
+    const word& fieldName,
+    const List<Field<Type>>& fieldValues
+)
+{
+    // Track writing
+
+    checkOpen();
+    if (coords_.empty())
+    {
+        return fileName::null;
+    }
+
+    fileName outputFile = path();
+
+    if (!wroteGeom_)
+    {
+        if (verbose_)
+        {
+            Info<< "Writing abaqus geometry to " << outputFile << endl;
+        }
+
+        if (!isDir(outputFile.path()))
+        {
+            mkDir(outputFile.path());
+        }
+
+        OFstream osGeom(outputFile.lessExt() + "." + fieldName + ".inp");
+
+        osGeom
+            << "** Geometry" << nl
+            << "**" << nl
+            << "** Points" << nl
+            << "**" << nl;
+
+        writeGeometry(osGeom, coords_.size());
+    }
+
+    return outputFile;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Field writing methods
+defineCoordSetWriterWriteFields(Foam::coordSetWriters::abaqusWriter);
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.H b/src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.H
new file mode 100644
index 0000000000..4e35084ea9
--- /dev/null
+++ b/src/meshTools/coordSet/writers/abaqus/abaqusCoordSetWriter.H
@@ -0,0 +1,251 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::coordSetWriters::abaqusWriter
+
+Description
+    Write coordSet(s) as Abaqus point fields
+
+    Example usage
+    \verbatim
+    T
+    {
+        type        sets;
+        setFormat   abaqus;
+        fields      (T);
+        sets
+        {
+            ...
+        }
+    }
+    \endverbatim
+
+    Optional format options
+    \verbatim
+    formatOptions
+    {
+        abaqus
+        {
+            format          ascii;
+
+            // Optional entries
+
+            // Custom header: $ entries are substituions
+            header
+            (
+                "** OpenFOAM abaqus output"
+                "** Project $FOAM_CASE"
+                "** File $FILE_NAME"
+                "** $FIELD_NAME Time t=$TIME"
+            );
+
+            // Write geometry in addition to field data
+            writeGeometry   yes;
+
+            // Null value when sample value is not found
+            // Default is scalar::min
+            nullValue       0;
+
+            // Insert additional time sub-directory in the output path
+            // - yes : postProcessing/<fo-name>/<time>/<file>
+            // - no  : postProcessing/<fo-name>/<file>
+            useTimeDir      no;
+
+            // Available when 'useTimeDir' is 'no' to disambiguate file names
+
+            // Time base for output file names:
+            // - 'time'      : <base>.inp_<field>.<time>
+            // - 'iteration' : <base>.inp_<field>.<iteration>
+            timeBase        iteration;
+
+            // Optional start counters when using timeBase iteration
+            writeIndex
+            (
+                T 1
+            );
+
+
+            ...
+        }
+    }
+    \endverbatim
+
+SourceFiles
+    abaqusCoordSetWriter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_coordSetWriters_abaqusWriter_H
+#define Foam_coordSetWriters_abaqusWriter_H
+
+#include "coordSetWriter.H"
+#include "Enum.H"
+#include "HashTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace coordSetWriters
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class abaqusWriter Declaration
+\*---------------------------------------------------------------------------*/
+
+class abaqusWriter
+:
+    public coordSetWriter
+{
+public:
+
+    // Public enumerations
+
+        //- Enumeration for time base
+        enum class timeBase { time, iter };
+
+
+private:
+
+    // Private Data
+
+        //- Time base names
+        static const Enum<timeBase> timeBaseNames_;
+
+        //- Optional user-defined header
+        List<string> outputHeader_;
+
+        //- User flag to write the geometry
+        bool writeGeometry_;
+
+        //- Null value; default = scalar::min
+        scalar nullValue_;
+
+        //- Optional override of localTimeDir
+        bool useLocalTimeDir_;
+
+        //- Optional time base when useLocalTimeDir_ = false
+        timeBase timeBase_;
+
+        //- Write index
+        HashTable<label> writeIndex_;
+
+
+    // Private Member Functions
+
+        //- Helper to replace $WORD entries in str
+        string replaceUserEntries
+        (
+            const string& str,
+            const dictionary& vars
+        ) const;
+
+        //- Write the formatted keyword to the output stream
+        Ostream& writeKeyword(Ostream& os, const word& keyword) const;
+
+        //- Append time name when useLocalTimeDir_ = false
+        void appendTimeName(const word& fieldName, fileName& fName) const;
+
+
+        //- Templated write operation
+        template<class Type>
+        fileName writeTemplate
+        (
+            const word& fieldName,      //!< Name of field
+            const Field<Type>& values   //!< Local field values to write
+        );
+
+        //- Templated write operation
+        template<class Type>
+        fileName writeTemplate
+        (
+            const word& fieldName,
+            const List<Field<Type>>& fieldValues
+        );
+
+        //- Write geometry to file.
+        void writeGeometry(Ostream& os, label nTracks) const;
+
+
+public:
+
+    //- Runtime type information (no debug)
+    TypeNameNoDebug("abaqus");
+
+
+    // Constructors
+
+        //- Default construct
+        abaqusWriter();
+
+        //- Default construct with specified options
+        explicit abaqusWriter(const dictionary& options);
+
+        //- Construct from components
+        abaqusWriter
+        (
+            const coordSet& coords,
+            const fileName& outputPath,
+            const dictionary& options = dictionary()
+        );
+
+        //- Construct from components
+        abaqusWriter
+        (
+            const UPtrList<coordSet>& tracks,
+            const fileName& outputPath,
+            const dictionary& options = dictionary()
+        );
+
+
+    //- Destructor. Calls close()
+    virtual ~abaqusWriter();
+
+
+    // Member Functions
+
+        //- Characteristic output file name - information only
+        virtual fileName path() const;  // override
+
+        declareCoordSetWriterWriteMethod(label);
+        declareCoordSetWriterWriteMethod(scalar);
+        declareCoordSetWriterWriteMethod(vector);
+        declareCoordSetWriterWriteMethod(sphericalTensor);
+        declareCoordSetWriterWriteMethod(symmTensor);
+        declareCoordSetWriterWriteMethod(tensor);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace coordSetWriters
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From a1819a5411aaf3f180622802b8f91fe1715d70f0 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Mon, 12 Jun 2023 09:16:57 +0100
Subject: [PATCH 3/3] ENH: Added new abaqus sampledSet

Generates sample positions from points specified in a file as Abaqus mesh
points.

Example usage:

    sets
    {
        cone25 // user-specified set name
        {
            type        abaqusMesh;
            file        "abaqusMesh.inp";

            // Optional entries

            // Scale, e.g. mm to m
            scale       0.001;

            // Search distance when the sample point is not located in a cell
            maxDist     0.25;

            ...
        }
    }
---
 src/sampling/Make/files                       |   1 +
 .../sampledSet/abaqus/abaqusMeshSet.C         | 346 ++++++++++++++++++
 .../sampledSet/abaqus/abaqusMeshSet.H         | 150 ++++++++
 3 files changed, 497 insertions(+)
 create mode 100644 src/sampling/sampledSet/abaqus/abaqusMeshSet.C
 create mode 100644 src/sampling/sampledSet/abaqus/abaqusMeshSet.H

diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 16e134d613..afa6baf0a3 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -1,6 +1,7 @@
 probes/probes.C
 probes/patchProbes.C
 
+sampledSet/abaqus/abaqusMeshSet.C
 sampledSet/circle/circleSet.C
 sampledSet/cloud/cloudSet.C
 sampledSet/patchCloud/patchCloudSet.C
diff --git a/src/sampling/sampledSet/abaqus/abaqusMeshSet.C b/src/sampling/sampledSet/abaqus/abaqusMeshSet.C
new file mode 100644
index 0000000000..a4495114fa
--- /dev/null
+++ b/src/sampling/sampledSet/abaqus/abaqusMeshSet.C
@@ -0,0 +1,346 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "abaqusMeshSet.H"
+#include "stringOps.H"
+#include "Fstream.H"
+#include "SpanStream.H"
+#include "meshSearch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(abaqusMeshSet, 0);
+    addToRunTimeSelectionTable(sampledSet, abaqusMeshSet, word);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::abaqusMeshSet::readCoord(ISstream& is, vector& coord) const
+{
+    string buffer;
+
+    do
+    {
+        //buffer.clear();
+
+        is.getLine(buffer);
+
+        const auto elems = buffer.find("*ELEMENT");
+
+        if (elems != std::string::npos)
+        {
+            buffer.clear();
+            break;
+        }
+
+        // Trim out any '*' comments
+        const auto pos = buffer.find('*');
+        if (pos != std::string::npos)
+        {
+            buffer.erase(pos);
+        }
+        stringOps::inplaceTrimRight(buffer);
+    }
+    while (buffer.empty() && is.good());
+
+    if (buffer.empty())
+    {
+        return false;
+    }
+
+    const auto strings = stringOps::split(buffer, ',');
+
+    if (strings.size() != 4)
+    {
+        FatalErrorInFunction
+            << "Read error: expected format int, float, float, float"
+            << " but read buffer " << buffer
+            << exit(FatalError);
+    }
+
+    for (int i = 0; i <= 2; ++i)
+    {
+        // Swallow i=0 for node label
+        const auto& s = strings[i+1].str();
+        ISpanStream buf(s.data(), s.length());
+        buf >> coord[i];
+    }
+
+    coord *= scale_;
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::abaqusMeshSet::calcSamples
+(
+    DynamicList<point>& samplingPts,
+    DynamicList<label>& samplingCells,
+    DynamicList<label>& samplingFaces,
+    DynamicList<label>& samplingSegments,
+    DynamicList<scalar>& samplingCurveDist
+) const
+{
+    DebugInfo
+        << "abaqusMeshSet : sampling " << sampleCoords_.size() << " points"
+        << endl;
+
+    List<bool> found(sampleCoords_.size(), false);
+
+    forAll(sampleCoords_, samplei)
+    {
+        const vector& pt = sampleCoords_[samplei];
+
+        label celli = searchEngine().findCell(pt);
+        if (celli != -1)
+        {
+            found[samplei] = true;
+
+            samplingPts.append(pt);
+            samplingCells.append(celli);
+            samplingFaces.append(-1);
+            samplingSegments.append(0);
+            samplingCurveDist.append(1.0*samplei);
+        }
+    }
+
+    Pstream::listCombineAllGather(found, orEqOp<bool>());
+
+    DynamicList<label> lost;
+    forAll(found, samplei)
+    {
+        if (!found[samplei]) lost.append(samplei);
+    }
+
+    const label nFound = sampleCoords_.size() - lost.size();
+    label nOutOfBounds = 0;
+
+    if (Pstream::parRun())
+    {
+        typedef Tuple2<label, Tuple2<label, scalar>> procCellDist;
+        List<procCellDist> pcd(lost.size(), procCellDist(-1, {-1, GREAT}));
+        forAll(pcd, i)
+        {
+            const label samplei = lost[i];
+            const vector& pt0 = sampleCoords_[samplei];
+            const label celli = searchEngine().findNearestCell(pt0);
+            const vector& pt1 = mesh().cellCentres()[celli];
+            const scalar distSqr = magSqr(pt1 - pt0);
+            pcd[i] = procCellDist(Pstream::myProcNo(), {celli, distSqr});
+        }
+
+        Pstream::listCombineAllGather
+        (
+            pcd,
+            [](procCellDist& x, const procCellDist& y)
+            {
+                // Select item with smallest distance
+                if (y.second().second() < x.second().second())
+                {
+                    x = y;
+                }
+            }
+        );
+
+        forAll(pcd, i)
+        {
+            const label samplei = lost[i];
+
+            if (pcd[i].second().second() < maxDistSqr_)
+            {
+                if (pcd[i].first() == Pstream::myProcNo())
+                {
+                    const label celli = pcd[i].second().first();
+                    const vector& pt1 = mesh().cellCentres()[celli];
+
+                    samplingPts.append(pt1);
+                    samplingCells.append(celli);
+                    samplingFaces.append(-1);
+                    samplingSegments.append(0);
+                    samplingCurveDist.append(1.0*samplei);
+                }
+            }
+            else
+            {
+                // Insert points that have not been found as null points
+                if (Pstream::master())
+                {
+                    samplingPts.append(sampleCoords_[samplei]);
+                    samplingCells.append(-1);
+                    samplingFaces.append(-1);
+                    samplingSegments.append(0);
+                    samplingCurveDist.append(1.0*samplei);
+                    ++nOutOfBounds;
+                }
+            }
+        }
+    }
+    else
+    {
+        // Serial running
+
+        forAll(lost, i)
+        {
+            const label samplei = lost[i];
+            const vector& pt0 = sampleCoords_[samplei];
+            const label celli = searchEngine().findNearestCell(pt0);
+            const vector& pt1 = mesh().cellCentres()[celli];
+            const scalar distSqr = magSqr(pt1 - pt0);
+
+            if (distSqr < maxDistSqr_)
+            {
+                samplingPts.append(pt1);
+                samplingCells.append(celli);
+                samplingFaces.append(-1);
+                samplingSegments.append(0);
+                samplingCurveDist.append(1.0*samplei);
+            }
+            else
+            {
+                // Insert points that have not been found as null points
+                samplingPts.append(sampleCoords_[samplei]);
+                samplingCells.append(-1);
+                samplingFaces.append(-1);
+                samplingSegments.append(0);
+                samplingCurveDist.append(1.0*samplei);
+                ++nOutOfBounds;
+            }
+        }
+    }
+
+    DebugInFunction
+        << "Sample size       : " << sampleCoords_.size() << nl
+        << "Lost samples      : " << nOutOfBounds << nl
+        << "Recovered samples : "
+        << (sampleCoords_.size() - nOutOfBounds - nFound) << nl
+        << endl;
+
+
+    if (nOutOfBounds)
+    {
+        WarningInFunction
+            << "Identified " << nOutOfBounds << " out-of-bounds points"
+            << endl;
+    }
+}
+
+
+void Foam::abaqusMeshSet::genSamples()
+{
+    // Storage for sample points
+    DynamicList<point> samplingPts;
+    DynamicList<label> samplingCells;
+    DynamicList<label> samplingFaces;
+    DynamicList<label> samplingSegments;
+    DynamicList<scalar> samplingCurveDist;
+
+    calcSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+
+    samplingPts.shrink();
+    samplingCells.shrink();
+    samplingFaces.shrink();
+    samplingSegments.shrink();
+    samplingCurveDist.shrink();
+
+    setSamples
+    (
+        std::move(samplingPts),
+        std::move(samplingCells),
+        std::move(samplingFaces),
+        std::move(samplingSegments),
+        std::move(samplingCurveDist)
+    );
+
+    if (debug > 1)
+    {
+        write(Info);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::abaqusMeshSet::abaqusMeshSet
+(
+    const word& name,
+    const polyMesh& mesh,
+    const meshSearch& searchEngine,
+    const dictionary& dict
+)
+:
+    sampledSet(name, mesh, searchEngine, dict),
+    scale_(dict.getOrDefault<scalar>("scale", 1)),
+    sampleCoords_(),
+    maxDistSqr_(sqr(dict.getOrDefault<scalar>("maxDist", 0)))
+{
+    if (Pstream::master())
+    {
+        const fileName inputFile(dict.get<fileName>("file").expand());
+        IFstream pointsFile(inputFile);
+
+        if (!pointsFile.good())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Unable to find file " << pointsFile.name()
+                << abort(FatalIOError);
+        }
+
+        // Read the points file
+        DynamicList<point> coords;
+        vector c;
+        while (readCoord(pointsFile, c))
+        {
+            coords.append(c);
+        }
+
+        sampleCoords_.transfer(coords);
+    }
+
+    Pstream::broadcast(sampleCoords_);
+
+    DebugInfo
+        << "Number of sample points: " << sampleCoords_.size() << nl
+        << "Sample points bounds: " << boundBox(sampleCoords_) << endl;
+
+    genSamples();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSet/abaqus/abaqusMeshSet.H b/src/sampling/sampledSet/abaqus/abaqusMeshSet.H
new file mode 100644
index 0000000000..78e8efdcf1
--- /dev/null
+++ b/src/sampling/sampledSet/abaqus/abaqusMeshSet.H
@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::abaqusMeshSet
+
+Description
+    Generates sample positions from points specified in a file as Abaqus mesh
+    points.
+
+    Example usage:
+
+    \verbatim
+    sets
+    {
+        cone25 // user-specified set name
+        {
+            type        abaqusMesh;
+            file        "abaqusMesh.inp";
+
+            // Optional entries
+
+            // Scale, e.g. mm to m
+            scale       0.001;
+
+            // Search distance when the sample point is not located in a cell
+            maxDist     0.25;
+
+            ...
+        }
+    }
+    \endverbatim
+
+    For a dictionary specification:
+    \table
+        Property | Description                             | Required | Default
+        type     | abaqusMesh                              | yes      |
+        file     | Path to Abaqus file                     | yes      |
+        scale    | scale input point positions             | no       | 1
+        maxDist  | Search distance for sample points       | no       | 1
+    \endtable
+
+SourceFiles
+    abaqusMeshSet.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_abaqusMeshSet_H
+#define Foam_abaqusMeshSet_H
+
+#include "sampledSet.H"
+#include "vectorList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class abaqusMeshSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class abaqusMeshSet
+:
+    public sampledSet
+{
+    // Private data
+
+        //- Scale
+        const scalar scale_;
+
+        //- Sample coordinates
+        vectorList sampleCoords_;
+
+        //- Maximum search distance-squared
+        const scalar maxDistSqr_;
+
+
+    // Private Member Functions
+
+        bool readCoord(ISstream& is, vector& coord) const;
+
+        //- Samples all points in sampleCoords.
+        void calcSamples
+        (
+            DynamicList<point>& samplingPts,
+            DynamicList<label>& samplingCells,
+            DynamicList<label>& samplingFaces,
+            DynamicList<label>& samplingSegments,
+            DynamicList<scalar>& samplingCurveDist
+        ) const;
+
+        //- Uses calcSamples to obtain samples. Copies them into *this.
+        void genSamples();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("abaqusMesh");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        abaqusMeshSet
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const meshSearch& searchEngine,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~abaqusMeshSet() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //