Commit 7cb51f5b authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: limit vtk floatField range (fixes #1055)

- for space-savings the VTK fields are normally written as 'float'
  rather than double. When a double field contains very large values,
  these can result in a overflow when converted to float.

  Now trap these with the appropriate numeric limits.
  No warning when these values are clipped: it should be readily
  apparent from the output.

ENH: handle symmTensor component swapping directly on VTK output.

- use VTK output routines in vtkSurfaceWriter to benefit from the
  above changes
parent 0ce7e364
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -25,6 +25,7 @@ License
#include "foamVtkAppendRawFormatter.H"
#include "foamVtkOutputOptions.H"
#include <limits>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -105,8 +106,21 @@ void Foam::vtk::appendRawFormatter::write(const float val)
void Foam::vtk::appendRawFormatter::write(const double val)
{
// std::cerr<<"double as float=" << val << '\n';
float copy(val);
write(copy);
// Limit range of double to float conversion
if (val >= std::numeric_limits<float>::max())
{
write(std::numeric_limits<float>::max());
}
else if (val <= std::numeric_limits<float>::lowest())
{
write(std::numeric_limits<float>::lowest());
}
else
{
float copy(val);
write(copy);
}
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -25,6 +25,7 @@ License
#include "foamVtkAsciiFormatter.H"
#include "foamVtkOutputOptions.H"
#include <limits>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -139,8 +140,20 @@ void Foam::vtk::asciiFormatter::write(const float val)
void Foam::vtk::asciiFormatter::write(const double val)
{
next();
os()<< float(val);
// Limit range of double to float conversion
if (val >= std::numeric_limits<float>::max())
{
write(std::numeric_limits<float>::max());
}
else if (val <= std::numeric_limits<float>::lowest())
{
write(std::numeric_limits<float>::lowest());
}
else
{
float copy(val);
write(copy);
}
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "foamVtkBase64Layer.H"
#include <limits>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -96,8 +97,21 @@ void Foam::vtk::foamVtkBase64Layer::write(const float val)
void Foam::vtk::foamVtkBase64Layer::write(const double val)
{
// std::cerr<<"double as float=" << val << '\n';
float copy(val);
write(copy);
// Limit range of double to float conversion
if (val >= std::numeric_limits<float>::max())
{
write(std::numeric_limits<float>::max());
}
else if (val <= std::numeric_limits<float>::lowest())
{
write(std::numeric_limits<float>::lowest());
}
else
{
float copy(val);
write(copy);
}
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -26,6 +26,7 @@ License
#include "foamVtkLegacyRawFormatter.H"
#include "foamVtkOutputOptions.H"
#include "endian.H"
#include <limits>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -98,10 +99,7 @@ void Foam::vtk::legacyRawFormatter::write
}
void Foam::vtk::legacyRawFormatter::write
(
const label val
)
void Foam::vtk::legacyRawFormatter::write(const label val)
{
// std::cerr<<"label is:" << sizeof(val) << '\n';
......@@ -121,10 +119,7 @@ void Foam::vtk::legacyRawFormatter::write
}
void Foam::vtk::legacyRawFormatter::write
(
const float val
)
void Foam::vtk::legacyRawFormatter::write(const float val)
{
// std::cerr<<"float is:" << sizeof(val) << '\n';
......@@ -142,15 +137,25 @@ void Foam::vtk::legacyRawFormatter::write
}
void Foam::vtk::legacyRawFormatter::write
(
const double val
)
void Foam::vtk::legacyRawFormatter::write(const double val)
{
// Legacy cannot support Float64 anyhow.
// std::cerr<<"write double as float:" << val << '\n';
float copy(val);
write(copy);
// Limit range of double to float conversion
if (val >= std::numeric_limits<float>::max())
{
write(std::numeric_limits<float>::max());
}
else if (val <= std::numeric_limits<float>::lowest())
{
write(std::numeric_limits<float>::lowest());
}
else
{
float copy(val);
write(copy);
}
}
......
......@@ -48,6 +48,7 @@ SourceFiles
#include "foamVtkCore.H"
#include "foamVtkFormatter.H"
#include "floatScalar.H"
#include "symmTensor.H"
#include "IOstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -231,6 +232,22 @@ namespace legacy
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Template specialization for symmTensor ordering
template<>
inline void write(vtk::formatter& fmt, const symmTensor& val)
{
// symmTensor ( XX, XY, XZ, YY, YZ, ZZ )
// VTK order ( XX, YY, ZZ, XY, YZ, XZ ) -> (0, 3, 5, 1, 4, 2)
fmt.write(component(val, 0)); // XX
fmt.write(component(val, 3)); // YY
fmt.write(component(val, 5)); // ZZ
fmt.write(component(val, 1)); // XY
fmt.write(component(val, 4)); // YZ
fmt.write(component(val, 2)); // XZ
}
} // End namespace vtk
} // End namespace Foam
......
......@@ -25,9 +25,8 @@ License
#include "vtkSurfaceWriter.H"
#include "foamVtkOutputOptions.H"
#include "OFstream.H"
#include "OSspecific.H"
#include <fstream>
#include "makeSurfaceWriterMethods.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -49,7 +48,7 @@ defineSurfaceWriterWriteFields(Foam::vtkSurfaceWriter);
void Foam::vtkSurfaceWriter::writeGeometry
(
Ostream& os,
vtk::formatter& format,
const meshedSurf& surf,
std::string title
)
......@@ -62,153 +61,40 @@ void Foam::vtkSurfaceWriter::writeGeometry
title = "sampleSurface";
}
// header
os
<< "# vtk DataFile Version 2.0" << nl
<< title.c_str() << nl
<< "ASCII" << nl
<< "DATASET POLYDATA" << nl;
vtk::legacy::fileHeader
(
format,
title,
vtk::fileTag::POLY_DATA
);
// Write vertex coords
os << "POINTS " << points.size() << " double" << nl;
for (const point& pt : points)
{
os << float(pt.x()) << ' '
<< float(pt.y()) << ' '
<< float(pt.z()) << nl;
}
os << nl;
vtk::legacy::beginPoints(format.os(), points.size());
vtk::writeList(format, points);
format.flush();
// Write faces
label nNodes = 0;
// connectivity count without additional storage (done internally)
label nConnectivity = 0;
for (const face& f : faces)
{
nNodes += f.size();
nConnectivity += f.size();
}
os << "POLYGONS " << faces.size() << ' '
<< faces.size() + nNodes << nl;
for (const face& f : faces)
{
os << f.size();
for (const label verti : f)
{
os << ' ' << verti;
}
os << nl;
}
}
namespace Foam
{
template<>
void Foam::vtkSurfaceWriter::writeData
(
Ostream& os,
const Field<scalar>& values
)
{
os << "1 " << values.size() << " double" << nl;
forAll(values, elemI)
{
if (elemI)
{
if (elemI % 10)
{
os << ' ';
}
else
{
os << nl;
}
}
os << float(values[elemI]);
}
os << nl;
}
template<>
void Foam::vtkSurfaceWriter::writeData
(
Ostream& os,
const Field<vector>& values
)
{
os << "3 " << values.size() << " double" << nl;
for (const vector& v : values)
{
os << float(v[0]) << ' '
<< float(v[1]) << ' '
<< float(v[2]) << nl;
}
}
template<>
void Foam::vtkSurfaceWriter::writeData
vtk::legacy::beginPolys
(
Ostream& os,
const Field<sphericalTensor>& values
)
{
os << "1 " << values.size() << " double" << nl;
for (const sphericalTensor& v : values)
{
os << float(v[0]) << nl;
}
}
format.os(),
faces.size(),
nConnectivity
);
template<>
void Foam::vtkSurfaceWriter::writeData
(
Ostream& os,
const Field<symmTensor>& values
)
for (const face& f : faces)
{
os << "6 " << values.size() << " double" << nl;
// symmTensor ( XX, XY, XZ, YY, YZ, ZZ )
// VTK order ( XX, YY, ZZ, XY, YZ, XZ ) -> (0, 3, 5, 1, 4, 2)
for (const symmTensor& v : values)
{
os << float(v[0]) << ' ' << float(v[3]) << ' ' << float(v[5])
<< ' '
<< float(v[1]) << ' ' << float(v[4]) << ' ' << float(v[2])
<< nl;
}
format.write(f.size()); // The size prefix
vtk::writeList(format, f);
}
template<>
void Foam::vtkSurfaceWriter::writeData
(
Ostream& os,
const Field<tensor>& values
)
{
os << "9 " << values.size() << " double" << nl;
for (const tensor& v : values)
{
os << float(v[0]) << ' ' << float(v[1]) << ' ' << float(v[2])
<< ' '
<< float(v[3]) << ' ' << float(v[4]) << ' ' << float(v[5])
<< ' '
<< float(v[6]) << ' ' << float(v[7]) << ' ' << float(v[8])
<< nl;
}
}
format.flush();
}
......@@ -235,11 +121,14 @@ Foam::vtkSurfaceWriter::vtkSurfaceWriter(const dictionary& options)
vtk::outputOptions opts(static_cast<vtk::formatType>(fmtType_));
opts.ascii
(
options.found("format")
&& (IOstream::formatEnum(options.get<word>("format")) == IOstream::ASCII)
);
const word formatName = options.lookupOrDefault<word>("format", "");
if (formatName.size())
{
opts.ascii
(
IOstream::formatEnum(formatName) == IOstream::ASCII
);
}
if (options.lookupOrDefault("legacy", false))
{
......@@ -272,22 +161,29 @@ Foam::fileName Foam::vtkSurfaceWriter::write
{
// geometry: rootdir/time/surfaceName.{vtk|vtp}
fileName outputFile(outputDir/surfaceName + ".vtk");
if (!isDir(outputDir))
{
mkDir(outputDir);
}
OFstream os(outputDir/surfaceName + ".vtk");
os.precision(precision_);
if (verbose)
{
Info<< "Writing geometry to " << os.name() << endl;
Info<< "Writing geometry to " << outputFile << endl;
}
writeGeometry(os, surf, surfaceName);
// As vtk::outputOptions
vtk::outputOptions opts(static_cast<vtk::formatType>(fmtType_));
opts.legacy(true);
opts.precision(precision_);
std::ofstream os(outputFile);
auto format = opts.newFormatter(os);
writeGeometry(*format, surf, surfaceName);
return os.name();
return outputFile;
}
......
......@@ -85,6 +85,12 @@ SourceFiles
namespace Foam
{
namespace vtk
{
// Forward declarations
class formatter;
}
/*---------------------------------------------------------------------------*\
Class vtkSurfaceWriter Declaration
\*---------------------------------------------------------------------------*/
......@@ -107,13 +113,17 @@ class vtkSurfaceWriter
static void writeGeometry
(
Ostream& os,
vtk::formatter& format,
const meshedSurf& surf,
std::string title = ""
);
template<class Type>
static void writeData(Ostream&, const Field<Type>&);
static void writeField
(
vtk::formatter& format,
const Field<Type>& values
);
//- Templated write operation
......
......@@ -23,28 +23,69 @@ License
\*---------------------------------------------------------------------------*/
#include "OFstream.H"
#include "foamVtkOutputOptions.H"
#include "OSspecific.H"
#include <fstream>
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type>
static inline void beginFloatField
(
vtk::formatter& format,
const word& fieldName,
const Field<Type>& values
)
{
// Use 'double' instead of 'float' for ASCII output (issue #891)
if (format.opts().ascii())
{
format.os()
<< fieldName << ' '
<< int(pTraits<Type>::nComponents) << ' '
<< values.size() << " double" << nl;
}
else
{
format.os()
<< fieldName << ' '
<< int(pTraits<Type>::nComponents) << ' '
<< values.size() << " float" << nl;
}
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// // Unspecialized (unknown) data type - map as zeros
// void Foam::vtkSurfaceWriter::writeZeros
// (
// vtk::formatter& format,
// const label count
// )
// {
// const float val(0);
//
// for (label i=0; i < count; ++i)
// {
// format.write(val);
// }
// format.flush();
// }
template<class Type>
void Foam::vtkSurfaceWriter::writeData
void Foam::vtkSurfaceWriter::writeField
(
Ostream& os,
vtk::formatter& format,
const Field<Type>& values
)
{
// Unspecialized (unknown) data type - map as zeros
const label len = values.size();
os << "1 " << len << " double" << nl;
for (label i=0; i < len; ++i)
{
os << float(0) << nl;
}
vtk::writeList(format, values);
format.flush();
}
......@@ -62,39 +103,48 @@ Foam::fileName Foam::vtkSurfaceWriter::writeTemplate
{
// field: rootdir/time/<field>_surfaceName.{vtk|vtp}
fileName outputFile(outputDir/fieldName + '_' + surfaceName + ".vtk");
if (!isDir(outputDir))
{
mkDir(outputDir);
}
OFstream os(outputDir/fieldName + '_' + surfaceName + ".vtk");
os.precision(precision_);
if (verbose)
{
Info<< "Writing field " << fieldName << " to " << os.name() << endl;
Info<< "Writing field " << fieldName << " to "
<< outputFile << endl;
}