ENH: initial support for abaqus surface sampled output (#1600)

- supports geometry and field-specific scaling, separate geometry and
  fields. Beta-feature for suppressing geometry output entirely.

src/surfMesh/Make/files

@@ -63,6 +63,7 @@ writers = writers
 
 $(writers)/surfaceWriter.C
 $(writers)/caching/surfaceWriterCaching.C
+$(writers)/abaqus/abaqusSurfaceWriter.C
 $(writers)/boundaryData/boundaryDataSurfaceWriter.C
 $(writers)/ensight/ensightSurfaceWriter.C
 $(writers)/foam/foamSurfaceWriter.C

src/surfMesh/writers/abaqus/abaqusSurfaceWriter.C

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 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 "abaqusSurfaceWriter.H"
+#include "ABAQUSCore.H"
+#include "IOmanip.H"
+#include "ListOps.H"
+#include "OSspecific.H"
+#include "surfaceWriterMethods.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceWriters
+{
+    defineTypeName(abaqusWriter);
+    addToRunTimeSelectionTable(surfaceWriter, abaqusWriter, word);
+    addToRunTimeSelectionTable(surfaceWriter, abaqusWriter, wordDict);
+}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Field writing implementation
+#include "abaqusSurfaceWriterImpl.C"
+
+// Field writing methods
+defineSurfaceWriterWriteFields(Foam::surfaceWriters::abaqusWriter);
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Write connectivity as CSV list
+inline static void writeConnectivity
+(
+    Ostream& os,
+    const label elemId,
+    const labelUList& elem
+)
+{
+    os  << "  " << elemId;
+
+    for (const label vert : elem)
+    {
+        os << ", " << (vert + 1);
+    }
+
+    os << nl;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::surfaceWriters::abaqusWriter::writeFace
+(
+    Ostream& os,
+    const labelUList& f,
+    const label elemId,
+    const label propId,
+    bool header
+) const
+{
+    // Only called with 3 or 4 points!
+
+    if (header)
+    {
+        os  << "*ELEMENT, TYPE=S" << f.size();
+
+        if (propId >= 0)
+        {
+            os  << ", ELSET=_" << propId;
+        }
+
+        os  << nl;
+    }
+
+    writeConnectivity(os, elemId, f);
+}
+
+
+void Foam::surfaceWriters::abaqusWriter::writeGeometry
+(
+    Ostream& os,
+    const meshedSurf& surf,
+    labelList& decompOffsets,
+    DynamicList<face>& decompFaces
+) const
+{
+    const pointField& points = surf.points();
+    const faceList&    faces = surf.faces();
+    const labelList&   zones = surf.zoneIds();
+    const labelList& elemIds = surf.faceIds();
+
+    // Possible to use faceIds?
+    bool useOrigFaceIds =
+    (
+        elemIds.size() == faces.size()
+     && !ListOps::found(elemIds, lessOp1<label>(0))
+    );
+
+    if (useOrigFaceIds)
+    {
+        // Not possible with on-the-fly face decomposition
+        for (const auto& f : faces)
+        {
+            if (f.size() > 4)
+            {
+                useOrigFaceIds = false;
+                break;
+            }
+        }
+    }
+
+
+    os  << "** Geometry" << nl;
+
+    os  << nl
+        << "**" << nl
+        << "** Points" << nl
+        << "**" << nl;
+
+    fileFormats::ABAQUSCore::writePoints(os, points, geometryScale_);
+
+
+    // Write faces, with on-the-fly decomposition (triangulation)
+    decompOffsets.resize(faces.size()+1);
+    decompFaces.clear();
+
+    decompOffsets[0] = 0; // The first offset is always zero
+
+    os  << "**" << nl
+        << "** Faces" << nl
+        << "**" << nl;
+
+    // Simple tracking for change of element type/set
+    labelPair prevOutput(-1, -1);
+
+    label elemId = 0;  // The element-id
+    forAll(faces, facei)
+    {
+        const face& f = faces[facei];
+
+        if (useOrigFaceIds)
+        {
+            // When available and not decomposed
+            elemId = elemIds[facei];
+        }
+
+        // 1-offset for PID
+        const label propId = 1 + (facei < zones.size() ? zones[facei] : 0);
+
+        const label n = f.size();
+
+        bool header =
+            (prevOutput.first() != n || prevOutput.second() != propId);
+
+        if (header)
+        {
+            // Update values
+            prevOutput.first() = n;
+            prevOutput.second() = propId;
+        }
+
+        if (n == 3 || n == 4)
+        {
+            writeFace(os, f, ++elemId, propId, header);
+        }
+        else
+        {
+            // Decompose into tris
+            prevOutput.first() = 3;
+
+            f.triangles(points, decompFaces);
+
+            for
+            (
+                label decompi = decompOffsets[facei];
+                decompi < decompFaces.size();
+                ++decompi
+            )
+            {
+                writeFace
+                (
+                    os,
+                    decompFaces[decompi],
+                    ++elemId,
+                    propId,
+                    header
+                );
+
+                header = false;
+            }
+        }
+
+        // The end offset, which is the next begin offset
+        decompOffsets[facei+1] = decompFaces.size();
+    }
+
+    os  << "**" << nl
+        << "**" << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter()
+:
+    surfaceWriter(),
+    geometryScale_(1),
+    fieldScale_(),
+    noGeometry_(false),
+    outputLayout_(outputLayoutType::BY_FIELD)
+{}
+
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter
+(
+    const dictionary& options
+)
+:
+    surfaceWriter(options),
+    geometryScale_(options.getOrDefault<scalar>("scale", 1)),
+    fieldScale_(options.subOrEmptyDict("fieldScale")),
+    noGeometry_(options.getOrDefault("noGeometry", false)),
+    outputLayout_(outputLayoutType::BY_FIELD)
+{}
+
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter
+(
+    const meshedSurf& surf,
+    const fileName& outputPath,
+    bool parallel,
+    const dictionary& options
+)
+:
+    abaqusWriter(options)
+{
+    open(surf, outputPath, parallel);
+}
+
+
+Foam::surfaceWriters::abaqusWriter::abaqusWriter
+(
+    const pointField& points,
+    const faceList& faces,
+    const fileName& outputPath,
+    bool parallel,
+    const dictionary& options
+)
+:
+    abaqusWriter(options)
+{
+    open(points, faces, outputPath, parallel);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::fileName Foam::surfaceWriters::abaqusWriter::write()
+{
+    checkOpen();
+
+    // Geometry:
+    // 1) rootdir/<TIME>/surfaceName.abq
+    // 2) rootdir/geometry/surfaceName_<TIME>.abq
+
+    fileName outputFile;
+
+    switch (outputLayout_)
+    {
+        case outputLayoutType::BY_TIME:
+        {
+            outputFile = outputPath_;
+            if (useTimeDir() && !timeName().empty())
+            {
+                // Splice in time-directory
+                outputFile =
+                    outputPath_.path() / timeName() / outputPath_.name();
+            }
+            break;
+        }
+        case outputLayoutType::BY_FIELD:
+        {
+            outputFile = outputPath_ / "geometry" / outputPath_.name();
+            if (!timeName().empty())
+            {
+                // Append time information to file name
+                outputFile += '_' + timeName();
+            }
+            break;
+        }
+    }
+    outputFile.ext("abq");
+
+    if (verbose_)
+    {
+        Info<< "Writing abaqus geometry to " << outputFile << endl;
+    }
+
+
+    const meshedSurf& surf = surface();
+
+    if (Pstream::master() || !parallel_)
+    {
+        if (!isDir(outputFile.path()))
+        {
+            mkDir(outputFile.path());
+        }
+
+        OFstream os(outputFile);
+
+        labelList decompOffsets;
+        DynamicList<face> decompFaces;
+
+        writeGeometry(os, surf, decompOffsets, decompFaces);
+    }
+
+    wroteGeom_ = true;
+    return outputFile;
+}
+
+
+// ************************************************************************* //

src/surfMesh/writers/abaqus/abaqusSurfaceWriter.H

@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 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::surfaceWriters::abaqusWriter
+
+Description
+    A surface writer for the ABAQUS file format - both surface mesh and fields
+
+    The formatOptions for abaqus:
+    \table
+        Property   | Description                            | Required | Default
+        scale      | output geometry scaling                | no  | 1
+        fieldScale | output field scaling (dictionary)      | no  | empty
+        noGeometry | Suppress geometry output (beta feature) | no  | false
+    \endtable
+
+    For example,
+    \verbatim
+    formatOptions
+    {
+        abaqus
+        {
+            scale   1000;     // [m] -> [mm]
+            fieldScale
+            {
+               "p.*"   0.01;  // [Pa] -> [mbar]
+            }
+        }
+    }
+    \endverbatim
+
+    \heading Output file locations
+
+    The \c rootdir normally corresponds to something like
+    \c postProcessing/\<name\>
+
+    \subheading Geometry
+    \verbatim
+    rootdir
+    `-- surfaceName0
+        `-- geometry
+            `-- surfaceName0_<time>.abq
+    \endverbatim
+
+    \subheading Fields
+    \verbatim
+    rootdir
+    `-- surfaceName0
+        `-- field0
+        |   `-- surfaceName0_<time>.{inp}
+        `-- field1
+            `-- surfaceName0_<time>.{inp}
+    \endverbatim
+
+Note
+    Output variable scaling does not apply to integer types such as Ids.
+
+SourceFiles
+    abaqusSurfaceWriter.C
+    abaqusSurfaceWriterImpl.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef abaqusSurfaceWriter_H
+#define abaqusSurfaceWriter_H
+
+#include "surfaceWriter.H"
+#include "ABAQUSCore.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceWriters
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class abaqusWriter Declaration
+\*---------------------------------------------------------------------------*/
+
+class abaqusWriter
+:
+    public surfaceWriter
+{
+    // Private Data
+
+        //- Output geometry scaling
+        const scalar geometryScale_;
+
+        //- Output field scaling
+        const dictionary fieldScale_;
+
+        //- BETA feature
+        bool noGeometry_;
+
+        //- Output directory layouts
+        enum class outputLayoutType
+        {
+            BY_TIME = 0,
+            BY_FIELD,
+        };
+
+        //- Output directory layout
+        outputLayoutType outputLayout_;
+
+
+    // Private Member Functions
+
+        //- Write face element (tri/quad)
+        void writeFace
+        (
+            Ostream& os,
+            const labelUList& f,
+            const label elemId,     //!< 1-based Element Id
+            const label propId,     //!< 1-based Property Id
+            const bool header = true
+        ) const;
+
+        //- Write the surface mesh geometry, tracking face decomposition
+        //
+        //  \param decompOffsets  begin/end offsets (size+1) into decompFaces
+        //  \param decompFaces    Non tri/quad decomposed into triangles
+        void writeGeometry
+        (
+            Ostream& os,
+            const meshedSurf& surf,
+            labelList& decompOffsets,
+            DynamicList<face>& decompFaces
+        ) const;
+
+        //- Write a face-based value
+        template<class Type>
+        Ostream& writeFaceValue
+        (
+            Ostream& os,
+            const Type& value,
+            const label elemId      //!< 0-based Element Id
+        ) const;
+
+        //- Templated write operation
+        template<class Type>
+        fileName writeTemplate
+        (
+            const word& fieldName,          //!< Name of field
+            const Field<Type>& localValues  //!< Local field values to write
+        );
+
+
+public:
+
+    //- Declare type-name, virtual type (with debug switch)
+    TypeNameNoDebug("abaqus");
+
+
+    // Constructors
+
+        //- Default construct
+        abaqusWriter();
+
+        //- Construct with some output options
+        explicit abaqusWriter(const dictionary& options);
+
+        //- Construct from components
+        abaqusWriter
+        (
+            const meshedSurf& surf,
+            const fileName& outputPath,
+            bool parallel = Pstream::parRun(),
+            const dictionary& options = dictionary()
+        );
+
+        //- Construct from components
+        abaqusWriter
+        (
+            const pointField& points,
+            const faceList& faces,
+            const fileName& outputPath,
+            bool parallel = Pstream::parRun(),
+            const dictionary& options = dictionary()
+        );
+
+
+    //- Destructor
+    virtual ~abaqusWriter() = default;
+
+
+    // Member Functions
+
+        //- Format uses faceIds as part of its output
+        virtual bool usesFaceIds() const // override
+        {
+            return true;
+        }
+
+        //- Write surface geometry to file.
+        virtual fileName write(); // override
+
+        declareSurfaceWriterWriteMethod(label);
+        declareSurfaceWriterWriteMethod(scalar);
+        declareSurfaceWriterWriteMethod(vector);
+        declareSurfaceWriterWriteMethod(sphericalTensor);
+        declareSurfaceWriterWriteMethod(symmTensor);
+        declareSurfaceWriterWriteMethod(tensor);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceWriters
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/surfMesh/writers/abaqus/abaqusSurfaceWriterImpl.C

@@ -0,0 +1,308 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 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 "OFstream.H"
+#include "IOmanip.H"
+#include "OSspecific.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::Ostream& Foam::surfaceWriters::abaqusWriter::writeFaceValue
+(
+    Ostream& os,
+    const Type& value,
+    const label elemId
+) const
+{
+    if (fileFormats::ABAQUSCore::isEncodedSolidId(elemId))
+    {
+        const label solidId =
+            fileFormats::ABAQUSCore::decodeSolidElementId(elemId);
+
+        const label sideNum =
+            fileFormats::ABAQUSCore::decodeSolidSideNum(elemId);
+
+        // Element-id: 0-based to 1-based
+        // Side-id: always 1-based
+        os  << (solidId + 1) << ", P" << sideNum;
+    }
+    else
+    {
+        // Element-id: 0-based to 1-based
+        os  << (elemId + 1) << ", P";
+    }
+
+    os << ", ";
+    if (pTraits<Type>::nComponents == 1)
+    {
+        os << value;
+    }
+    else
+    {
+        os << mag(value);
+    }
+    os  << nl;
+
+    return os;
+}
+
+
+template<class Type>
+Foam::fileName Foam::surfaceWriters::abaqusWriter::writeTemplate
+(
+    const word& fieldName,
+    const Field<Type>& localValues
+)
+{
+    checkOpen();
+
+    // Field:
+    // 1) rootdir/<TIME>/<field>_surfaceName.inp
+    // 2) rootdir/<field>/surfaceName_<TIME>.inp
+
+    fileName outputFile;
+
+    switch (outputLayout_)
+    {
+        case outputLayoutType::BY_TIME:
+        {
+            outputFile = outputPath_;
+            if (useTimeDir() && !timeName().empty())
+            {
+                // Splice in time-directory
+                outputFile =
+                    outputPath_.path() / timeName() / outputPath_.name();
+            }
+
+            // Append <field>_surfaceName
+            outputFile /= fieldName + '_' + outputPath_.name();
+            break;
+        }
+        case outputLayoutType::BY_FIELD:
+        {
+            outputFile = outputPath_ / fieldName / outputPath_.name();
+            if (!timeName().empty())
+            {
+                // Append time information to file name
+                outputFile += '_' + timeName();
+            }
+            break;
+        }
+    }
+    outputFile.ext("inp");
+
+
+    // Output scaling for the variable, but not for integer types.
+    // could also solve with clever templating
+
+    const scalar varScale =
+    (
+        std::is_integral<Type>::value
+      ? scalar(1)
+      : fieldScale_.getOrDefault<scalar>(fieldName, 1)
+    );
+
+    if (verbose_)
+    {
+        Info<< "Writing field " << fieldName;
+        if (!equal(varScale, 1))
+        {
+            Info<< " (scaling " << varScale << ')';
+        }
+        Info<< " to " << outputFile << endl;
+    }
+
+
+    // Implicit geometry merge()
+    tmp<Field<Type>> tfield = mergeField(localValues) * varScale;
+
+    const meshedSurf& surf = surface();
+
+    if (Pstream::master() || !parallel_)
+    {
+        const auto& values = tfield();
+
+        if (!isDir(outputFile.path()))
+        {
+            mkDir(outputFile.path());
+        }
+
+        // const scalar timeValue(0);
+
+
+        // Additional bookkeeping for decomposing non tri/quad
+        labelList decompOffsets;
+        DynamicList<face> decompFaces;
+
+
+        OFstream os(outputFile);
+
+        if (noGeometry_ || wroteGeom_)
+        {
+            // Geometry already written (or suppressed)
+            // - still need decomposition information
+
+            fileFormats::ABAQUSCore::faceDecomposition
+            (
+                surf.points(),
+                surf.faces(),
+                decompOffsets,
+                decompFaces
+            );
+        }
+        else
+        {
+            // Write geometry (separate file)
+
+            OFstream osGeom(outputFile.lessExt().ext("abq"));
+            writeGeometry(osGeom, surf, decompOffsets, decompFaces);
+        }
+
+        // Write field
+        // *DLOAD
+        // Element-based: // elemId, P, 1000
+        // Surface-based: // elemId, P4, 1000
+
+        os  << "**" << nl
+            << "** field = " << fieldName << nl
+            << "** type = " << pTraits<Type>::typeName << nl;
+
+        if (useTimeDir() && !timeName().empty())
+        {
+            os  << "** time = " << timeName() << nl;
+        }
+
+        os  << "**" << nl
+            << "*DLOAD" << nl;
+
+
+        // Regular (undecomposed) faces
+        const faceList&    faces = surf.faces();
+        const labelList& elemIds = surf.faceIds();
+
+        // Possible to use faceIds?
+        const bool useOrigFaceIds =
+        (
+            elemIds.size() == faces.size()
+         && decompFaces.empty()
+        );
+
+
+        label elemId = 0;
+
+        if (this->isPointData())
+        {
+            forAll(faces, facei)
+            {
+                if (useOrigFaceIds)
+                {
+                    // When available and not decomposed
+                    elemId = elemIds[facei];
+                }
+
+                const label beginElemId = elemId;
+
+                // Any face decomposition
+                for
+                (
+                    label decompi = decompOffsets[facei];
+                    decompi < decompOffsets[facei+1];
+                    ++decompi
+                )
+                {
+                    const face& f = decompFaces[decompi];
+
+                    Type v = Zero;
+                    for (const label verti : f)
+                    {
+                        v += values[verti];
+                    }
+                    v /= f.size();
+
+                    writeFaceValue(os, v, elemId);  // 0-based
+                    ++elemId;
+                }
+
+
+                // Face not decomposed
+                if (beginElemId == elemId)
+                {
+                    const face& f = faces[facei];
+
+                    Type v = Zero;
+                    for (const label verti : f)
+                    {
+                        v += values[verti];
+                    }
+                    v /= f.size();
+
+                    writeFaceValue(os, v, elemId);  // 0-based
+                    ++elemId;
+                }
+            }
+        }
+        else
+        {
+            auto valIter = values.cbegin();
+
+            forAll(faces, facei)
+            {
+                if (useOrigFaceIds)
+                {
+                    // When available and not decomposed
+                    elemId = elemIds[facei];
+                }
+
+                const Type v(*valIter);
+                ++valIter;
+
+                label nValues =
+                    max
+                    (
+                        label(1),
+                        (decompOffsets[facei+1] - decompOffsets[facei])
+                    );
+
+                while (nValues--)
+                {
+                    writeFaceValue(os, v, elemId);  // 0-based
+                    ++elemId;
+                }
+            }
+        }
+
+        os  << "**" << nl
+            << "**" << nl;
+    }
+
+    wroteGeom_ = true;
+    return outputFile;
+}
+
+
+// ************************************************************************* //