zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-noise' into 'develop'

Browse Source 
Feature noise

New functionality includes:
- run-time selectable noise models: point|surface
- run-time selectable window models: Hanning (+ options symmetric, extended), uniform
- calculates PSD (Pa^2/Hz) and dB/HZ; SPL (Pa^2) and dB
- calculates 1/3 octave data, with centre frequency 1kHz

surfaceNoise only:
- reads ascii/binary ensight surface data (requires collateTimes option)
- generates graphs for surface average quantities
- operates in parallel

See merge request !50
This commit is contained in:
Andrew Heather
2016-06-29 21:44:36 +01:00
parent 038356cdda 5062c50ebc
commit 4b266671a3
41 changed files with 4709 additions and 902 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
applications/solvers/DNS/dnsFoam/Allwmake Executable file
View File

@ -0,0 +1,14 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \
[ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ]
then
wmake
else
echo
echo "Skipping dnsFoam solver (no FFTW)"
echo
fi
#------------------------------------------------------------------------------

14
applications/utilities/postProcessing/noise/Allwmake Executable file
View File

@ -0,0 +1,14 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \
[ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ]
then
wmake
else
echo
echo "Skipping noise utility (no FFTW)"
echo
fi
#------------------------------------------------------------------------------

3
applications/utilities/postProcessing/noise/Make/options
View File

@ -1,6 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/randomProcesses/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
EXE_LIBS = \
-lrandomProcesses \

184
applications/utilities/postProcessing/noise/noise.C
View File

@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -28,157 +28,103 @@ Group
grpPostProcessingUtilities
Description
Utility to perform noise analysis of pressure data using the noiseFFT
library.
Utility to perform noise analysis of pressure data.
Control settings are read from the $FOAM_CASE/system/noiseDict dictionary,
or user-specified dictionary using the -dict option. Pressure data is
read using a CSV reader:
The utility provides a light wrapper around the run-time selectable
noise model. Current options include:
- point, and
- surface noise.
\heading Usage
\heading Example usage
\verbatim
pRef 101325;
N 65536;
nw 100;
f1 25;
fU 10000;
graphFormat raw;
noiseModel surfaceNoise; // pointNoise
pressureData
surfaceNoiseCoeffs
{
fileName "pressureData"
nHeaderLine 1; // number of header lines
refColumn 0; // reference column index
componentColumns (1); // component column indices
separator " "; // optional (defaults to ",")
mergeSeparators no; // merge multiple separators
outOfBounds clamp; // optional out-of-bounds handling
interpolationScheme linear; // optional interpolation scheme
windowModel Hanning;
HanningCoeffs
{
// Window overlap percentage
overlapPercent 50;
symmetric yes;
extended yes;
// Optional number of windows, default = all available
nWindow 5;
}
// Input file
inputFile "postProcessing/faceSource1/surface/patch/patch.case";
// Surface reader
reader ensight;
// Surface writer
writer ensight;
// Collate times for ensight output - ensures geometry is only written once
writeOptions
{
ensight
{
collateTimes 1;
}
}
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096;
// Write interval for FFT data, default = 1
fftWriteInterval 100;
}
\endverbatim
where
\table
Property | Description | Required | Default value
pRef | Reference pressure | no | 0
N | Number of samples in sampling window | no | 65536
nw | Number of sampling windows | no | 100
fl | Lower frequency band | no | 25
fU | Upper frequency band | no | 10000
graphFormat | Output graph format | no | raw
\endtable
Current graph outputs include:
- FFT of the pressure data
- narrow-band PFL (pressure-fluctuation level) spectrum
- one-third-octave-band PFL spectrum
- one-third-octave-band pressure spectrum
SeeAlso
CSV.H
noiseFFT.H
noiseModel.H
windowModel.H
\*---------------------------------------------------------------------------*/
#include "noiseFFT.H"
#include "argList.H"
#include "Time.H"
#include "functionObjectFile.H"
#include "CSV.H"
#include "noiseModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::scalar checkUniformTimeStep(const scalarField& t)
{
// check that a uniform time step has been applied
scalar deltaT = -1.0;
if (t.size() > 1)
{
for (label i = 1; i < t.size(); i++)
{
scalar dT = t[i] - t[i-1];
if (deltaT < 0)
{
deltaT = dT;
}
if (mag(deltaT - dT) > SMALL)
{
FatalErrorInFunction
<< "Unable to process data with a variable time step"
<< exit(FatalError);
}
}
}
else
{
FatalErrorInFunction
<< "Unable to create FFT with a single value"
<< exit(FatalError);
}
return deltaT;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::noParallel();
#include "addDictOption.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createFields.H"
Info<< "Reading data file" << endl;
Function1Types::CSV<scalar> pData("pressure", dict, "Data");
// time history data
const scalarField t(pData.x());
// pressure data
const scalarField p(pData.y());
if (t.size() < N)
fileName dictName(runTime.system()/"noiseDict");
if (args.optionFound("dict"))
{
FatalErrorInFunction
<< "Block size N = " << N
<< " is larger than number of data = " << t.size()
<< exit(FatalError);
dictName = args["dict"];
}
Info<< " read " << t.size() << " values" << nl << endl;
IOdictionary dict
(
IOobject
(
dictName,
runTime,
IOobject::MUST_READ
)
);
Info<< "Creating noise FFT" << endl;
noiseFFT nfft(checkUniformTimeStep(t), p);
nfft -= pRef;
fileName baseFileName(pData.fName().lessExt());
graph Pf(nfft.RMSmeanPf(N, min(nfft.size()/N, nw)));
Info<< " Creating graph for " << Pf.title() << endl;
Pf.write(baseFileName + graph::wordify(Pf.title()), graphFormat);
graph Lf(nfft.Lf(Pf));
Info<< " Creating graph for " << Lf.title() << endl;
Lf.write(baseFileName + graph::wordify(Lf.title()), graphFormat);
graph Ldelta(nfft.Ldelta(Lf, f1, fU));
Info<< " Creating graph for " << Ldelta.title() << endl;
Ldelta.write(baseFileName + graph::wordify(Ldelta.title()), graphFormat);
graph Pdelta(nfft.Pdelta(Pf, f1, fU));
Info<< " Creating graph for " << Pdelta.title() << endl;
Pdelta.write(baseFileName + graph::wordify(Pdelta.title()), graphFormat);
autoPtr<noiseModel> model(noiseModel::New(dict));
model->calculate();
Info<< nl << "End\n" << endl;

136
applications/utilities/postProcessing/noise/noiseDict Normal file
View File

@ -0,0 +1,136 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object noiseDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
noiseModel surfaceNoise;
surfaceNoiseCoeffs
{
windowModel Hanning;
HanningCoeffs
{
// Window overlap percentage
overlapPercent 50;
symmetric yes;
extended yes;
// Optional number of windows, default = all available
// nWindow 1;
}
/*
windowModel uniform;
uniformCoeffs
{
// Window overlap percentage
overlapPercent 50;
value 1;
// Optional number of windows, default = all available
// nWindow 1;
}
*/
// Input file
inputFile "postProcessing/faceSource1/surface/patch_motorBike_rider-helmet%65/patch_motorBike_rider-helmet%65.case";
// Surface reader
reader ensight;
// Surface writer
writer ensight;
// Collate times for ensight output - ensures geometry is only written once
writeOptions
{
ensight
{
collateTimes 1;
}
}
// Reference density (to convert from kinematic to static pressure)
rhoRef 1.205;
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096; // 8192; // 4096;
// Lower frequency limit, default = 25Hz
//fl 25;
// Upper frequency limit, default = 10kHz
fu 15000;
// Start time, default = 0s
//startTime 0;
// Write interval for FFT data, default = 1
// fftWriteInterval 100;
}
pointNoiseCoeffs
{
csvFileData
{
fileName "pressureData";
nHeaderLine 1;
refColumn 0;
componentColumns (1);
separator " ";
mergeSeparators yes;
}
HanningCoeffs
{
// Window overlap percentage
overlapPercent 50;
symmetric yes;
extended yes;
// Optional number of windows, default = all available
//nWindow 5;
}
// Graph format, default = raw
graphFormat raw;
// Reference density (to convert from kinematic to static pressure)
rhoRef 1.2;
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096;
// Lower frequency limit, default = 25Hz
//fl 25;
// Upper frequency limit, default = 10kHz
//fu 10000;
// Start time, default = 0s
//startTime 0;
// Write interval for FFT data, default = 1
fftWriteInterval 100;
}
// ************************************************************************* //

14
applications/utilities/preProcessing/boxTurb/Allwmake Executable file
View File

@ -0,0 +1,14 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \
[ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ]
then
wmake
else
echo
echo "Skipping boxTurb utility (no FFTW)"
echo
fi
#------------------------------------------------------------------------------

2
etc/cshrc
View File

@ -31,7 +31,7 @@
#------------------------------------------------------------------------------
setenv WM_PROJECT OpenFOAM
setenv WM_PROJECT_VERSION plus
setenv WM_PROJECT_VERSION stage
################################################################################
# USER EDITABLE PART: Changes made here may be lost with the next upgrade

1
src/conversion/Make/files
View File

@ -1,5 +1,6 @@
ensight/file/ensightFile.C
ensight/file/ensightGeoFile.C
ensight/readFile/ensightReadFile.C
ensight/part/ensightPart.C
ensight/part/ensightPartIO.C
ensight/part/ensightPartCells.C

157
src/conversion/ensight/readFile/ensightReadFile.C Normal file
View File

@ -0,0 +1,157 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "ensightReadFile.H"
#include <sstream>
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ensightReadFile::ensightReadFile
(
const fileName& pathname,
IOstream::streamFormat format
)
:
IFstream(pathname, format)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::ensightReadFile::~ensightReadFile()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::Istream& Foam::ensightReadFile::read
(
char* buf,
std::streamsize count
)
{
stdStream().read(buf, count);
return *this;
}
Foam::Istream& Foam::ensightReadFile::read(string& value)
{
if (format() == IOstream::BINARY)
{
char buf[80];
read(reinterpret_cast<char*>(buf), sizeof(buf));
string strBuf(value);
const size_t iEnd = strBuf.find('\0', 0);
if (iEnd == string::npos)
{
value = buf;
}
else
{
value = strBuf.substr(0, iEnd - 1);
}
}
else
{
value = "";
while (value.empty() && !eof())
{
getLine(value);
}
}
return *this;
}
Foam::Istream& Foam::ensightReadFile::read(label& value)
{
int ivalue;
if (format() == IOstream::BINARY)
{
read
(
reinterpret_cast<char*>(&ivalue),
sizeof(ivalue)
);
}
else
{
stdStream() >> ivalue;
}
value = ivalue;
return *this;
}
Foam::Istream& Foam::ensightReadFile::read(scalar& value)
{
float fvalue;
if (format() == IOstream::BINARY)
{
read
(
reinterpret_cast<char*>(&fvalue),
sizeof(fvalue)
);
value = fvalue;
}
else
{
stdStream() >> value;
}
return *this;
}
Foam::Istream& Foam::ensightReadFile::readKeyword(string& key)
{
read(key);
return *this;
}
Foam::Istream& Foam::ensightReadFile::readBinaryHeader()
{
if (format() == IOstream::BINARY)
{
string buffer;
read(buffer);
}
return *this;
}
// ************************************************************************* //

110
src/conversion/ensight/readFile/ensightReadFile.H Normal file
View File

@ -0,0 +1,110 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::ensightReadFile
Description
Ensight output with specialized read() for strings, integers and floats.
Correctly handles binary read as well.
\*---------------------------------------------------------------------------*/
#ifndef ensightReadFile_H
#define ensightReadFile_H
#include "IFstream.H"
#include "IOstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ensightReadFile Declaration
\*---------------------------------------------------------------------------*/
class ensightReadFile
:
public IFstream
{
// Private Member Functions
//- Disallow default bitwise assignment
void operator=(const ensightReadFile&);
//- Disallow default copy constructor
ensightReadFile(const ensightReadFile&);
public:
// Constructors
//- Construct from pathname
ensightReadFile
(
const fileName& pathname,
IOstream::streamFormat format=IOstream::BINARY
);
//- Destructor
~ensightReadFile();
// Output
//- Inherit read from Istream
using Istream::read;
//- Binary read
virtual Istream& read(char* buf, std::streamsize count);
//- Read string as "%80s" or as binary
Istream& read(string& value);
//- Read integer as "%10d" or as binary
Istream& read(label& value);
//- Read float as "%12.5e" or as binary
Istream& read(scalar& value);
//- Read element keyword
virtual Istream& readKeyword(string& key);
//- Read "C Binary" for binary files (eg, geometry/measured)
Istream& readBinaryHeader();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

18
src/randomProcesses/Allwmake Executable file
View File

@ -0,0 +1,18 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
# Parse arguments for library compilation
targetType=libso
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
if [ -f "$FFTW_ARCH_PATH/include/fftw3.h" ] || \
[ "${FFTW_ARCH_PATH##*-}" = system -a -f "/usr/include/fftw3.h" ]
then
wmake $targetType
else
echo
echo "Skipping randomProcesses library (no FFTW)"
echo
fi
#------------------------------------------------------------------------------

32
src/randomProcesses/Make/files
View File

@ -1,25 +1,31 @@
Kmesh = Kmesh
fft = fft
processes = processes
UOprocess = $(processes)/UOprocess
turbulence = turbulence
noise = noise
Kmesh = Kmesh
$(Kmesh)/Kmesh.C
fft = fft
$(fft)/fft.C
$(fft)/fftRenumber.C
$(fft)/calcEk.C
$(fft)/kShellIntegration.C
processes = processes
UOprocess = $(processes)/UOprocess
$(UOprocess)/UOprocess.C
turbulence = turbulence
$(turbulence)/turbGen.C
$(noise)/noiseFFT.C
noise = noise
$(noise)/noiseFFT/noiseFFT.C
$(noise)/noiseModels/noiseModel/noiseModel.C
$(noise)/noiseModels/noiseModel/noiseModelNew.C
$(noise)/noiseModels/pointNoise/pointNoise.C
$(noise)/noiseModels/surfaceNoise/surfaceNoise.C
windowModels = windowModels
$(windowModels)/windowModel/windowModel.C
$(windowModels)/windowModel/windowModelNew.C
$(windowModels)/Hanning/Hanning.C
$(windowModels)/uniform/uniform.C
LIB = $(FOAM_LIBBIN)/librandomProcesses

10
src/randomProcesses/Make/options
View File

@ -1,5 +1,11 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(FFTW_ARCH_PATH)/include \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
LIB_LIBS = \
-lfiniteVolume
-L$(FFTW_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -lfftw3 \
-lfiniteVolume \
-lsampling \
-lsurfMesh

174
src/randomProcesses/fft/fft.C
View File

@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -24,7 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "fft.H"
#include "fftRenumber.H"
#include <fftw3.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -33,16 +34,11 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr
#define TWOPI 6.28318530717959
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void fft::transform
(
complexField& field,
const labelList& nn,
transformDirection isign
transformDirection dir
)
{
forAll(nn, idim)
@ -59,128 +55,58 @@ void fft::transform
}
}
const label ndim = nn.size();
label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
label ibit, k1, k2, n, nprev, nrem, idim;
scalar tempi, tempr;
scalar theta, wi, wpi, wpr, wr, wtemp;
scalar* data = reinterpret_cast<scalar*>(field.begin()) - 1;
// if inverse transform : renumber before transform
if (isign == REVERSE_TRANSFORM)
{
fftRenumber(field, nn);
}
label ntot = 1;
forAll(nn, idim)
{
ntot *= nn[idim];
}
nprev = 1;
for (idim=ndim; idim>=1; idim--)
{
n = nn[idim-1];
nrem = ntot/(n*nprev);
ip1 = nprev << 1;
ip2 = ip1*n;
ip3 = ip2*nrem;
i2rev = 1;
for (i2=1; i2<=ip2; i2+=ip1)
{
if (i2 < i2rev)
{
for (i1=i2; i1<=i2 + ip1 - 2; i1+=2)
{
for (i3=i1; i3<=ip3; i3+=ip2)
{
i3rev = i2rev + i3 - i2;
SWAP(data[i3], data[i3rev]);
SWAP(data[i3 + 1], data[i3rev + 1]);
}
}
}
ibit = ip2 >> 1;
while (ibit >= ip1 && i2rev > ibit)
{
i2rev -= ibit;
ibit >>= 1;
}
i2rev += ibit;
}
ifp1 = ip1;
while (ifp1 < ip2)
{
ifp2 = ifp1 << 1;
theta = isign*TWOPI/(ifp2/ip1);
wtemp = sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for (i3 = 1; i3 <= ifp1; i3 += ip1)
{
for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
{
for (i2 = i1; i2 <= ip3; i2 += ifp2)
{
k1 = i2;
k2 = k1 + ifp1;
tempr = scalar(wr*data[k2]) - scalar(wi*data[k2 + 1]);
tempi = scalar(wr*data[k2 + 1]) + scalar(wi*data[k2]);
data[k2] = data[k1] - tempr;
data[k2 + 1] = data[k1 + 1] - tempi;
data[k1] += tempr;
data[k1 + 1] += tempi;
}
}
wr = (wtemp = wr)*wpr - wi*wpi + wr;
wi = wi*wpr + wtemp*wpi + wi;
}
ifp1 = ifp2;
}
nprev *= n;
}
// if forward transform : renumber after transform
if (isign == FORWARD_TRANSFORM)
{
fftRenumber(field, nn);
}
// scale result (symmetric scale both forward and inverse transform)
scalar recRootN = 1.0/sqrt(scalar(ntot));
// Copy field into fftw containers
label N = field.size();
fftw_complex in[N], out[N];
forAll(field, i)
{
field[i] *= recRootN;
in[i][0] = field[i].Re();
in[i][1] = field[i].Im();
}
// Create the plan
// FFTW_FORWARD = -1
// FFTW_BACKWARD = 1
// 1-D plan
// fftw_plan plan =
// fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
// Generic 1..3-D plan
label rank = nn.size();
fftw_plan plan =
fftw_plan_dft(rank, nn.begin(), in, out, dir, FFTW_ESTIMATE);
// Compute the FFT
fftw_execute(plan);
forAll(field, i)
{
field[i].Re() = out[i][0];
field[i].Im() = out[i][1];
}
fftw_destroy_plan(plan);
/*
fftw_real in[N], out[N];
// Create a plan for real data input
// - generates 1-sided FFT
// - direction given by FFTW_R2HC or FFTW_HC2R.
// - in[N] is the array of real input val ues
// - out[N] is the array of N/2 real valuea followed by N/2 complex values
// - 0 component = DC component
fftw_plan plan = fftw_plan_r2r_2d(N, in, out, FFTW_R2HC, FFTW_ESTIMATE)
// Compute the FFT
fftw_execute(plan);
fftw_destroy_plan(plan);
*/
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#undef SWAP
#undef TWOPI
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
tmp<complexField> fft::forwardTransform

11
src/randomProcesses/fft/fft.H
View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -25,11 +25,10 @@ Class
Foam::fft
Description
Fast fourier transform derived from the Numerical
Recipes in C routine.
Fast fourier transform using the fftw library.
The complex transform field is returned in the field supplied. The
direction of transform is supplied as an argument (1 = forward, -1 =
direction of transform is supplied as an argument (-1 = forward, 1 =
reverse). The dimensionality and organisation of the array of values
in space is supplied in the nn indexing array.
@ -56,8 +55,8 @@ public:
enum transformDirection
{
FORWARD_TRANSFORM = 1,
REVERSE_TRANSFORM = -1
FORWARD_TRANSFORM = -1,
REVERSE_TRANSFORM = 1
};

124
src/randomProcesses/fft/fftRenumber.C
View File

@ -1,124 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Multi-dimensional renumbering used in the Numerical Recipes
fft routine. This version is recursive, so works in n-d :
determined by the length of array nn
\*---------------------------------------------------------------------------*/
#include "fftRenumber.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::fftRenumberRecurse
(
List<complex>& data,
List<complex>& renumData,
const labelList& nn,
label nnprod,
label ii,
label l1,
label l2
)
{
if (ii == nn.size())
{
// we've worked out the renumbering scheme. Now copy
// the components across
data[l1] = complex(renumData[l2].Re(),renumData[l2].Im());
}
else
{
// do another level of folding. First work out the
// multiplicative value of the index
nnprod /= nn[ii];
label i_1(0);
for (label i=0; i<nn[ii]; i++)
{
// now evaluate the indices (both from array 1 and to
// array 2). These get multiplied by nnprod to (cumulatively)
// find the real position in the list corresponding to
// this set of indices.
if (i<nn[ii]/2)
{
i_1 = i + nn[ii]/2;
}
else
{
i_1 = i - nn[ii]/2;
}
// go to the next level of recursion.
fftRenumberRecurse
(
data,
renumData,
nn,
nnprod,
ii+1,
l1+i*nnprod,
l2+i_1*nnprod
);
}
}
}
void Foam::fftRenumber
(
List<complex>& data,
const labelList& nn
)
{
List<complex> renumData(data);
label nnprod(1);
forAll(nn, i)
{
nnprod *= nn[i];
}
label ii(0), l1(0), l2(0);
fftRenumberRecurse
(
data,
renumData,
nn,
nnprod,
ii,
l1,
l2
);
}
// ************************************************************************* //

454
src/randomProcesses/noise/noiseFFT.C
View File

@ -1,454 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "noiseFFT.H"
#include "IFstream.H"
#include "DynamicList.H"
#include "fft.H"
#include "SubField.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
Foam::scalar Foam::noiseFFT::p0 = 2e-5;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::noiseFFT::noiseFFT
(
const scalar deltat,
const scalarField& pressure
)
:
scalarField(pressure),
deltat_(deltat)
{}
Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
:
scalarField(),
deltat_(0.0)
{
// Construct pressure data file
IFstream pFile(pFileName);
// Check pFile stream is OK
if (!pFile.good())
{
FatalErrorInFunction
<< "Cannot read file " << pFileName
<< exit(FatalError);
}
if (skip)
{
scalar dummyt, dummyp;
for (label i=0; i<skip; i++)
{
pFile >> dummyt;
if (!pFile.good() || pFile.eof())
{
FatalErrorInFunction
<< "Number of points in file " << pFileName
<< " is less than the number to be skipped = " << skip
<< exit(FatalError);
}
pFile >> dummyp;
}
}
scalar t = 0, T = 0;
DynamicList<scalar> pData(100000);
label i = 0;
while (!(pFile >> t).eof())
{
T = t;
pFile >> pData(i++);
}
deltat_ = T/pData.size();
this->transfer(pData);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::graph Foam::noiseFFT::pt() const
{
scalarField t(size());
forAll(t, i)
{
t[i] = i*deltat_;
}
return graph
(
"p(t)",
"t [s]",
"p(t) [Pa]",
t,
*this
);
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
(
const label N,
const label ni
) const
{
label windowOffset = N;
if ((N + ni*windowOffset) > size())
{
FatalErrorInFunction
<< "Requested window is outside set of data" << endl
<< "number of data = " << size() << endl
<< "size of window = " << N << endl
<< "window = " << ni
<< exit(FatalError);
}
tmp<scalarField> tpw(new scalarField(N));
scalarField& pw = tpw.ref();
label offset = ni*windowOffset;
forAll(pw, i)
{
pw[i] = operator[](i + offset);
}
return tpw;
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
{
scalarField t(N);
forAll(t, i)
{
t[i] = i*deltat_;
}
scalar T = N*deltat_;
return 2*(0.5 - 0.5*cos(constant::mathematical::twoPi*t/T));
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
(
const tmp<scalarField>& tpn
) const
{
tmp<scalarField> tPn2
(
mag
(
fft::reverseTransform
(
ReComplexField(tpn),
labelList(1, tpn().size())
)
)
);
tpn.clear();
tmp<scalarField> tPn
(
new scalarField
(
scalarField::subField(tPn2(), tPn2().size()/2)
)
);
scalarField& Pn = tPn.ref();
Pn *= 2.0/sqrt(scalar(tPn2().size()));
Pn[0] /= 2.0;
return tPn;
}
Foam::graph Foam::noiseFFT::meanPf
(
const label N,
const label nw
) const
{
if (N > size())
{
FatalErrorInFunction
<< "Requested window is outside set of data" << nl
<< "number of data = " << size() << nl
<< "size of window = " << N << nl
<< "window = " << nw
<< exit(FatalError);
}
scalarField MeanPf(N/2, 0.0);
scalarField Hwf(Hanning(N));
for (label wi=0; wi<nw; ++wi)
{
MeanPf += Pf(Hwf*window(N, wi));
}
MeanPf /= nw;
scalarField f(MeanPf.size());
scalar deltaf = 1.0/(N*deltat_);
forAll(f, i)
{
f[i] = i*deltaf;
}
return graph
(
"P(f)",
"f [Hz]",
"P(f) [Pa]",
f,
MeanPf
);
}
Foam::graph Foam::noiseFFT::RMSmeanPf
(
const label N,
const label nw
) const
{
if (N > size())
{
FatalErrorInFunction
<< "Requested window is outside set of data" << endl
<< "number of data = " << size() << endl
<< "size of window = " << N << endl
<< "window = " << nw
<< exit(FatalError);
}
scalarField RMSMeanPf(N/2, 0.0);
scalarField Hwf(Hanning(N));
for (label wi=0; wi<nw; ++wi)
{
RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
}
RMSMeanPf = sqrt(RMSMeanPf/nw);
scalarField f(RMSMeanPf.size());
scalar deltaf = 1.0/(N*deltat_);
forAll(f, i)
{
f[i] = i*deltaf;
}
return graph
(
"P(f)",
"f [Hz]",
"P(f) [Pa]",
f,
RMSMeanPf
);
}
Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
{
return graph
(
"L(f)",
"f [Hz]",
"L(f) [dB]",
gPf.x(),
20*log10(gPf.y()/p0)
);
}
Foam::graph Foam::noiseFFT::Ldelta
(
const graph& gLf,
const scalar f1,
const scalar fU
) const
{
const scalarField& f = gLf.x();
const scalarField& Lf = gLf.y();
scalarField ldelta(Lf.size(), 0.0);
scalarField fm(ldelta.size());
scalar fratio = cbrt(2.0);
scalar deltaf = 1.0/(2*Lf.size()*deltat_);
scalar fl = f1/sqrt(fratio);
scalar fu = fratio*fl;
label istart = label(fl/deltaf);
label j = 0;
for (label i = istart; i<Lf.size(); i++)
{
scalar fmi = sqrt(fu*fl);
if (fmi > fU + 1) break;
if (f[i] >= fu)
{
fm[j] = fmi;
ldelta[j] = 10*log10(ldelta[j]);
j++;
fl = fu;
fu *= fratio;
}
ldelta[j] += pow(10, Lf[i]/10.0);
}
fm.setSize(j);
ldelta.setSize(j);
return graph
(
"Ldelta",
"fm [Hz]",
"Ldelta [dB]",
fm,
ldelta
);
}
Foam::graph Foam::noiseFFT::Pdelta
(
const graph& gPf,
const scalar f1,
const scalar fU
) const
{
const scalarField& f = gPf.x();
const scalarField& Pf = gPf.y();
scalarField pdelta(Pf.size(), 0.0);
scalarField fm(pdelta.size());
scalar fratio = cbrt(2.0);
scalar deltaf = 1.0/(2*Pf.size()*deltat_);
scalar fl = f1/sqrt(fratio);
scalar fu = fratio*fl;
label istart = label(fl/deltaf + 1.0 - SMALL);
label j = 0;
for (label i = istart; i<Pf.size(); i++)
{
scalar fmi = sqrt(fu*fl);
if (fmi > fU + 1) break;
if (f[i] >= fu)
{
fm[j] = fmi;
pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
j++;
fl = fu;
fu *= fratio;
}
pdelta[j] += sqr(Pf[i]);
}
fm.setSize(j);
pdelta.setSize(j);
return graph
(
"Pdelta",
"fm [Hz]",
"Pdelta [dB]",
fm,
pdelta
);
}
Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
{
const scalarField& Lf = gLf.y();
scalar lsum = 0.0;
forAll(Lf, i)
{
lsum += pow(10, Lf[i]/10.0);
}
lsum = 10*log10(lsum);
return lsum;
}
Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
{
return p0*pow(10.0, db/20.0);
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
(
const tmp<scalarField>& db
) const
{
return p0*pow(10.0, db/20.0);
}
// ************************************************************************* //

474
src/randomProcesses/noise/noiseFFT/noiseFFT.C Normal file
View File

@ -0,0 +1,474 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 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 "noiseFFT.H"
#include "IFstream.H"
#include "DynamicList.H"
#include "SubField.H"
#include "mathematicalConstants.H"
#include "HashSet.H"
#include "fft.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * //
Foam::scalar Foam::noiseFFT::p0 = 2e-5;
Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
(
const label N,
const scalar deltaT
)
{
tmp<scalarField> tf(new scalarField(N/2, 0));
scalarField& f = tf.ref();
scalar deltaf = 1.0/(N*deltaT);
forAll(f, i)
{
f[i] = i*deltaf;
}
return tf;
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::PSD(const scalarField& PSDf)
{
return 10*log10(PSDf/sqr(p0));
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::SPL(const scalarField& Prms2)
{
return 10*log10(Prms2/sqr(p0));
}
void Foam::noiseFFT::octaveBandInfo
(
const scalarField& f,
const scalar fLower,
const scalar fUpper,
const scalar octave,
labelList& fBandIDs,
scalarField& fCentre
)
{
// Set the indices of to the lower frequency bands for the input frequency
// range. Ensure that the centre frequency passes though 1000 Hz
// Low frequency bound given by:
// fLow = f0*(2^(0.5*bandI/octave))
// Initial (lowest centre frequency)
scalar fTest = 15.625;
const scalar fRatio = pow(2, 1.0/octave);
const scalar fRatioL2C = pow(2, 0.5/octave);
// IDs of band IDs
labelHashSet bandIDs(f.size());
// Centre frequencies
DynamicList<scalar> fc;
// Convert to lower band limit
fTest /= fRatioL2C;
forAll(f, i)
{
if (f[i] >= fTest)
{
// Advance band if appropriate
while (f[i] > fTest)
{
fTest *= fRatio;
}
fTest /= fRatio;
if (bandIDs.insert(i))
{
// Also store (next) centre frequency
fc.append(fTest*fRatioL2C);
}
fTest *= fRatio;
if (fTest > fUpper)
{
break;
}
}
}
fBandIDs = bandIDs.toc();
// Remove the last centre frequency (beyond upper frequency limit)
fc.remove();
fCentre.transfer(fc);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::noiseFFT::noiseFFT
(
const scalar deltaT,
const scalarField& pressure
)
:
scalarField(pressure),
deltaT_(deltaT)
{
scalarField& p = *this;
p -= average(p);
}
Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
:
scalarField(),
deltaT_(0.0)
{
// Construct pressure data file
IFstream pFile(pFileName);
// Check pFile stream is OK
if (!pFile.good())
{
FatalErrorInFunction
<< "Cannot read file " << pFileName
<< exit(FatalError);
}
if (skip)
{
scalar dummyt, dummyp;
for (label i = 0; i < skip; i++)
{
pFile >> dummyt;
if (!pFile.good() || pFile.eof())
{
FatalErrorInFunction
<< "Number of points in file " << pFileName
<< " is less than the number to be skipped = " << skip
<< exit(FatalError);
}
pFile >> dummyp;
}
}
scalar t = 0, T0 = 0, T1 = 0;
DynamicList<scalar> pData(100000);
label i = 0;
while (!(pFile >> t).eof())
{
if (i == 0)
{
T0 = t;
}
T1 = t;
pFile >> pData(i++);
}
// Note: assumes fixed time step
deltaT_ = (T1 - T0)/pData.size();
this->transfer(pData);
scalarField& p = *this;
p -= average(p);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::graph Foam::noiseFFT::pt() const
{
scalarField t(size());
forAll(t, i)
{
t[i] = i*deltaT_;
}
return graph
(
"p(t)",
"t [s]",
"p(t) [Pa]",
t,
*this
);
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
(
const tmp<scalarField>& tpn
) const
{
// Calculate the 2-sided fft
// Note: result not scaled
tmp<scalarField> tPn2
(
mag
(
fft::reverseTransform
(
ReComplexField(tpn),
labelList(1, tpn().size())
)
)
);
tpn.clear();
// Only storing the positive half
// Note: storing (N/2) values, DC component at position (0)
tmp<scalarField> tPn
(
new scalarField
(
scalarField::subField(tPn2(), tPn2().size()/2 + 1)
)
);
return tPn;
}
Foam::graph Foam::noiseFFT::meanPf(const windowModel& window) const
{
const label N = window.nSamples();
const label nWindow = window.nWindow();
scalarField meanPf(N/2 + 1, 0.0);
for (label windowI = 0; windowI < nWindow; ++windowI)
{
meanPf += Pf(window.apply<scalar>(*this, windowI));
}
meanPf /= scalar(nWindow);
scalar deltaf = 1.0/(N*deltaT_);
scalarField f(meanPf.size());
forAll(f, i)
{
f[i] = i*deltaf;
}
return graph
(
"P(f)",
"f [Hz]",
"P(f) [Pa]",
f,
meanPf
);
}
Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
{
const label N = window.nSamples();
const label nWindow = window.nWindow();
scalarField RMSMeanPf(N/2 + 1, 0.0);
for (label windowI = 0; windowI < nWindow; ++windowI)
{
RMSMeanPf += sqr(Pf(window.apply<scalar>(*this, windowI)));
}
RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow))/scalar(N);
scalar deltaf = 1.0/(N*deltaT_);
scalarField f(RMSMeanPf.size());
forAll(f, i)
{
f[i] = i*deltaf;
}
return graph
(
"Prms(f)",
"f [Hz]",
"Prms(f) [Pa]",
f,
RMSMeanPf
);
}
Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const
{
const label N = window.nSamples();
const label nWindow = window.nWindow();
scalarField psd(N/2 + 1, 0.0);
for (label windowI = 0; windowI < nWindow; ++windowI)
{
psd += sqr(Pf(window.apply<scalar>(*this, windowI)));
}
scalar deltaf = 1.0/(N*deltaT_);
scalar fs = 1.0/deltaT_;
psd /= scalar(nWindow)*fs*N;
// Scaling due to use of 1-sided FFT
// Note: do not scale DC component
psd *= 2;
psd.first() /= 2;
psd.last() /= 2;
scalarField f(psd.size());
if (0) // if (debug)
{
Pout<< "Average PSD: " << average(psd) << endl;
}
forAll(f, i)
{
f[i] = i*deltaf;
}
return graph
(
"PSD(f)",
"f [Hz]",
"PSD(f) [PaPa_Hz]",
f,
psd
);
}
Foam::graph Foam::noiseFFT::PSD(const graph& gPSDf) const
{
return graph
(
"PSD(f)",
"f [Hz]",
"PSD_dB(f) [dB_Hz]",
gPSDf.x(),
10*log10(gPSDf.y()/sqr(p0))
);
}
Foam::graph Foam::noiseFFT::octaves
(
const graph& g,
const labelList& freqBandIDs,
bool integrate
) const
{
if (freqBandIDs.size() < 2)
{
WarningInFunction
<< "Octave frequency bands are not defined "
<< "- skipping octaves calculation"
<< endl;
return graph
(
"octave",
"x",
"y",
scalarField(),
scalarField()
);
}
const scalarField& f = g.x();
const scalarField& data = g.y();
scalarField octData(freqBandIDs.size() - 1, 0.0);
scalarField fm(freqBandIDs.size() -1, 0.0);
for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
{
label fb0 = freqBandIDs[bandI];
label fb1 = freqBandIDs[bandI+1];
fm[bandI] = f[fb0];
if (fb0 == fb1) continue;
if (integrate)
{
for (label freqI = fb0; freqI < fb1; freqI++)
{
label f0 = f[freqI];
label f1 = f[freqI + 1];
scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
octData[bandI] += dataAve*(f1 - f0);
}
}
else
{
for (label freqI = fb0; freqI < fb1; freqI++)
{
octData[bandI] += data[freqI];
}
}
}
return graph
(
"octaves(f)",
"fm [Hz]",
"octave data",
fm,
octData
);
}
Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
{
return p0*pow(10.0, db/20.0);
}
Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
(
const tmp<scalarField>& db
) const
{
return p0*pow(10.0, db/20.0);
}
// ************************************************************************* //

105
src/randomProcesses/noise/noiseFFT.H → src/randomProcesses/noise/noiseFFT/noiseFFT.H
View File

@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -25,11 +25,26 @@ Class
Foam::noiseFFT
Description
FFT of the pressure field
Performs FFT of pressure field to generate noise data.
General functionality:
- Pf: fft of the pressure data
- meanPf: multi-window mean fft
- RMSmeanPf: multi-window RMS mean fft
- PSDf: multi-window power spectral density (PSD) in frequency domain
- PSD: power spectral density in dB/Hz
- SPL: sound pressure level in dB
Octave-based data:
- PSD spectrum
- SPL spectrum
SourceFiles
noiseFFT.C
SeeAlso
windowModel.H
\*---------------------------------------------------------------------------*/
#ifndef noiseFFT_H
@ -37,6 +52,7 @@ SourceFiles
#include "scalarField.H"
#include "graph.H"
#include "windowModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,7 +60,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class noiseFFT Declaration
Class noiseFFT Declaration
\*---------------------------------------------------------------------------*/
class noiseFFT
@ -54,7 +70,18 @@ class noiseFFT
// Private data
//- Time spacing of the raw data
scalar deltat_;
scalar deltaT_;
struct octaveBandInfo
{
label octave;
// IDs of bin boundaries in pressure data
labelList binIDs;
// Centre frequencies for each bin
scalarField centreFreq;
};
public:
@ -81,37 +108,61 @@ public:
// Member Functions
//- Return the graph of p(t)
//- Return the FFT frequencies
static tmp<scalarField> frequencies
(
const label N,
const scalar deltaT
);
//- Return the PSD [dB/Hz]
// Input PSD in [Pa^2/Hz]
static tmp<scalarField> PSD(const scalarField& PSDf);
//- Return the SPL [dB]
// Input P(rms)^2 in [Pa^2]
static tmp<scalarField> SPL(const scalarField& Prms2);
//- Return a list of the frequency indices wrt f field that
// correspond to the bands limits for a given octave
static void octaveBandInfo
(
const scalarField& f,
const scalar fLower,
const scalar fUpper,
const scalar octave,
labelList& fBandIDs,
scalarField& fCentre
);
//- Return the graph of pressure as a function of time
graph pt() const;
//- Return the nth window
tmp<scalarField> window(const label N, const label n) const;
//- Return the Hanning window function
tmp<scalarField> Hanning(const label N) const;
//- Return the fft of the given pressure data
tmp<scalarField> Pf(const tmp<scalarField>& pn) const;
//- Return the multi-window mean fft of the complete pressure data
graph meanPf(const label N, const label nw) const;
//- Return the multi-window mean fft of the complete pressure data [Pa]
graph meanPf(const windowModel& window) const;
//- Return the multi-window RMS mean fft of the complete pressure data
graph RMSmeanPf(const label N, const label nw) const;
//- Return the multi-window RMS mean fft of the complete pressure
// data [Pa]
graph RMSmeanPf(const windowModel& window) const;
//- Return the narrow-band PFL (pressure-fluctuation level) spectrum
graph Lf(const graph& gPf) const;
//- Return the multi-window PSD (power spectral density) of the complete
// pressure data [Pa^2/Hz]
graph PSDf(const windowModel& window) const;
//- Return the one-third-octave-band PFL spectrum
// starting at octave with mean frequency f1
graph Ldelta(const graph& gLf, const scalar f1, const scalar fU) const;
//- Return the PSD [dB/Hz]
// Takes PSD in [Pa^2/Hz]
graph PSD(const graph& gPSDf) const;
//- Return the one-third-octave-band pressure spectrum
// starting at octave with mean frequency f1
graph Pdelta(const graph& gLf, const scalar f1, const scalar fU) const;
//- Return the total PFL as the sum of Lf over all frequencies
scalar Lsum(const graph& gLf) const;
//- Generate octave data
graph octaves
(
const graph& g,
const labelList& freqBandIDs,
bool integrate
) const;
//- Convert the db into Pa
scalar dbToPa(const scalar db) const;

151
src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C Normal file
View File

@ -0,0 +1,151 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "noiseModel.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(noiseModel, 0);
defineRunTimeSelectionTable(noiseModel, dictionary);
}
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
Foam::scalar Foam::noiseModel::checkUniformTimeStep
(
const scalarList& times
) const
{
scalar deltaT = -1.0;
if (times.size() > 1)
{
for (label i = 1; i < times.size(); i++)
{
scalar dT = times[i] - times[i-1];
if (deltaT < 0)
{
deltaT = dT;
}
if (mag(deltaT - dT) > SMALL)
{
FatalErrorInFunction
<< "Unable to process data with a variable time step"
<< exit(FatalError);
}
}
}
else
{
FatalErrorInFunction
<< "Unable to create FFT with a single value"
<< exit(FatalError);
}
return deltaT;
}
Foam::label Foam::noiseModel::findStartTimeIndex
(
const instantList& allTimes,
const scalar startTime
) const
{
forAll(allTimes, timeI)
{
const instant& i = allTimes[timeI];
if (i.value() >= startTime)
{
return timeI;
}
}
return 0;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::noiseModel::noiseModel(const dictionary& dict)
:
dict_(dict),
rhoRef_(dict.lookupOrDefault<scalar>("rhoRef", 1)),
nSamples_(dict.lookupOrDefault<label>("N", 65536)),
fLower_(dict.lookupOrDefault<scalar>("fl", 25)),
fUpper_(dict.lookupOrDefault<scalar>("fu", 10000)),
startTime_(dict.lookupOrDefault<scalar>("startTime", 0)),
windowModelPtr_(windowModel::New(dict, nSamples_)),
graphFormat_(dict.lookupOrDefault<word>("graphFormat", "raw"))
{
// Check number of samples - must be a power of 2 for our FFT
bool powerOf2 = ((nSamples_ != 0) && !(nSamples_ & (nSamples_ - 1)));
if (!powerOf2)
{
FatalIOErrorInFunction(dict)
<< "N: Number of samples in sampling windows must be a "
<< "power of 2"
<< exit(FatalIOError);
}
if (fLower_ < 0)
{
FatalIOErrorInFunction(dict)
<< "fl: lower frequency bound must be greater than zero"
<< exit(FatalIOError);
}
if (fUpper_ < 0)
{
FatalIOErrorInFunction(dict)
<< "fu: upper frequency bound must be greater than zero"
<< exit(FatalIOError);
}
if (fUpper_ < fLower_)
{
FatalIOErrorInFunction(dict)
<< "fu: upper frequency bound must be greater than lower "
<< "frequency bound (fl)"
<< exit(FatalIOError);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::noiseModel::~noiseModel()
{}
// ************************************************************************* //

179
src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H Normal file
View File

@ -0,0 +1,179 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::noiseModel
Description
Base class for noise models.
Data is read from a dictionary, e.g.
\verbatim
rhoRef 0;
N 4096;
fl 25;
fu 25;
startTime 0;
\endverbatim
where
\table
Property | Description | Required | Default value
rhoRef | Reference density | no | 1
N | Number of samples in sampling window | no | 65536 (2^16)
fl | Lower frequency bounds | no | 25
fu | Upper frequency bounds | no | 10000
startTime | Start time | no | 0
graphFormat | Graph format | no | raw
\endtable
Note
The number of samples in the sampling window must be a power of 2
SourceFiles
noiseModel.C
\*---------------------------------------------------------------------------*/
#ifndef noiseModel_H
#define noiseModel_H
#include "dictionary.H"
#include "scalarList.H"
#include "instantList.H"
#include "windowModel.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class noiseModel Declaration
\*---------------------------------------------------------------------------*/
class noiseModel
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
noiseModel(const noiseModel&);
//- Disallow default bitwise assignment
void operator=(const noiseModel&);
protected:
// Protected Data
//- Copy of dictionary used for construction
const dictionary dict_;
//- Reference density (to convert from kinematic to static pressure)
scalar rhoRef_;
//- Number of samples in sampling window, default = 2^16
label nSamples_;
//- Lower frequency limit, default = 25Hz
scalar fLower_;
//- Upper frequency limit, default = 10kHz
scalar fUpper_;
//- Start time, default = 0s
scalar startTime_;
//- Window model
autoPtr<windowModel> windowModelPtr_;
//- Graph format
word graphFormat_;
// Protected Member Functions
//- Check and return uniform time step
scalar checkUniformTimeStep
(
const scalarList& times
) const;
//- Find and return start time index
label findStartTimeIndex
(
const instantList& allTimes,
const scalar startTime
) const;
public:
//- Runtime type information
TypeName("noiseModel");
//- Run time selection table
declareRunTimeSelectionTable
(
autoPtr,
noiseModel,
dictionary,
(
const dictionary& dict
),
(dict)
);
//- Selector
static autoPtr<noiseModel> New(const dictionary& dict);
//- Constructor
noiseModel(const dictionary& dict);
//- Destructor
virtual ~noiseModel();
// Public Member Functions
//- Abstract call to calculate
virtual void calculate() = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

53
src/randomProcesses/noise/noiseModels/noiseModel/noiseModelNew.C Normal file
View File

@ -0,0 +1,53 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "noiseModel.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::noiseModel> Foam::noiseModel::New(const dictionary& dict)
{
const word modelType(dict.lookup("noiseModel"));
Info<< "Selecting noiseModel " << modelType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(modelType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorInFunction
<< "Unknown noiseModel type "
<< modelType << nl << nl
<< "Valid noiseModel types are:" << nl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<noiseModel>(cstrIter()(dict.subDict(modelType + "Coeffs")));
}
// ************************************************************************* //

207
src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C Normal file
View File

@ -0,0 +1,207 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "pointNoise.H"
#include "noiseFFT.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace noiseModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(pointNoise, 0);
addToRunTimeSelectionTable(noiseModel, pointNoise, dictionary);
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void pointNoise::filterTimeData
(
const Function1Types::CSV<scalar>& pData,
scalarField& t,
scalarField& p
)
{
const scalarField t0(pData.x());
const scalarField p0(pData.y());
DynamicList<scalar> tf(t0.size());
DynamicList<scalar> pf(t0.size());
forAll(t0, timeI)
{
if (t0[timeI] >= startTime_)
{
tf.append(t0[timeI]);
pf.append(p0[timeI]);
}
}
t.transfer(tf);
p.transfer(pf);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void pointNoise::calculate()
{
// Point data only handled by master
if (!Pstream::master())
{
return;
}
Info<< "Reading data file" << endl;
Function1Types::CSV<scalar> pData("pressure", dict_, "Data");
// Time and pressure history data
scalarField t, p;
filterTimeData(pData, t, p);
p *= rhoRef_;
Info<< " read " << t.size() << " values" << nl << endl;
Info<< "Creating noise FFT" << endl;
const scalar deltaT = checkUniformTimeStep(t);
// Determine the windowing
windowModelPtr_->validate(t.size());
const windowModel& win = windowModelPtr_();
const scalar deltaf = 1.0/(deltaT*win.nSamples());
fileName outDir(fileName("postProcessing")/"noise"/typeName);
// Create the fft
noiseFFT nfft(deltaT, p);
// Narrow band data
// ----------------
// RMS pressure [Pa]
graph Prmsf(nfft.RMSmeanPf(win));
Info<< " Creating graph for " << Prmsf.title() << endl;
Prmsf.write(outDir, graph::wordify(Prmsf.title()), graphFormat_);
// PSD [Pa^2/Hz]
graph PSDf(nfft.PSDf(win));
Info<< " Creating graph for " << PSDf.title() << endl;
PSDf.write(outDir, graph::wordify(PSDf.title()), graphFormat_);
// PSD [dB/Hz]
graph PSDg
(
"PSD_dB_Hz(f)",
"f [Hz]",
"PSD(f) [dB_Hz]",
Prmsf.x(),
noiseFFT::PSD(PSDf.y())
);
Info<< " Creating graph for " << PSDg.title() << endl;
PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_);
// SPL [dB]
graph SPLg
(
"SPL_dB(f)",
"f [Hz]",
"SPL(f) [dB]",
Prmsf.x(),
noiseFFT::SPL(PSDf.y()*deltaf)
);
Info<< " Creating graph for " << SPLg.title() << endl;
SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
labelList octave13BandIDs;
scalarField octave13FreqCentre;
noiseFFT::octaveBandInfo
(
Prmsf.x(),
fLower_,
fUpper_,
3,
octave13BandIDs,
octave13FreqCentre
);
// 1/3 octave data
// ---------------
// PSD [Pa^2/Hz]
graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false));
// Integrated PSD = P(rms)^2 [Pa^2]
graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true));
graph PSD13g
(
"PSD13_dB_Hz(fm)",
"fm [Hz]",
"PSD(fm) [dB_Hz]",
octave13FreqCentre,
noiseFFT::PSD(PSD13f.y())
);
Info<< " Creating graph for " << PSD13g.title() << endl;
PSD13g.write(outDir, graph::wordify(PSD13g.title()), graphFormat_);
graph SPL13g
(
"SPL13_dB(fm)",
"fm [Hz]",
"SPL(fm) [dB]",
octave13FreqCentre,
noiseFFT::SPL(Prms13f2.y())
);
Info<< " Creating graph for " << SPL13g.title() << endl;
SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
pointNoise::pointNoise(const dictionary& dict)
:
noiseModel(dict)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
pointNoise::~pointNoise()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace noiseModels
} // End namespace Foam
// ************************************************************************* //

142
src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H Normal file
View File

@ -0,0 +1,142 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::noiseModels::pointNoise
Description
Perform noise analysis on point-based pressure data.
Input data is read from a dictionary, e.g.
\verbatim
// Pressure reference
pRef 0;
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096;
// Lower frequency bounds
fl 25;
// Upper frequency bounds
fu 25;
// Start time
startTime 0;
windowModel <modelType>
<modelType>Coeffs
{
...
}
// Pressure data supplied in CSV file format
csvFileData
{
fileName "pressureData";
nHeaderLine 1;
refColumn 0;
componentColumns (1);
separator " ";
mergeSeparators yes;
}
graphFormat raw;
\endverbatim
SourceFiles
pointNoise.C
SeeAlso
noiseModel.H
windowModel.H
\*---------------------------------------------------------------------------*/
#ifndef pointNoise_H
#define pointNoise_H
#include "noiseModel.H"
#include "CSV.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace noiseModels
{
/*---------------------------------------------------------------------------*\
Class pointNoise Declaration
\*---------------------------------------------------------------------------*/
class pointNoise
:
public noiseModel
{
protected:
// Protected Member Functions
void filterTimeData
(
const Function1Types::CSV<scalar>& pData,
scalarField& t,
scalarField& p
);
public:
//- Runtime type information
TypeName("pointNoise");
//- Constructor
pointNoise(const dictionary& dict);
//- Destructor
virtual ~pointNoise();
// Public Member Functions
//- Calculate
virtual void calculate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace noiseModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

705
src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C Normal file
View File

@ -0,0 +1,705 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "surfaceNoise.H"
#include "surfaceReader.H"
#include "surfaceWriter.H"
#include "noiseFFT.H"
#include "graph.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace noiseModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(surfaceNoise, 0);
addToRunTimeSelectionTable(noiseModel, surfaceNoise, dictionary);
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void surfaceNoise::initialise(const dictionary& dict)
{
label nAvailableTimes = 0;
// All reading performed on the master processor only
if (Pstream::master())
{
// Create the surface reader
const word readerType(dict.lookup("reader"));
readerPtr_.reset(surfaceReader::New(readerType, inputFileName_).ptr());
// Find the index of the pressure data
const word pName(dict.lookupOrDefault<word>("pName", "p"));
const List<word> fieldNames(readerPtr_->fieldNames(0));
pIndex_ = findIndex(fieldNames, pName);
if (pIndex_ == -1)
{
FatalErrorInFunction
<< "Unable to find pressure field name " << pName
<< " in list of available fields: " << fieldNames
<< exit(FatalError);
}
// Create the surface writer
// - Could be done later, but since this utility can process a lot of
// data we can ensure that the user-input is correct prior to doing
// the heavy lifting
const word writerType(dict.lookup("writer"));
dictionary optDict
(
dict.subOrEmptyDict("writeOptions").subOrEmptyDict(writerType)
);
writerPtr_.reset(surfaceWriter::New(writerType, optDict).ptr());
// Set the time range
const instantList allTimes = readerPtr_->times();
startTimeIndex_ = findStartTimeIndex(allTimes, startTime_);
// Determine the windowing
nAvailableTimes = allTimes.size() - startTimeIndex_;
}
Pstream::scatter(pIndex_);
Pstream::scatter(startTimeIndex_);
Pstream::scatter(nAvailableTimes);
// Note: all processors should call the windowing validate function
label nRequiredTimes = windowModelPtr_->validate(nAvailableTimes);
if (Pstream::master())
{
// Restrict times
const instantList allTimes = readerPtr_->times();
times_.setSize(nRequiredTimes);
forAll(times_, timeI)
{
times_[timeI] = allTimes[timeI + startTimeIndex_].value();
}
deltaT_ = checkUniformTimeStep(times_);
// Read the surface geometry
const meshedSurface& surf = readerPtr_->geometry();
nFace_ = surf.size();
}
Pstream::scatter(times_);
Pstream::scatter(deltaT_);
Pstream::scatter(nFace_);
}
void surfaceNoise::readSurfaceData
(
const labelList& procFaceOffset,
List<scalarField>& pData
)
{
// Data is stored as pressure on surface at a given time. Now we need to
// pivot the data so that we have the complete pressure time history per
// surface face. In serial mode, this results in all pressure data being
// loaded into memory (!)
Info << "Reading pressure data" << endl;
if (Pstream::parRun())
{
PstreamBuffers pBufs(Pstream::nonBlocking);
// Procedure:
// 1. Master processor reads pressure data for all faces for all times
// 2. Master sends each processor a subset of faces
// 3. Local processor reconstructs the full time history for the subset
// of faces
// Note: reading all data on master to avoid potential NFS problems...
const label myProcI = Pstream::myProcNo();
const label nLocalFace =
procFaceOffset[myProcI + 1] - procFaceOffset[myProcI];
// Complete pressure time history data for subset of faces
pData.setSize(nLocalFace);
const label nTimes = times_.size();
forAll(pData, faceI)
{
pData[faceI].setSize(nTimes);
}
// Read and send data
forAll(times_, i)
{
pBufs.clear();
if (Pstream::master())
{
label timeI = i + startTimeIndex_;
Info<< " time: " << times_[i] << endl;
// Read pressure at all faces for time timeI
scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
// Apply conversions
p *= rhoRef_;
// Send subset of faces to each processor
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
label face0 = procFaceOffset[procI];
label nLocalFace =
procFaceOffset[procI + 1] - procFaceOffset[procI];
UOPstream toProc(procI, pBufs);
toProc << SubList<scalar>(p, nLocalFace, face0);
}
}
pBufs.finishedSends();
// Receive data from the master
UIPstream fromProc(0, pBufs);
scalarList pSlice(fromProc);
forAll(pSlice, faceI)
{
pData[faceI][i] = pSlice[faceI];
}
}
}
else
{
const label nLocalFace = procFaceOffset[0];
pData.setSize(nLocalFace);
forAll(times_, timeI)
{
forAll(pData, faceI)
{
pData[faceI].setSize(times_.size());
}
}
forAll(times_, i)
{
label timeI = i + startTimeIndex_;
Info<< " time: " << times_[i] << endl;
const scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
forAll(p, faceI)
{
pData[faceI][i] = p[faceI]*rhoRef_;
}
}
}
Info<< "Read "
<< returnReduce(pData.size(), sumOp<label>())
<< " pressure traces with "
<< times_.size()
<< " time values" << nl << endl;
}
Foam::scalar surfaceNoise::writeSurfaceData
(
const word& fName,
const word& groupName,
const word& title,
const scalar freq,
const scalarField& data,
const labelList& procFaceOffset
) const
{
Info<< " processing " << title << " for frequency " << freq << endl;
fileName outDir
(
fileName("postProcessing")/"noise"/groupName/Foam::name(freq)
);
if (Pstream::parRun())
{
// Collect the surface data so that we can output the surfaces
PstreamBuffers pBufs(Pstream::nonBlocking);
if (!Pstream::master())
{
UOPstream toProc(0, pBufs);
toProc << data;
}
pBufs.finishedSends();
scalar areaAverage = 0;
if (Pstream::master())
{
const meshedSurface& surf = readerPtr_->geometry();
scalarField allData(surf.size());
forAll(data, faceI)
{
// Master procFaceOffset is zero...
allData[faceI] = data[faceI];
}
for (label procI = 1; procI < Pstream::nProcs(); procI++)
{
UIPstream fromProc(procI, pBufs);
scalarList dataSlice(fromProc);
forAll(dataSlice, i)
{
label faceI = procFaceOffset[procI] + i;
allData[faceI] = dataSlice[i];
}
}
fileName outFileName = writerPtr_->write
(
outDir,
fName,
surf.points(),
surf.faces(),
title,
allData,
false
);
// TODO: Move faceAreas to demand-driven function in MeshedSurface
// scalarField faceAreas(surf.faces().size());
// forAll(faceAreas, i)
// {
// faceAreas[i] = surf.faces()[i].mag(surf.points());
// }
//
// areaAverage = sum(allData*faceAreas)/sum(faceAreas);
areaAverage = sum(allData)/allData.size();
}
Pstream::scatter(areaAverage);
return areaAverage;
}
else
{
const meshedSurface& surf = readerPtr_->geometry();
writerPtr_->write
(
outDir,
fName,
surf.points(),
surf.faces(),
title,
data,
false
);
// TODO: Move faceAreas to demand-driven function in MeshedSurface
// scalarField faceAreas(surf.faces().size());
// forAll(faceAreas, i)
// {
// faceAreas[i] = surf.faces()[i].mag(surf.points());
// }
//
// return sum(data*faceAreas)/sum(faceAreas);
return sum(data)/data.size();
}
}
Foam::scalar surfaceNoise::surfaceAverage
(
const scalarField& data,
const labelList& procFaceOffset
) const
{
if (Pstream::parRun())
{
// Collect the surface data so that we can output the surfaces
PstreamBuffers pBufs(Pstream::nonBlocking);
if (!Pstream::master())
{
UOPstream toProc(0, pBufs);
toProc << data;
}
pBufs.finishedSends();
scalar areaAverage = 0;
if (Pstream::master())
{
const meshedSurface& surf = readerPtr_->geometry();
scalarField allData(surf.size());
forAll(data, faceI)
{
// Master procFaceOffset is zero...
allData[faceI] = data[faceI];
}
for (label procI = 1; procI < Pstream::nProcs(); procI++)
{
UIPstream fromProc(procI, pBufs);
scalarList dataSlice(fromProc);
forAll(dataSlice, i)
{
label faceI = procFaceOffset[procI] + i;
allData[faceI] = dataSlice[i];
}
}
// TODO: Move faceAreas to demand-driven function in MeshedSurface
scalarField faceAreas(surf.faces().size());
forAll(faceAreas, i)
{
faceAreas[i] = surf.faces()[i].mag(surf.points());
}
// areaAverage = sum(allData*faceAreas)/sum(faceAreas);
areaAverage = sum(allData)/allData.size();
}
Pstream::scatter(areaAverage);
return areaAverage;
}
else
{
const meshedSurface& surf = readerPtr_->geometry();
// TODO: Move faceAreas to demand-driven function in MeshedSurface
scalarField faceAreas(surf.faces().size());
forAll(faceAreas, i)
{
faceAreas[i] = surf.faces()[i].mag(surf.points());
}
// return sum(data*faceAreas)/sum(faceAreas);
return sum(data)/data.size();
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
surfaceNoise::surfaceNoise(const dictionary& dict)
:
noiseModel(dict),
inputFileName_(dict.lookup("inputFile")),
pIndex_(0),
times_(),
deltaT_(0),
startTimeIndex_(0),
nFace_(0),
fftWriteInterval_(dict.lookupOrDefault("fftWriteInterval", 1))
{
initialise(dict);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
surfaceNoise::~surfaceNoise()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void surfaceNoise::calculate()
{
// Container for pressure time history data per face
List<scalarField> pData;
// Processor procFaceOffsets
labelList procFaceOffset;
if (Pstream::parRun())
{
const label nProcs = Pstream::nProcs();
const label nFacePerProc = floor(nFace_/nProcs) + 1;
procFaceOffset.setSize(nProcs + 1, 0);
for (label i = 1; i < procFaceOffset.size(); i++)
{
procFaceOffset[i] = min(i*nFacePerProc, nFace_);
}
}
else
{
procFaceOffset.setSize(1, nFace_);
}
// Read pressure data from file
readSurfaceData(procFaceOffset, pData);
// Process the pressure data, and store results as surface values per
// frequency so that it can be output using the surface writer
Info<< "Creating noise FFTs" << endl;
// Storage for FFT data
const label nLocalFace = pData.size();
const scalarField freq1(noiseFFT::frequencies(nSamples_, deltaT_));
const label nFFT = freq1.size()/fftWriteInterval_;
List<scalarField> surfPrmsf(nFFT);
List<scalarField> surfPSDf(nFFT);
forAll(surfPrmsf, freqI)
{
surfPrmsf[freqI].setSize(nLocalFace);
surfPSDf[freqI].setSize(nLocalFace);
}
// Storage for 1/3 octave data
labelList octave13BandIDs;
scalarField octave13FreqCentre;
noiseFFT::octaveBandInfo
(
freq1,
fLower_,
fUpper_,
3,
octave13BandIDs,
octave13FreqCentre
);
List<scalarField> surfPSD13f(octave13BandIDs.size() - 1);
List<scalarField> surfPrms13f2(octave13BandIDs.size() - 1);
forAll(surfPSD13f, freqI)
{
surfPSD13f[freqI].setSize(nLocalFace);
surfPrms13f2[freqI].setSize(nLocalFace);
}
const windowModel& win = windowModelPtr_();
forAll(pData, faceI)
{
const scalarField& p = pData[faceI];
noiseFFT nfft(deltaT_, p);
graph Prmsf(nfft.RMSmeanPf(win));
graph PSDf(nfft.PSDf(win));
// Store the frequency results in slot for face of surface
forAll(surfPrmsf, i)
{
label freqI = (i + 1)*fftWriteInterval_ - 1;
surfPrmsf[i][faceI] = Prmsf.y()[freqI];
surfPSDf[i][faceI] = PSDf.y()[freqI];
}
// PSD [Pa^2/Hz]
graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false));
// Integrated PSD = P(rms)^2 [Pa^2]
graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true));
// Store the 1/3 octave results in slot for face of surface
forAll(surfPSD13f, freqI)
{
surfPSD13f[freqI][faceI] = PSD13f.y()[freqI];
surfPrms13f2[freqI][faceI] = Prms13f2.y()[freqI];
}
// Free the storage for p
// p.clear();
}
// Output directory for graphs
fileName outDir(fileName("postProcessing")/"noise"/typeName);
const scalar deltaf = 1.0/(deltaT_*win.nSamples());
Info<< "Writing fft surface data" << endl;
{
scalarField PrmsfAve(surfPrmsf.size(), 0);
scalarField PSDfAve(surfPrmsf.size(), 0);
scalarField fOut(surfPrmsf.size(), 0);
forAll(surfPrmsf, i)
{
label freqI = i*fftWriteInterval_;
fOut[i] = freq1[freqI];
const word& fName = inputFileName_.name(true);
const word gName = "fft";
PrmsfAve[i] = writeSurfaceData
(
fName,
gName,
"Prmsf",
freq1[freqI],
surfPrmsf[i],
procFaceOffset
);
PSDfAve[i] = writeSurfaceData
(
fName,
gName,
"PSDf",
freq1[freqI],
surfPSDf[i],
procFaceOffset
);
writeSurfaceData
(
fName,
gName,
"PSD",
freq1[freqI],
noiseFFT::PSD(surfPSDf[i]),
procFaceOffset
);
writeSurfaceData
(
fName,
gName,
"SPL",
freq1[freqI],
noiseFFT::SPL(surfPSDf[i]*deltaf),
procFaceOffset
);
}
graph Prmsfg
(
"Average Prms(f)",
"f [Hz]",
"P(f) [Pa]",
fOut,
PrmsfAve
);
Prmsfg.write(outDir, graph::wordify(Prmsfg.title()), graphFormat_);
graph PSDfg
(
"Average PSD_f(f)",
"f [Hz]",
"PSD(f) [PaPa_Hz]",
fOut,
PSDfAve
);
PSDfg.write(outDir, graph::wordify(PSDfg.title()), graphFormat_);
graph PSDg
(
"Average PSD_dB_Hz(f)",
"f [Hz]",
"PSD(f) [dB_Hz]",
fOut,
noiseFFT::PSD(PSDfAve)
);
PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_);
graph SPLg
(
"Average SPL_dB(f)",
"f [Hz]",
"SPL(f) [dB]",
fOut,
noiseFFT::SPL(PSDfAve*deltaf)
);
SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
}
Info<< "Writing one-third octave surface data" << endl;
{
scalarField PSDfAve(surfPSD13f.size(), 0);
scalarField Prms13f2Ave(surfPSD13f.size(), 0);
forAll(surfPSD13f, i)
{
const word& fName = inputFileName_.name(true);
const word gName = "oneThirdOctave";
PSDfAve[i] = writeSurfaceData
(
fName,
gName,
"PSD13f",
octave13FreqCentre[i],
surfPSD13f[i],
procFaceOffset
);
writeSurfaceData
(
fName,
gName,
"PSD13",
octave13FreqCentre[i],
noiseFFT::PSD(surfPSD13f[i]),
procFaceOffset
);
writeSurfaceData
(
fName,
gName,
"SPL13",
octave13FreqCentre[i],
noiseFFT::SPL(surfPrms13f2[i]),
procFaceOffset
);
Prms13f2Ave[i] = surfaceAverage(surfPrms13f2[i], procFaceOffset);
}
graph PSD13g
(
"Average PSD13_dB_Hz(fm)",
"fm [Hz]",
"PSD(fm) [dB_Hz]",
octave13FreqCentre,
noiseFFT::PSD(PSDfAve)
);
PSD13g.write(outDir, graph::wordify(PSD13g.title()), graphFormat_);
graph SPL13g
(
"Average SPL13_dB(fm)",
"fm [Hz]",
"SPL(fm) [dB]",
octave13FreqCentre,
noiseFFT::SPL(Prms13f2Ave)
);
SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace noiseModels
} // End namespace Foam
// ************************************************************************* //

208
src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H Normal file
View File

@ -0,0 +1,208 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::noiseModels::surfaceNoise
Description
Perform noise analysis on surface-based pressure data.
Input data is read from a dictionary, e.g.
\verbatim
// Pressure reference
pRef 0;
// Number of samples in sampling window
// Must be a power of 2, default = 2^16 (=65536)
N 4096;
// Lower frequency bounds
fl 25;
// Upper frequency bounds
fu 25;
// Start time
startTime 0;
windowModel <modelType>
<modelType>Coeffs
{
...
}
// Input file
inputFile "postProcessing/faceSource1/surface/patch/patch.case";
// Write interval for FFT data, default = 1
fftWriteInterval 100;
// Surface reader
reader ensight;
// Surface writer
writer ensight;
// Collate times for ensight output - ensures geometry is only written once
writeOptions
{
ensight
{
collateTimes 1;
}
}
\endverbatim
SourceFiles
surfaceNoise.C
SeeAlso
noiseModel.H
\*---------------------------------------------------------------------------*/
#ifndef surfaceNoise_H
#define surfaceNoise_H
#include "noiseModel.H"
#include "labelList.H"
#include "scalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class surfaceReader;
class surfaceWriter;
namespace noiseModels
{
/*---------------------------------------------------------------------------*\
Class surfaceNoise Declaration
\*---------------------------------------------------------------------------*/
class surfaceNoise
:
public noiseModel
{
protected:
// Protected Data
//- Input file name
const fileName inputFileName_;
//- Index of pressure field in reader field list
label pIndex_;
//- Sample times
scalarList times_;
//- Time step (constant)
scalar deltaT_;
//- Start time index
label startTimeIndex_;
//- Number of surface faces
label nFace_;
//- Frequency data output interval, default = 1
// nSamples/2 data points are returned from the FFT, which can
// result in a very large number of output files (1 per frequency)
label fftWriteInterval_;
//- Pointer to the surface reader
mutable autoPtr<surfaceReader> readerPtr_;
//- Pointer to the surface writer
autoPtr<surfaceWriter> writerPtr_;
// Protected Member Functions
//- Initialise
void initialise(const dictionary& dict);
//- Read surface data
void readSurfaceData
(
const labelList& procFaceOffset,
List<scalarField>& pData
);
//- Write surface data to file
// Returns the area average value
scalar writeSurfaceData
(
const word& fName,
const word& groupName,
const word& title,
const scalar freq,
const scalarField& data,
const labelList& procFaceOffset
) const;
//- Calculate the area average value
scalar surfaceAverage
(
const scalarField& data,
const labelList& procFaceOffset
) const;
public:
//- Runtime type information
TypeName("surfaceNoise");
//- Constructor
surfaceNoise(const dictionary& dict);
//- Destructor
virtual ~surfaceNoise();
// Public Member Functions
//- Calculate
virtual void calculate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace noiseModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

122
src/randomProcesses/windowModels/Hanning/Hanning.C Normal file
View File

@ -0,0 +1,122 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "Hanning.H"
#include "addToRunTimeSelectionTable.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace windowModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Hanning, 0);
addToRunTimeSelectionTable(windowModel, Hanning, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Hanning::Hanning(const dictionary& dict, const label nSamples)
:
windowModel(dict, nSamples),
symmetric_(readBool(dict.lookup("symmetric"))),
extended_(readBool(dict.lookup("extended"))),
alpha_(dict.lookupOrDefault("alpha", 0.5)) // Hamming = 0.54
{
// Extend range if required
label offset = extended_ ? 1 : 0;
scalar m = nSamples - 1 + 2*offset;
scalarField t(nSamples);
forAll(t, i)
{
t[i] = i + offset;
}
scalarField& wf = *this;
wf = alpha_ - (1 - alpha_)*cos(constant::mathematical::twoPi*t/m);
// Reset second half of window if symmetric
if (symmetric_)
{
label midPointI = 0;
if (nSamples % 2 == 0)
{
midPointI = nSamples/2;
}
else
{
midPointI = (nSamples + 1)/2;
}
for (label i = 0; i < midPointI; i++)
{
wf[nSamples - i - 1] = wf[i];
}
}
scalar sumSqr = sum(sqr(wf));
// Normalisation
wf *= sqrt(nSamples/sumSqr);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Hanning::~Hanning()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Hanning::symmetric() const
{
return symmetric_;
}
bool Hanning::extended() const
{
return extended_;
}
Foam::scalar Hanning::alpha() const
{
return alpha_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace windowModels
} // End namespace Foam
// ************************************************************************* //

125
src/randomProcesses/windowModels/Hanning/Hanning.H Normal file
View File

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::windowModels::Hanning
Description
Hanning window
The window is described by the function
\f[
wf = (1 - \alpha) - \alpha cos(2 \pi t/m);
\f]
Where:
\vartable
\alpha | Coefficient with a default value of 0.5
t | time
m | window width
\endvartable
The window can be further manipulated by the controls:
- \c symmetric: force the window to be symmetric
- \c extended: extend the window by 1 element at start and end to produce
non-zero values at the start and end positions. Note: window is
normalised to preserve energy content
SourceFiles
Hanning.C
\*---------------------------------------------------------------------------*/
#ifndef Hanning_H
#define Hanning_H
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "windowModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace windowModels
{
/*---------------------------------------------------------------------------*\
Class Hanning Declaration
\*---------------------------------------------------------------------------*/
class Hanning
:
public windowModel
{
protected:
// Protected Member Data
//- Symmetric switch
bool symmetric_;
//- Extended switch
bool extended_;
//- Window coefficient, default = 0.5
scalar alpha_;
public:
//- Runtime type information
TypeName("Hanning");
//- Construct from dictionary
Hanning(const dictionary& dict, const label nSamples);
//- Destuctor
virtual ~Hanning();
// Public Member Functions
//- Return the symmetric flag
bool symmetric() const;
//- Return the extended flag
bool extended() const;
//- Return the window coefficient
scalar alpha() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace windowModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

66
src/randomProcesses/windowModels/uniform/uniform.C Normal file
View File

@ -0,0 +1,66 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "uniform.H"
#include "addToRunTimeSelectionTable.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace windowModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(uniform, 0);
addToRunTimeSelectionTable(windowModel, uniform, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
uniform::uniform(const dictionary& dict, const label nSamples)
:
windowModel(dict, nSamples),
value_(readScalar(dict.lookup("value")))
{
scalarField& wf = *this;
wf = value_;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
uniform::~uniform()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace windowModels
} // End namespace Foam
// ************************************************************************* //

70
src/randomProcesses/fft/fftRenumber.H → src/randomProcesses/windowModels/uniform/uniform.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,57 +21,65 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
fftRenumber
Class
Foam::windowModels::uniform
Description
Multi-dimensional renumbering used in the Numerical Recipes
fft routine.
A window that applies uniform scaling.
SourceFiles
fftRenumber.C
uniform.C
\*---------------------------------------------------------------------------*/
#ifndef fftRenumber_H
#define fftRenumber_H
#ifndef uniform_H
#define uniform_H
#include "complex.H"
#include "List.H"
#include "labelList.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "windowModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace windowModels
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
/*---------------------------------------------------------------------------*\
Class uniform Declaration
\*---------------------------------------------------------------------------*/
// Recursively evaluate the indexing necessary to do the folding of the fft
// data. We recurse until we have the indexing ncessary for the folding in all
// directions.
void fftRenumberRecurse
(
List<complex>& data,
List<complex>& renumData,
const labelList& nn,
label nnprod,
label ii,
label l1,
label l2
);
class uniform
:
public windowModel
{
protected:
// Protected data
//- Uniform value
scalar value_;
// Fold the n-d data array to get the fft components in the right places.
void fftRenumber
(
List<complex>& data,
const labelList& nn
);
public:
//- Runtime type information
TypeName("uniform");
//- Construct from dictionary
uniform(const dictionary& dict, const label nSamples);
//- Destuctor
virtual ~uniform();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace windowModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

125
src/randomProcesses/windowModels/windowModel/windowModel.C Normal file
View File

@ -0,0 +1,125 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "windowModel.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(windowModel, 0);
defineRunTimeSelectionTable(windowModel, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::windowModel::windowModel(const dictionary& dict, const label nSamples)
:
scalarField(nSamples),
nOverlapSamples_(0),
nWindow_(dict.lookupOrDefault("nWindow", -1))
{
scalar prc = readScalar(dict.lookup("overlapPercent"));
nOverlapSamples_ = floor(prc/scalar(100)*nSamples);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::windowModel::~windowModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::windowModel::nSamples() const
{
return size();
}
Foam::label Foam::windowModel::nWindow() const
{
return nWindow_;
}
Foam::label Foam::windowModel::nWindowsTotal(label nSamplesTotal) const
{
const label nSamples = this->nSamples();
return floor((nSamplesTotal - nSamples)/(nSamples - nOverlapSamples_)) + 1;
}
Foam::label Foam::windowModel::validate(const label nSamplesTotal)
{
label nSamples = this->nSamples();
if (nSamplesTotal < nSamples)
{
FatalErrorInFunction
<< "Block size N = " << nSamples
<< " is larger than total number of data points = " << nSamplesTotal
<< exit(FatalError);
}
label nWindowAvailable = nWindowsTotal(nSamplesTotal);
if (nWindow_ == -1)
{
nWindow_ = nWindowAvailable;
}
if (nWindow_ > nWindowAvailable)
{
FatalErrorInFunction
<< "Number of data points calculated with " << nWindow_
<< " windows greater than the total number of data points"
<< nl
<< " Block size, N = " << nSamples << nl
<< " Total number of data points = " << nSamplesTotal << nl
<< " Maximum number of windows = " << nWindowAvailable << nl
<< " Requested number of windows = " << nWindow_
<< exit(FatalError);
}
label nRequiredSamples =
nWindow_*nSamples - (nWindow_ - 1)*nOverlapSamples_;
Info<< "Windowing:" << nl
<< " Total samples : " << nSamplesTotal << nl
<< " Samples per window : " << nSamples << nl
<< " Number of windows : " << nWindow_ << nl
<< " Overlap size : " << nOverlapSamples_ << nl
<< " Required number of samples : " << nRequiredSamples
<< endl;
return nRequiredSamples;
}
// ************************************************************************* //

143
src/randomProcesses/windowModels/windowModel/windowModel.H Normal file
View File

@ -0,0 +1,143 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 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::windowModel
Description
Base class for windowing models
SourceFiles
noiseFFT.C
\*---------------------------------------------------------------------------*/
#ifndef windowModel_H
#define windowModel_H
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "scalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class windowModel Declaration
\*---------------------------------------------------------------------------*/
class windowModel
:
public scalarField
{
protected:
// Protected Data
//- Number of overlap samples per window
label nOverlapSamples_;
//- Number of windows
label nWindow_;
public:
//- Runtime type information
TypeName("windowModel");
// Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
windowModel,
dictionary,
(
const dictionary& dict,
const label nSamples
),
(dict, nSamples)
);
//- Construct from dictionary
windowModel(const dictionary& dict, const label nSamples);
// Selectors
//- Return a reference to the selected window model
static autoPtr<windowModel> New
(
const dictionary& dict,
const label nSamples
);
//- Destuctor
virtual ~windowModel();
// Public Member Functions
//- Return the number of samples in the window
label nSamples() const;
//- Return the number of windows
label nWindow() const;
//- Return the total number of windows for a given number of samples
label nWindowsTotal(label nSamplesTotal) const;
//- Validate that the window is applicable to the data set size, and
// return the number of required data points
label validate(label n);
//- Return the windowed data
template<class Type>
tmp<Field<Type> > apply
(
const Field<Type>& fld,
const label windowI
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "windowModelTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

63
src/randomProcesses/windowModels/windowModel/windowModelNew.C Normal file
View File

@ -0,0 +1,63 @@
#/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "windowModel.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::windowModel> Foam::windowModel::New
(
const dictionary& dict,
const label nSamples
)
{
const word modelType(dict.lookup("windowModel"));
Info<< "Selecting windowModel " << modelType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(modelType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"windowModel::New(const dictionary&, const label)"
) << "Unknown windowModel type "
<< modelType << nl << nl
<< "Valid windowModel types are:" << nl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<windowModel>
(
cstrIter()(dict.subDict(modelType + "Coeffs"), nSamples)
);
}
// ************************************************************************* //

68
src/randomProcesses/windowModels/windowModel/windowModelTemplates.C Normal file
View File

@ -0,0 +1,68 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::windowModel::apply
(
const Field<Type>& fld,
const label windowI
) const
{
label nSamples = this->nSamples();
if (nSamples > fld.size())
{
FatalErrorInFunction
<< "Number of samples in sampling window is greater than the "
<< "size of the input field" << nl
<< " input field size = " << fld.size() << nl
<< " window size = " << nSamples << nl
<< " requested window index = " << windowI
<< exit(FatalError);
}
tmp<Field<Type> > tresult(new Field<Type>(nSamples, pTraits<Type>::zero));
Field<Type>& result = tresult.ref();
label nWindow = nWindowsTotal(fld.size());
if (windowI >= nWindow)
{
FatalErrorInFunction
<< "Requested window " << windowI << " outside of range. "
<< "Number of available windows is " << nWindow
<< abort(FatalError);
}
label windowOffset = windowI*(nSamples - nOverlapSamples_);
const scalarField& wf = *this;
result = wf*SubField<Type>(fld, nSamples, windowOffset);
return tresult;
}
// ************************************************************************* //

6
src/sampling/Make/files
View File

@ -51,6 +51,12 @@ $(surfWriters)/starcd/starcdSurfaceWriter.C
$(surfWriters)/vtk/vtkSurfaceWriter.C
$(surfWriters)/boundaryData/boundaryDataSurfaceWriter.C
surfReaders = sampledSurface/readers
$(surfReaders)/surfaceReader.C
$(surfReaders)/surfaceReaderNew.C
$(surfReaders)/ensight/ensightSurfaceReader.C
graphField/writePatchGraph.C
graphField/writeCellGraph.C
graphField/makeGraph.C

515
src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C Normal file
View File

@ -0,0 +1,515 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "ensightSurfaceReader.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(ensightSurfaceReader, 0);
addToRunTimeSelectionTable(surfaceReader, ensightSurfaceReader, fileName);
}
void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
{
label i = 0;
token t;
while (is.good() && (i < n))
{
is >> t;
i++;
if (debug)
{
Info<< "Skipping token " << t << endl;
}
}
if (i != n)
{
WarningInFunction
<< "Requested to skip " << n << "tokens, but stream exited after "
<< i << " tokens. Last token read: " << t
<< endl;
}
}
void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
{
// Binary flag string if applicable
is.readBinaryHeader();
string buffer;
// Ensight Geometry File
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
// Description - 1
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
// Node info
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
// Element info
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
// Part
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
// Part number
label ibuffer;
is.read(ibuffer);
if (debug) Info<< "ibuffer: " << ibuffer << endl;
// Description - 2
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
// Co-ordinates
is.read(buffer);
if (debug) Info<< "buffer: " << buffer << endl;
}
void Foam::ensightSurfaceReader::debugSection
(
const word& expected,
IFstream& is
) const
{
word actual(is);
if (expected != actual)
{
FatalIOErrorInFunction(is)
<< "Expected section header '" << expected
<< "' but read the word '" << actual << "'"
<< exit(FatalIOError);
}
if (debug)
{
Info<< "Read section header: " << expected << endl;
}
}
void Foam::ensightSurfaceReader::readCase(IFstream& is)
{
if (debug)
{
InfoInFunction<< endl;
}
if (!is.good())
{
FatalErrorInFunction
<< "Cannot read file " << is.name()
<< exit(FatalError);
}
string buffer;
// Read the file
debugSection("FORMAT", is);
skip(3, is); // type: ensight gold
debugSection("GEOMETRY", is);
readSkip(is, 2, meshFileName_);
debugSection("VARIABLE", is);
DynamicList<word> fieldNames(10);
DynamicList<word> fieldFileNames(10);
word fieldName;
word fieldFileName;
while (is.good())
{
word primitiveType(is); // scalar, vector
if (primitiveType == "TIME")
{
break;
}
readSkip(is, 3, fieldName); // p, U etc
fieldNames.append(fieldName);
is >> fieldFileName; // surfaceName.****.fieldName
fieldFileNames.append(fieldFileName);
}
fieldNames_.transfer(fieldNames);
fieldFileNames_.transfer(fieldFileNames);
if (debug)
{
Info<< "fieldNames: " << fieldNames_ << nl
<< "fieldFileNames: " << fieldFileNames_ << endl;
}
// Start reading time information
skip(3, is); // time set: 1
readSkip(is, 3, nTimeSteps_);
readSkip(is, 3, timeStartIndex_);
readSkip(is, 2, timeIncrement_);
if (debug)
{
Info<< "nTimeSteps: " << nTimeSteps_ << nl
<< "timeStartIndex: " << timeStartIndex_ << nl
<< "timeIncrement: " << timeIncrement_ << endl;
}
// Read the time values
skip(2, is);
timeValues_.setSize(nTimeSteps_);
for (label i = 0; i < nTimeSteps_; i++)
{
scalar t(readScalar(is));
timeValues_[i].value() = t;
// TODO: use character representation of t directly instead of
// regenerating from scalar value
timeValues_[i].name() = Foam::name(t);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ensightSurfaceReader::ensightSurfaceReader(const fileName& fName)
:
surfaceReader(fName),
streamFormat_(IOstream::ASCII),
baseDir_(fName.path()),
meshFileName_(),
fieldNames_(),
fieldFileNames_(),
nTimeSteps_(0),
timeStartIndex_(0),
timeIncrement_(1),
timeValues_(),
surfPtr_(NULL)
{
IFstream is(fName);
readCase(is);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::ensightSurfaceReader::~ensightSurfaceReader()
{}
// * * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * //
const Foam::meshedSurface& Foam::ensightSurfaceReader::geometry()
{
if (debug)
{
InfoInFunction<< endl;
}
if (!surfPtr_.valid())
{
IFstream isBinary(baseDir_/meshFileName_, IOstream::BINARY);
if (!isBinary.good())
{
FatalErrorInFunction
<< "Cannot read file " << isBinary.name()
<< exit(FatalError);
}
streamFormat_ = IOstream::BINARY;
{
istream& is = isBinary.stdStream();
char buffer[80];
is.read(buffer, 80);
char test[80];
label nChar = 0;
for (label i = 0; i < 80; ++i)
{
if (buffer[i] == '\0')
{
break;
}
test[i] = buffer[i];
nChar++;
}
string testStr(test, nChar);
if
(
(testStr.find("binary", 0) == string::npos)
&& (testStr.find("Binary", 0) == string::npos)
)
{
streamFormat_ = IOstream::ASCII;
}
}
if (debug)
{
Info<< "stream format: ";
if (streamFormat_ == IOstream::ASCII)
{
Info<< "ascii" << endl;
}
else
{
Info<< "binary" << endl;
}
}
ensightReadFile is(baseDir_/meshFileName_, streamFormat_);
if (debug)
{
Info<< "File: " << is.name() << endl;
}
readGeometryHeader(is);
label nPoints;
is.read(nPoints);
if (debug)
{
Info<< "nPoints: " << nPoints << endl;
}
pointField points(nPoints);
{
scalarField x(nPoints);
for (label dir = 0; dir < 3; dir++)
{
forAll(points, pointI)
{
is.read(x[pointI]);
}
points.replace(dir, x);
}
}
// Read faces - may be a mix of tris, quads and polys
DynamicList<face> faces(ceil(nPoints/3));
DynamicList<Tuple2<string, label> > schema(faces.size());
string faceType = "";
label nFace = 0;
while (is.good()) // (is.peek() != EOF)
{
is.read(faceType);
if (!is.good())
{
break;
}
if (debug)
{
Info<< "faceType: " << faceType << endl;
}
if (faceType == "tria3")
{
is.read(nFace);
label np = 3;
for (label faceI = 0; faceI < nFace; ++faceI)
{
face f(np);
for (label fpI = 0; fpI < np; fpI++)
{
is.read(f[fpI]);
}
faces.append(f);
}
}
else if (faceType == "quad4")
{
is.read(nFace);
label np = 4;
for (label faceI = 0; faceI < nFace; ++faceI)
{
face f(np);
for (label fpI = 0; fpI < np; fpI++)
{
is.read(f[fpI]);
}
faces.append(f);
}
}
else if (faceType == "nsided")
{
is.read(nFace);
labelList np(nFace);
for (label faceI = 0; faceI < nFace; ++faceI)
{
is.read(np[faceI]);
}
for (label faceI = 0; faceI < nFace; ++faceI)
{
face f(np[faceI]);
for (label fpI = 0; fpI < f.size(); ++fpI)
{
is.read(f[fpI]);
}
faces.append(f);
}
}
else
{
if (debug)
{
WarningInFunction
<< "Unknown face type: " << faceType
<< ". Aborting read and continuing with current "
<< "elements only" << endl;
}
break;
}
schema.append(Tuple2<string, label>(faceType, nFace));
}
schema_.transfer(schema);
if (debug)
{
Info<< "read nFaces: " << faces.size() << nl
<< "file schema: " << schema_ << endl;
}
// Convert from 1-based Ensight addressing to 0-based OF addressing
forAll(faces, faceI)
{
face& f = faces[faceI];
forAll(f, fpI)
{
f[fpI]--;
}
}
surfPtr_.reset(new meshedSurface(xferMove(points), faces.xfer()));
}
return surfPtr_();
}
Foam::instantList Foam::ensightSurfaceReader::times() const
{
return timeValues_;
}
Foam::wordList Foam::ensightSurfaceReader::fieldNames
(
const label timeIndex
) const
{
return fieldNames_;
}
Foam::tmp<Foam::Field<Foam::scalar> > Foam::ensightSurfaceReader::field
(
const label timeIndex,
const label fieldIndex,
const scalar& refValue
) const
{
return readField<scalar>(timeIndex, fieldIndex);
}
Foam::tmp<Foam::Field<Foam::vector> > Foam::ensightSurfaceReader::field
(
const label timeIndex,
const label fieldIndex,
const vector& refValue
) const
{
return readField<vector>(timeIndex, fieldIndex);
}
Foam::tmp<Foam::Field<Foam::sphericalTensor> >
Foam::ensightSurfaceReader::field
(
const label timeIndex,
const label fieldIndex,
const sphericalTensor& refValue
) const
{
return readField<sphericalTensor>(timeIndex, fieldIndex);
}
Foam::tmp<Foam::Field<Foam::symmTensor> > Foam::ensightSurfaceReader::field
(
const label timeIndex,
const label fieldIndex,
const symmTensor& refValue
) const
{
return readField<symmTensor>(timeIndex, fieldIndex);
}
Foam::tmp<Foam::Field<Foam::tensor> > Foam::ensightSurfaceReader::field
(
const label timeIndex,
const label fieldIndex,
const tensor& refValue
) const
{
return readField<tensor>(timeIndex, fieldIndex);
}
// ************************************************************************* //

201
src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H Normal file
View File

@ -0,0 +1,201 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::ensightensightSurfaceReader
Description
Ensight format surface reader
SourceFiles
ensightSurfaceReader.C
\*---------------------------------------------------------------------------*/
#ifndef ensightSurfaceReader_H
#define ensightSurfaceReader_H
#include "surfaceReader.H"
#include "ensightReadFile.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ensightSurfaceReader Declaration
\*---------------------------------------------------------------------------*/
class ensightSurfaceReader
:
public surfaceReader
{
protected:
// Protected Data
//- Format flag
IOstream::streamFormat streamFormat_;
//- Base directory
fileName baseDir_;
//- Name of mesh file
word meshFileName_;
//- Field names
List<word> fieldNames_;
//- Field file names
List<word> fieldFileNames_;
//- Number of time steps
label nTimeSteps_;
//- Start time index
label timeStartIndex_;
//- Time increment
label timeIncrement_;
//- Times
instantList timeValues_;
//- Pointer to the surface
autoPtr<meshedSurface> surfPtr_;
List<Tuple2<string, label> > schema_;
// Protected Member Functions
//- Read and check a section header
void debugSection(const word& expected, IFstream& is) const;
//- Read the case file
void readCase(IFstream& is);
//- Helper function to skip forward n steps in stream
void skip(const label n, IFstream& is) const;
//- Read (and throw away) geometry file header
void readGeometryHeader(ensightReadFile& is) const;
//- Helper function to return Type after skipping n tokens
template<class Type>
void readSkip(IFstream& is, const label nSkip, Type& value) const;
//- Helper function to return a field
template<class Type>
tmp<Field<Type> > readField
(
const label timeIndex,
const label fieldIndex
) const;
public:
//- Runtime type information
TypeName("ensight");
// Constructors
//- Construct from fileName
ensightSurfaceReader(const fileName& fName);
//- Destructor
virtual ~ensightSurfaceReader();
// Member Functions
//- Return a reference to the surface geometry
virtual const meshedSurface& geometry();
//- Return a list of the available times
virtual instantList times() const;
//- Return a list of the available fields at a given time
virtual wordList fieldNames(const label timeIndex) const;
//- Return a scalar field at a given time
virtual tmp<Field<scalar> > field
(
const label timeIndex,
const label fieldIndex,
const scalar& refValue = pTraits<scalar>::zero
) const;
//- Return a scalar field at a given time
virtual tmp<Field<vector> > field
(
const label timeIndex,
const label fieldIndex,
const vector& refValue = pTraits<vector>::zero
) const;
//- Return a sphericalTensor field at a given time
virtual tmp<Field<sphericalTensor> > field
(
const label timeIndex,
const label fieldIndex,
const sphericalTensor& reValue = pTraits<sphericalTensor>::zero
) const;
//- Return a symmTensor field at a given time
virtual tmp<Field<symmTensor> > field
(
const label timeIndex,
const label fieldIndex,
const symmTensor& reValue = pTraits<symmTensor>::zero
) const;
//- Return a tensor field at a given time
virtual tmp<Field<tensor> > field
(
const label timeIndex,
const label fieldIndex,
const tensor& reValue = pTraits<tensor>::zero
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "ensightSurfaceReaderTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

161
src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C Normal file
View File

@ -0,0 +1,161 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include <iomanip>
#include <sstream>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
void Foam::ensightSurfaceReader::readSkip
(
IFstream& is,
const label nSkip,
Type& value
) const
{
skip(nSkip, is);
is >> value;
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
(
const label timeIndex,
const label fieldIndex
) const
{
if (debug)
{
InfoInFunction<< endl;
}
const word& fieldName(fieldNames_[fieldIndex]);
const label fileIndex = timeStartIndex_ + timeIndex*timeIncrement_;
fileName fieldFileName(fieldFileNames_[fieldIndex]);
std::ostringstream oss;
oss << std::setfill('0') << std::setw(4) << fileIndex;
const word indexStr = oss.str();
fieldFileName.replace("****", indexStr);
ensightReadFile is(baseDir_/fieldFileName, streamFormat_);
if (!is.good())
{
FatalErrorInFunction
<< "Cannot read file " << is.name()
<< " for field " << fieldName
<< exit(FatalError);
}
// Check that data type is as expected
string primitiveType;
is.read(primitiveType);
if (debug)
{
Info<< "primitiveType: " << primitiveType << endl;
}
if (primitiveType != pTraits<Type>::typeName)
{
FatalIOErrorInFunction(is)
<< "Expected " << pTraits<Type>::typeName << "values "
<< "but found type " << primitiveType
<< exit(FatalIOError);
}
scalar value;
string strValue;
label iValue;
// Read header info: part index, e.g. part 1
is.read(strValue);
is.read(iValue);
// Allocate storage for data as a list per component
List<DynamicList<scalar> > values(pTraits<Type>::nComponents);
label n = surfPtr_->size();
forAll(values, cmptI)
{
values.setSize(n);
}
// Read data file using schema generated while reading the surface
forAll(schema_, i)
{
if (debug)
{
const string& faceType = schema_[i].first();
Info<< "Reading face type " << faceType << " data" << endl;
}
const label nFace = schema_[i].second();
if (nFace != 0)
{
is.read(strValue);
for
(
direction cmptI=0;
cmptI < pTraits<Type>::nComponents;
++cmptI
)
{
for (label faceI = 0; faceI < nFace; ++faceI)
{
is.read(value);
values[cmptI].append(value);
}
}
}
}
tmp<Field<Type> > tField(new Field<Type>(n, pTraits<Type>::zero));
Field<Type>& field = tField.ref();
for
(
direction cmptI=0;
cmptI < pTraits<Type>::nComponents;
++cmptI
)
{
field.replace(cmptI, values[cmptI]);
values[cmptI].clear();
}
return tField;
}
// ************************************************************************* //

51
src/sampling/sampledSurface/readers/surfaceReader.C Normal file
View File

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "surfaceReader.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(surfaceReader, 0);
defineRunTimeSelectionTable(surfaceReader, fileName);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::surfaceReader::surfaceReader(const fileName& fName)
:
fileName_(fName)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::surfaceReader::~surfaceReader()
{}
// ************************************************************************* //

161
src/sampling/sampledSurface/readers/surfaceReader.H Normal file
View File

@ -0,0 +1,161 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::surfaceReader
Description
Base class for surface readers
SourceFiles
surfaceReader.C
\*---------------------------------------------------------------------------*/
#ifndef surfaceReader_H
#define surfaceReader_H
#include "typeInfo.H"
#include "autoPtr.H"
#include "MeshedSurfaces.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class surfaceReader Declaration
\*---------------------------------------------------------------------------*/
class surfaceReader
{
protected:
//- File name
fileName fileName_;
public:
//- Runtime type information
TypeName("surfaceReader");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
surfaceReader,
fileName,
(
const fileName& fName
),
(fName)
);
// Selectors
//- Return a reference to the selected surfaceReader
static autoPtr<surfaceReader> New
(
const word& readType,
const fileName& fName
);
// Constructors
//- Construct from fileName
surfaceReader(const fileName& fName);
//- Destructor
virtual ~surfaceReader();
// Member Functions
//- Return a reference to the surface geometry
virtual const meshedSurface& geometry() = 0;
//- Return a list of the available times
virtual instantList times() const = 0;
//- Return a list of the available fields at a given time
virtual wordList fieldNames(const label timeIndex) const = 0;
//- Return a scalar field at a given time
virtual tmp<Field<scalar> > field
(
const label timeIndex,
const label fieldIndex,
const scalar& refValue = pTraits<scalar>::zero
) const = 0;
//- Return a vector field at a given time
virtual tmp<Field<vector> > field
(
const label timeIndex,
const label fieldIndex,
const vector& refValue = pTraits<vector>::zero
) const = 0;
//- Return a sphericalTensor field at a given time
virtual tmp<Field<sphericalTensor> > field
(
const label timeIndex,
const label fieldIndex,
const sphericalTensor& reValue = pTraits<sphericalTensor>::zero
) const = 0;
//- Return a symmTensor field at a given time
virtual tmp<Field<symmTensor> > field
(
const label timeIndex,
const label fieldIndex,
const symmTensor& reValue = pTraits<symmTensor>::zero
) const = 0;
//- Return a tensor field at a given time
virtual tmp<Field<tensor> > field
(
const label timeIndex,
const label fieldIndex,
const tensor& reValue = pTraits<tensor>::zero
) const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

52
src/sampling/sampledSurface/readers/surfaceReaderNew.C Normal file
View File

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "surfaceReader.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
Foam::autoPtr<Foam::surfaceReader> Foam::surfaceReader::New
(
const word& readerType,
const fileName& fName
)
{
fileNameConstructorTable::iterator cstrIter =
fileNameConstructorTablePtr_->find(readerType);
if (cstrIter == fileNameConstructorTablePtr_->end())
{
FatalErrorInFunction
<< "Unknown reader type \"" << readerType << "\"\n\n"
<< "Valid reader types: "
<< fileNameConstructorTablePtr_->sortedToc() << nl
<< exit(FatalError);
}
return autoPtr<surfaceReader>(cstrIter()(fName));
}
// ************************************************************************* //