Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2012-04-18 16:34:07 +01:00
26 changed files with 2424 additions and 123 deletions

View File

@ -0,0 +1,206 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
Application
testExtendedStencil2
Description
Test app for determining extended cell-to-cell stencil.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "Time.H"
#include "OFstream.H"
#include "meshTools.H"
#include "CFCCellToCellStencil.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void writeStencilOBJ
(
const fileName& fName,
const point& fc,
const List<point>& stencilCc
)
{
OFstream str(fName);
label vertI = 0;
meshTools::writeOBJ(str, fc);
vertI++;
forAll(stencilCc, i)
{
meshTools::writeOBJ(str, stencilCc[i]);
vertI++;
str << "l 1 " << vertI << nl;
}
}
// Stats
void writeStencilStats(const labelListList& stencil)
{
label sumSize = 0;
label nSum = 0;
label minSize = labelMax;
label maxSize = labelMin;
forAll(stencil, i)
{
const labelList& sCells = stencil[i];
if (sCells.size() > 0)
{
sumSize += sCells.size();
nSum++;
minSize = min(minSize, sCells.size());
maxSize = max(maxSize, sCells.size());
}
}
reduce(sumSize, sumOp<label>());
reduce(nSum, sumOp<label>());
sumSize /= nSum;
reduce(minSize, minOp<label>());
reduce(maxSize, maxOp<label>());
Info<< "Stencil size :" << nl
<< " average : " << sumSize << nl
<< " min : " << minSize << nl
<< " max : " << maxSize << nl
<< endl;
}
// Main program:
int main(int argc, char *argv[])
{
# include "addTimeOptions.H"
# include "setRootCase.H"
# include "createTime.H"
// Get times list
instantList Times = runTime.times();
# include "checkTimeOptions.H"
runTime.setTime(Times[startTime], startTime);
# include "createMesh.H"
//---- CELL CENTRED STENCIL -----
// Centred, cell stencil
// ~~~~~~~~~~~~~~~~~~~~~
{
// Full stencil. This is per local cell a set of global indices, either
// into cells or also boundary faces.
CFCCellToCellStencil stencil(mesh);
// Construct exchange map and renumber stencil
List<Map<label> > compactMap(Pstream::nProcs());
mapDistribute map
(
stencil.globalNumbering(),
stencil,
compactMap
);
// Print some stats
Info<< "cellFaceCell:" << endl;
writeStencilStats(stencil);
// Collect the data in stencil form
List<List<vector> > stencilPoints;
{
const volVectorField& fld = mesh.C();
// 1. Construct cell data in compact addressing
List<point> compactFld(map.constructSize(), pTraits<point>::zero);
// Insert my internal values
forAll(fld, cellI)
{
compactFld[cellI] = fld[cellI];
}
// Insert my boundary values
label nCompact = fld.size();
forAll(fld.boundaryField(), patchI)
{
const fvPatchField<vector>& pfld = fld.boundaryField()[patchI];
forAll(pfld, i)
{
compactFld[nCompact++] = pfld[i];
}
}
// Do all swapping
map.distribute(compactFld);
// 2. Pull to stencil
stencilPoints.setSize(stencil.size());
forAll(stencil, cellI)
{
const labelList& compactCells = stencil[cellI];
stencilPoints[cellI].setSize(compactCells.size());
forAll(compactCells, i)
{
stencilPoints[cellI][i] = compactFld[compactCells[i]];
}
}
}
forAll(stencilPoints, cellI)
{
writeStencilOBJ
(
runTime.path()/"centredCell" + Foam::name(cellI) + ".obj",
mesh.cellCentres()[cellI],
stencilPoints[cellI]
);
}
}
Pout<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,6 +37,8 @@ Description
#include "surfaceFields.H"
#include "pointFields.H"
#include "ReadFields.H"
#include "interpolationWeights.H"
#include "uniformInterpolate.H"
using namespace Foam;
@ -48,6 +50,8 @@ class fieldInterpolator
const HashSet<word>& selectedFields_;
instant ti_;
instant ti1_;
const interpolationWeights& interpolator_;
const wordList& timeNames_;
int divisions_;
public:
@ -60,6 +64,8 @@ public:
const HashSet<word>& selectedFields,
const instant& ti,
const instant& ti1,
const interpolationWeights& interpolator,
const wordList& timeNames,
int divisions
)
:
@ -69,6 +75,8 @@ public:
selectedFields_(selectedFields),
ti_(ti),
ti1_(ti1),
interpolator_(interpolator),
timeNames_(timeNames),
divisions_(divisions)
{}
@ -98,34 +106,6 @@ void fieldInterpolator::interpolate()
{
Info<< " " << fieldIter()->name() << '(';
GeoFieldType fieldi
(
IOobject
(
fieldIter()->name(),
ti_.name(),
fieldIter()->db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh_
);
GeoFieldType fieldi1
(
IOobject
(
fieldIter()->name(),
ti1_.name(),
fieldIter()->db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh_
);
scalar deltaT = (ti1_.value() - ti_.value())/(divisions_ + 1);
for (int j=0; j<divisions_; j++)
@ -141,20 +121,51 @@ void fieldInterpolator::interpolate()
Info<< " ";
}
scalar lambda = scalar(j + 1)/scalar(divisions_ + 1);
// Calculate times to read and weights
labelList indices;
scalarField weights;
interpolator_.valueWeights
(
runTime_.value(),
indices,
weights
);
const wordList selectedTimeNames
(
UIndirectList<word>(timeNames_, indices)()
);
//Info<< "For time " << runTime_.value()
// << " need times " << selectedTimeNames
// << " need weights " << weights << endl;
// Read on the objectRegistry all the required fields
ReadFields<GeoFieldType>
(
fieldIter()->name(),
mesh_,
selectedTimeNames
);
GeoFieldType fieldj
(
IOobject
uniformInterpolate<GeoFieldType>
(
IOobject
(
fieldIter()->name(),
runTime_.timeName(),
fieldIter()->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
fieldIter()->name(),
timej.name(),
fieldIter()->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
(1.0 - lambda)*fieldi + lambda*fieldi1
selectedTimeNames,
weights
)
);
fieldj.write();
@ -188,6 +199,12 @@ int main(int argc, char *argv[])
"integer",
"specify number of temporal sub-divisions to create (default = 1)."
);
argList::addOption
(
"interpolationType",
"word",
"specify type of interpolation (linear or spline)"
);
#include "setRootCase.H"
#include "createTime.H"
@ -198,15 +215,51 @@ int main(int argc, char *argv[])
{
args.optionLookup("fields")() >> selectedFields;
}
if (selectedFields.size())
{
Info<< "Interpolating fields " << selectedFields << nl << endl;
}
else
{
Info<< "Interpolating all fields" << nl << endl;
}
int divisions = 1;
if (args.optionFound("divisions"))
{
args.optionLookup("divisions")() >> divisions;
}
Info<< "Using " << divisions << " per time interval" << nl << endl;
const word interpolationType = args.optionLookupOrDefault<word>
(
"interpolationType",
"linear"
);
Info<< "Using interpolation " << interpolationType << nl << endl;
instantList timeDirs = timeSelector::select0(runTime, args);
scalarField timeVals(timeDirs.size());
wordList timeNames(timeDirs.size());
forAll(timeDirs, i)
{
timeVals[i] = timeDirs[i].value();
timeNames[i] = timeDirs[i].name();
}
autoPtr<interpolationWeights> interpolatorPtr
(
interpolationWeights::New
(
interpolationType,
timeVals
)
);
#include "createMesh.H"
Info<< "Interpolating fields for times:" << endl;
@ -226,6 +279,8 @@ int main(int argc, char *argv[])
selectedFields,
timeDirs[timei],
timeDirs[timei+1],
interpolatorPtr(),
timeNames,
divisions
);

View File

@ -576,6 +576,13 @@ $(interpolations)/interpolationTable/tableReaders/tableReaders.C
$(interpolations)/interpolationTable/tableReaders/openFoam/openFoamTableReaders.C
$(interpolations)/interpolationTable/tableReaders/csv/csvTableReaders.C
interpolationWeights = $(interpolations)/interpolationWeights
$(interpolationWeights)/interpolationWeights/interpolationWeights.C
$(interpolationWeights)/linearInterpolationWeights/linearInterpolationWeights.C
$(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
algorithms/indexedOctree/indexedOctreeName.C
algorithms/indexedOctree/treeDataCell.C

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -119,6 +119,12 @@ public:
inline void operator=(const T&);
// STL type definitions
//- Type of values the UList contains.
typedef T value_type;
// Ostream operator
//- Write UIndirectList to Ostream

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -151,9 +151,25 @@ Foam::wordList Foam::objectRegistry::sortedNames(const word& ClassName) const
const Foam::objectRegistry& Foam::objectRegistry::subRegistry
(
const word& name
const word& name,
const bool forceCreate
) const
{
if (forceCreate && !foundObject<objectRegistry>(name))
{
objectRegistry* fieldsCachePtr = new objectRegistry
(
IOobject
(
name,
time().constant(),
*this,
IOobject::NO_READ,
IOobject::NO_WRITE
)
);
fieldsCachePtr->store();
}
return lookupObject<objectRegistry>(name);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -147,8 +147,13 @@ public:
template<class Type>
wordList names() const;
//- Lookup and return a const sub-objectRegistry
const objectRegistry& subRegistry(const word& name) const;
//- Lookup and return a const sub-objectRegistry. Optionally create
// it if it does not exist.
const objectRegistry& subRegistry
(
const word& name,
const bool forceCreate = false
) const;
//- Lookup and return all objects of the given Type
template<class Type>

View File

@ -0,0 +1,142 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class GeoField>
Foam::tmp<GeoField> Foam::uniformInterpolate
(
const HashPtrTable<GeoField, label, Hash<label> >& fields,
const labelList& indices,
const scalarField& weights
)
{
const GeoField& field0 = *(*fields.begin());
// Interpolate
tmp<GeoField> tfld
(
new GeoField
(
IOobject
(
"uniformInterpolate(" + field0.name() + ')',
field0.time().timeName(),
field0.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
weights[0]*(*fields[indices[0]])
)
);
GeoField& fld = tfld();
for (label i = 1; i < indices.size(); ++i)
{
fld += weights[i]*(*fields[indices[i]]);
}
return tfld;
}
template<class GeoField>
Foam::tmp<GeoField> Foam::uniformInterpolate
(
const IOobject& fieldIO,
const word& fieldName,
const wordList& times,
const scalarField& weights,
const objectRegistry& fieldsCache
)
{
// Look up the first field
const objectRegistry& time0Fields = fieldsCache.lookupObject
<
const objectRegistry
>
(
times[0]
);
const GeoField& field0 = time0Fields.lookupObject
<
const GeoField
>
(
fieldName
);
// Interpolate
tmp<GeoField> tfld(new GeoField(fieldIO, weights[0]*field0));
GeoField& fld = tfld();
for (label i = 1; i < times.size(); ++i)
{
const objectRegistry& timeIFields = fieldsCache.lookupObject
<
const objectRegistry
>
(
times[i]
);
const GeoField& fieldI = timeIFields.lookupObject
<
const GeoField
>
(
fieldName
);
fld += weights[i]*fieldI;
}
return tfld;
}
template<class GeoField>
Foam::tmp<GeoField> Foam::uniformInterpolate
(
const IOobject& fieldIO,
const word& fieldName,
const wordList& times,
const scalarField& weights,
const word& registryName
)
{
return uniformInterpolate<GeoField>
(
fieldIO,
fieldName,
times,
weights,
fieldIO.db().subRegistry(registryName, true)
);
}
// ************************************************************************* //

View File

@ -0,0 +1,83 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "GeometricField.H"
#include "HashPtrTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Global functions * * * * * * * * * * * * * //
//- Interpolate selected fields (given by indices and corresponding weights)
// (boundary type becomes calculated). Fields stored per index. Field gets name
// "uniformInterpolate(" + fld.name() + ')'.
template<class GeoField>
tmp<GeoField> uniformInterpolate
(
const HashPtrTable<GeoField, label, Hash<label> >& fields,
const labelList& indices,
const scalarField& weights
);
//- Interpolate fields. fieldsCache contains per timeName all loaded fields.
// Resulting field gets properties according to fieldIO
template<class GeoField>
tmp<GeoField> uniformInterpolate
(
const IOobject& fieldIO,
const word& fieldName,
const wordList& times,
const scalarField& weights,
const objectRegistry& fieldsCache
);
//- Interpolate fields. fieldsCache contains per timeName all loaded fields.
// Resulting field gets properties according to fieldIO
template<class GeoField>
tmp<GeoField> uniformInterpolate
(
const IOobject& fieldIO,
const word& fieldName,
const wordList& times,
const scalarField& weights,
const word& registryName = "fieldsCache"
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "uniformInterpolate.C"
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -127,4 +127,134 @@ Foam::wordList Foam::ReadFields
}
template<class GeoField>
void Foam::ReadFields
(
const word& fieldName,
const typename GeoField::Mesh& mesh,
const wordList& timeNames,
objectRegistry& fieldsCache
)
{
// Collect all times that are no longer used
{
HashSet<word> usedTimes(timeNames);
DynamicList<word> unusedTimes(fieldsCache.size());
forAllIter(objectRegistry, fieldsCache, timeIter)
{
const word& tm = timeIter.key();
if (!usedTimes.found(tm))
{
unusedTimes.append(tm);
}
}
//Info<< "Unloading times " << unusedTimes << endl;
forAll(unusedTimes, i)
{
objectRegistry& timeCache = const_cast<objectRegistry&>
(
fieldsCache.lookupObject<objectRegistry>(unusedTimes[i])
);
fieldsCache.checkOut(timeCache);
}
}
// Load any new fields
forAll(timeNames, i)
{
const word& tm = timeNames[i];
// Create if not found
if (!fieldsCache.found(tm))
{
//Info<< "Creating registry for time " << tm << endl;
// Create objectRegistry if not found
objectRegistry* timeCachePtr = new objectRegistry
(
IOobject
(
tm,
tm,
fieldsCache,
IOobject::NO_READ,
IOobject::NO_WRITE
)
);
timeCachePtr->store();
}
// Obtain cache for current time
const objectRegistry& timeCache =
fieldsCache.lookupObject<objectRegistry>
(
tm
);
// Store field if not found
if (!timeCache.found(fieldName))
{
//Info<< "Loading field " << fieldName
// << " for time " << tm << endl;
GeoField loadedFld
(
IOobject
(
fieldName,
tm,
mesh.thisDb(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh
);
// Transfer to timeCache (new objectRegistry and store flag)
GeoField* fldPtr = new GeoField
(
IOobject
(
fieldName,
tm,
timeCache,
IOobject::NO_READ,
IOobject::NO_WRITE
),
loadedFld
);
fldPtr->store();
}
}
}
template<class GeoField>
void Foam::ReadFields
(
const word& fieldName,
const typename GeoField::Mesh& mesh,
const wordList& timeNames,
const word& registryName
)
{
ReadFields<GeoField>
(
fieldName,
mesh,
timeNames,
const_cast<objectRegistry&>
(
mesh.thisDb().subRegistry(registryName, true)
)
);
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -55,6 +55,28 @@ wordList ReadFields
const bool syncPar = true
);
//- Helper routine to read GeometricFields. The fieldsCache is per time
// an objectRegistry of all stored fields
template<class GeoField>
static void ReadFields
(
const word& fieldName,
const typename GeoField::Mesh& mesh,
const wordList& timeNames,
objectRegistry& fieldsCache
);
//- Helper routine to read GeometricFields. The fieldsCache is per time
// an objectRegistry of all stored fields
template<class GeoField>
static void ReadFields
(
const word& fieldName,
const typename GeoField::Mesh& mesh,
const wordList& timeNames,
const word& registryName = "fieldsCache"
);
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,121 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "interpolationWeights.H"
#include "addToRunTimeSelectionTable.H"
#include "Time.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defineTypeNameAndDebug(interpolationWeights, 0);
defineRunTimeSelectionTable(interpolationWeights, word);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
interpolationWeights::interpolationWeights
(
const scalarField& samples
)
:
samples_(samples)
{}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
autoPtr<interpolationWeights> interpolationWeights::New
(
const word& type,
const scalarField& samples
)
{
Info<< nl << "Selecting interpolationWeights "
<< type << endl;
wordConstructorTable::iterator cstrIter =
wordConstructorTablePtr_->find(type);
if (cstrIter == wordConstructorTablePtr_->end())
{
FatalErrorIn
(
"interpolationWeights::New(const word&, "
"const scalarField&)"
) << "Unknown interpolationWeights type "
<< type
<< endl << endl
<< "Valid interpolationWeights types are :" << endl
<< wordConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<interpolationWeights>(cstrIter()(samples));
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
interpolationWeights::~interpolationWeights()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//objectRegistry& interpolationWeights::registry
//(
// const objectRegistry& obr,
// const word& name
//)
//{
// if (!obr.foundObject<objectRegistry>(name))
// {
// objectRegistry* fieldsCachePtr = new objectRegistry
// (
// IOobject
// (
// name,
// obr.time().constant(),
// obr,
// IOobject::NO_READ,
// IOobject::NO_WRITE
// )
// );
// fieldsCachePtr->store();
// }
// return const_cast<objectRegistry&>(obr.subRegistry(name));
//}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,159 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
Class
Foam::interpolationWeights
Description
Abstract base class for interpolating in 1D
SourceFiles
interpolationWeights.C
interpolationWeightsTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef interpolationWeights_H
#define interpolationWeights_H
#include "scalarField.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "pointField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class fvMesh;
class objectRegistry;
/*---------------------------------------------------------------------------*\
Class interpolationWeights Declaration
\*---------------------------------------------------------------------------*/
class interpolationWeights
{
private:
// Private Member Functions
//- Disallow default bitwise copy construct
interpolationWeights(const interpolationWeights&);
//- Disallow default bitwise assignment
void operator=(const interpolationWeights&);
protected:
const scalarField& samples_;
public:
//- Runtime type information
TypeName("interpolationWeights");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
interpolationWeights,
word,
(
const scalarField& samples
),
(samples)
);
// Constructors
//- Construct from components
interpolationWeights(const scalarField& samples);
// Selectors
//- Return a reference to the selected interpolationWeights
static autoPtr<interpolationWeights> New
(
const word& type,
const scalarField& samples
);
//- Destructor
virtual ~interpolationWeights();
// Member Functions
//- Calculate weights and indices to calculate t from samples.
// Returns true if indices changed.
virtual bool valueWeights
(
const scalar t,
labelList& indices,
scalarField& weights
) const = 0;
//- Calculate weights and indices to calculate integrand of t1..t2
// from samples. Returns true if indices changed.
virtual bool integrationWeights
(
const scalar t1,
const scalar t2,
labelList& indices,
scalarField& weights
) const = 0;
//- Helper: weighted sum
template<class ListType1, class ListType2>
static typename outerProduct
<
typename ListType1::value_type,
typename ListType2::value_type
>::type
weightedSum(const ListType1& f1, const ListType2& f2);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "interpolationWeightsTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,77 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "interpolationWeights.H"
#include "ListOps.H"
#include "IOobject.H"
#include "HashSet.H"
#include "objectRegistry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ListType1, class ListType2>
typename Foam::outerProduct
<
typename ListType1::value_type,
typename ListType2::value_type
>::type
Foam::interpolationWeights::weightedSum
(
const ListType1& f1,
const ListType2& f2
)
{
typedef typename outerProduct
<
typename ListType1::value_type,
typename ListType2::value_type
>::type returnType;
if (f1.size())
{
returnType SumProd = f1[0]*f2[0];
for (label i = 1; i < f1.size(); ++i)
{
SumProd += f1[i]*f2[i];
}
return SumProd;
}
else
{
return pTraits<returnType>::zero;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,249 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "linearInterpolationWeights.H"
#include "addToRunTimeSelectionTable.H"
#include "ListOps.H"
#include "Pair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(linearInterpolationWeights, 0);
addToRunTimeSelectionTable
(
interpolationWeights,
linearInterpolationWeights,
word
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::Pair<Foam::scalar> linearInterpolationWeights::integrationWeights
(
const label i,
const scalar t
) const
{
// t is in range samples_[i] .. samples_[i+1]
scalar s = (t-samples_[i])/(samples_[i+1]-samples_[i]);
if (s < -SMALL || s > 1+SMALL)
{
FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
<< "Value " << t << " outside range " << samples_[i]
<< " .. " << samples_[i+1]
<< exit(FatalError);
}
scalar d = samples_[i+1]-t;
return Pair<scalar>(d*0.5*(1-s), d*0.5*(1+s));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
linearInterpolationWeights::linearInterpolationWeights
(
const scalarField& samples
)
:
interpolationWeights(samples),
index_(-1)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool linearInterpolationWeights::valueWeights
(
const scalar t,
labelList& indices,
scalarField& weights
) const
{
bool indexChanged = false;
// Check if current timeIndex is still valid
if
(
index_ >= 0
&& index_ < samples_.size()
&& (
samples_[index_] <= t
&& (index_ == samples_.size()-1 || t <= samples_[index_+1])
)
)
{
// index_ still at correct slot
}
else
{
// search for correct index
index_ = findLower(samples_, t);
indexChanged = true;
}
if (index_ == -1)
{
// Use first element only
indices.setSize(1);
weights.setSize(1);
indices[0] = 0;
weights[0] = 1.0;
}
else if (index_ == samples_.size()-1)
{
// Use last element only
indices.setSize(1);
weights.setSize(1);
indices[0] = samples_.size()-1;
weights[0] = 1.0;
}
else
{
// Interpolate
indices.setSize(2);
weights.setSize(2);
indices[0] = index_;
indices[1] = index_+1;
scalar t0 = samples_[indices[0]];
scalar t1 = samples_[indices[1]];
scalar deltaT = t1-t0;
weights[0] = (t1-t)/deltaT;
weights[1] = 1.0-weights[0];
}
return indexChanged;
}
bool linearInterpolationWeights::integrationWeights
(
const scalar t1,
const scalar t2,
labelList& indices,
scalarField& weights
) const
{
if (t2 < t1-VSMALL)
{
FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
<< "Integration should be in positive direction."
<< " t1:" << t1 << " t2:" << t2
<< exit(FatalError);
}
// Currently no fancy logic on cached index
label i1 = findLower(samples_, t1);
label i2 = findLower(samples_, t2);
// For now just fail if any outside table
if (i1 == -1 || i2 == samples_.size()-1)
{
FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
<< "Integrating outside table " << samples_[0] << ".."
<< samples_.last() << " not implemented."
<< " t1:" << t1 << " t2:" << t2 << exit(FatalError);
}
label nIndices = i2-i1+2;
// Determine if indices already correct
bool anyChanged = false;
if (nIndices != indices.size())
{
anyChanged = true;
}
else
{
// Closer check
label index = i1;
forAll(indices, i)
{
if (indices[i] != index)
{
anyChanged = true;
break;
}
index++;
}
}
indices.setSize(nIndices);
weights.setSize(nIndices);
weights = 0.0;
// Sum from i1+1 to i2+1
for (label i = i1+1; i <= i2; i++)
{
scalar d = samples_[i+1]-samples_[i];
indices[i-i1] = i;
weights[i-i1] += 0.5*d;
indices[i+1-i1] = i+1;
weights[i+1-i1] += 0.5*d;
}
// Add from i1 to t1
{
Pair<scalar> i1Tot1 = integrationWeights(i1, t1);
indices[0] = i1;
weights[0] += i1Tot1.first();
indices[1] = i1+1;
weights[1] += i1Tot1.second();
}
// Subtract from t2 to i2+1
{
Pair<scalar> wghts = integrationWeights(i2, t2);
indices[i2-i1] = i2;
weights[i2-i1] += -wghts.first();
indices[i2-i1+1] = i2+1;
weights[i2-i1+1] += -wghts.second();
}
return anyChanged;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,120 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
Class
Foam::linearInterpolationWeights
Description
SourceFiles
linearInterpolationWeights.C
\*---------------------------------------------------------------------------*/
#ifndef linearInterpolationWeights_H
#define linearInterpolationWeights_H
#include "interpolationWeights.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class linearInterpolationWeights Declaration
\*---------------------------------------------------------------------------*/
class linearInterpolationWeights
:
public interpolationWeights
{
private:
// Private data
//- Cached index in samples from previous invocation
mutable label index_;
// Private Member Functions
//- Get weights of i and i+1 to calculate integration from t to
// samples_[i+1]
Pair<scalar> integrationWeights
(
const label i,
const scalar t
) const;
public:
//- Runtime type information
TypeName("linear");
// Constructors
//- Construct from components
linearInterpolationWeights
(
const scalarField& samples
);
//- Destructor
virtual ~linearInterpolationWeights()
{}
// Member Functions
//- Calculate weights and indices to calculate t from samples.
// Returns true if indices changed.
virtual bool valueWeights
(
const scalar t,
labelList& indices,
scalarField& weights
) const;
//- Calculate weights and indices to calculate integrand of t1..t2
// from samples. Returns true if indices changed.
virtual bool integrationWeights
(
const scalar t1,
const scalar t2,
labelList& indices,
scalarField& weights
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,228 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "splineInterpolationWeights.H"
#include "addToRunTimeSelectionTable.H"
#include "ListOps.H"
#include "linearInterpolationWeights.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(splineInterpolationWeights, 0);
addToRunTimeSelectionTable
(
interpolationWeights,
splineInterpolationWeights,
word
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
splineInterpolationWeights::splineInterpolationWeights
(
const scalarField& samples,
const bool checkEqualDistance
)
:
interpolationWeights(samples),
index_(-1)
{
if (checkEqualDistance && samples_.size() > 2)
{
const scalar interval = samples_[1]-samples[0];
for (label i = 2; i < samples_.size(); i++)
{
scalar d = samples_[i]-samples[i-1];
if (mag(d-interval) > SMALL)
{
WarningIn
(
"splineInterpolationWeights::splineInterpolationWeights"
"(const scalarField&)"
) << "Spline interpolation only valid for constant intervals."
<< nl
<< "Interval 0-1 : " << interval << nl
<< "Interval " << i-1 << '-' << i << " : "
<< d << endl;
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool splineInterpolationWeights::valueWeights
(
const scalar t,
labelList& indices,
scalarField& weights
) const
{
bool indexChanged = false;
// linear interpolation
if (samples_.size() <= 2)
{
return linearInterpolationWeights(samples_).valueWeights
(
t,
indices,
weights
);
}
// Check if current timeIndex is still valid
if
(
index_ >= 0
&& index_ < samples_.size()
&& (
samples_[index_] <= t
&& (index_ == samples_.size()-1 || t <= samples_[index_+1])
)
)
{
// index_ still at correct slot
}
else
{
// search for correct index
index_ = findLower(samples_, t);
indexChanged = true;
}
// Clamp if outside table
if (index_ == -1)
{
indices.setSize(1);
weights.setSize(1);
indices[0] = 0;
weights[0] = 1;
return indexChanged;
}
else if (index_ == samples_.size()-1)
{
indices.setSize(1);
weights.setSize(1);
indices[0] = samples_.size()-1;
weights[0] = 1;
return indexChanged;
}
label lo = index_;
label hi = index_+1;
// weighting
scalar mu = (t - samples_[lo])/(samples_[hi] - samples_[lo]);
scalar w0 = 0.5*(mu*(-1+mu*(2-mu))); // coeff of lo-1
scalar w1 = 0.5*(2+mu*(mu*(-5 + mu*(3)))); // coeff of lo
scalar w2 = 0.5*(mu*(1 + mu*(4 + mu*(-3)))); // coeff of hi
scalar w3 = 0.5*(mu*mu*(-1 + mu)); // coeff of hi+1
if (lo > 0)
{
if (hi < samples_.size()-1)
{
// Four points available
indices.setSize(4);
weights.setSize(4);
indices[0] = lo-1;
indices[1] = lo;
indices[2] = hi;
indices[3] = hi+1;
weights[0] = w0;
weights[1] = w1;
weights[2] = w2;
weights[3] = w3;
}
else
{
// No y3 available. Extrapolate: y3=3*y2-y1
indices.setSize(3);
weights.setSize(3);
indices[0] = lo-1;
indices[1] = lo;
indices[2] = hi;
weights[0] = w0;
weights[1] = w1 - w3;
weights[2] = w2 + 2*w3;
}
}
else
{
// No y0 available. Extrapolate: y0=2*y1-y2;
if (hi < samples_.size()-1)
{
indices.setSize(3);
weights.setSize(3);
indices[0] = lo;
indices[1] = hi;
indices[2] = hi+1;
weights[0] = w1 + 2*w0;
weights[1] = w2 - w0;
weights[2] = w3;
}
else
{
indices.setSize(2);
weights.setSize(2);
indices[0] = lo;
indices[1] = hi;
weights[0] = w1 + 2*w0 - w3;
weights[1] = w2 - w0 + 2*w3;
}
}
return indexChanged;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,121 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
Class
Foam::splineInterpolationWeights
Description
Catmull-Rom spline interpolation.
SourceFiles
splineInterpolationWeights.C
\*---------------------------------------------------------------------------*/
#ifndef splineInterpolationWeights_H
#define splineInterpolationWeights_H
#include "interpolationWeights.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class splineInterpolationWeights Declaration
\*---------------------------------------------------------------------------*/
class splineInterpolationWeights
:
public interpolationWeights
{
private:
// Private data
//- Cached index in samples from previous invocation
mutable label index_;
public:
//- Runtime type information
TypeName("spline");
// Constructors
//- Construct from components. By default make sure samples are
// equidistant.
splineInterpolationWeights
(
const scalarField& samples,
const bool checkEqualDistance = true
);
//- Destructor
virtual ~splineInterpolationWeights()
{}
// Member Functions
//- Calculate weights and indices to calculate t from samples.
// Returns true if indices changed.
virtual bool valueWeights
(
const scalar t,
labelList& indices,
scalarField& weights
) const;
//- Calculate weights and indices to calculate integrand of t1..t2
// from samples. Returns true if indices changed.
virtual bool integrationWeights
(
const scalar t1,
const scalar t2,
labelList& indices,
scalarField& weights
) const
{
notImplemented
(
"splineInterpolationWeights::integrationWeights(..)"
);
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -26,6 +26,29 @@ License
#include "TableBase.H"
#include "Time.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
const Foam::interpolationWeights& Foam::TableBase<Type>::interpolator() const
{
if (interpolatorPtr_.empty())
{
// Re-work table into linear list
tableSamples_.setSize(table_.size());
forAll(table_, i)
{
tableSamples_[i] = table_[i].first();
}
interpolatorPtr_ = interpolationWeights::New
(
interpolationScheme_,
tableSamples_
);
}
return interpolatorPtr_();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
@ -39,6 +62,10 @@ Foam::TableBase<Type>::TableBase(const word& name, const dictionary& dict)
dict.lookupOrDefault<word>("outOfBounds", "clamp")
)
),
interpolationScheme_
(
dict.lookupOrDefault<word>("interpolationScheme", "linear")
),
table_(),
dimensions_(dimless)
{}
@ -49,8 +76,11 @@ Foam::TableBase<Type>::TableBase(const TableBase<Type>& tbl)
:
name_(tbl.name_),
boundsHandling_(tbl.boundsHandling_),
interpolationScheme_(tbl.interpolationScheme_),
table_(tbl.table_),
dimensions_(tbl.dimensions_)
dimensions_(tbl.dimensions_),
tableSamples_(tbl.tableSamples_),
interpolatorPtr_(tbl.interpolatorPtr_)
{}
@ -307,6 +337,9 @@ void Foam::TableBase<Type>::convertTimeBase(const Time& t)
scalar value = table_[i].first();
table_[i].first() = t.userTimeToTime(value);
}
tableSamples_.clear();
interpolatorPtr_.clear();
}
@ -325,88 +358,104 @@ Type Foam::TableBase<Type>::value(const scalar x) const
return table_.last().second();
}
// Find i such that x(i) < xDash < x(i+1)
label i = 0;
while ((table_[i+1].first() < xDash) && (i+1 < table_.size()))
// Use interpolator
interpolator().valueWeights(x, currentIndices_, currentWeights_);
Type t = currentWeights_[0]*table_[currentIndices_[0]].second();
for (label i = 1; i < currentIndices_.size(); i++)
{
i++;
t += currentWeights_[i]*table_[currentIndices_[i]].second();
}
return t;
Info <<
(xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
* (table_[i+1].second() - table_[i].second())
+ table_[i].second() << endl;
// Linear interpolation to find value
return Type
(
(xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
* (table_[i+1].second() - table_[i].second())
+ table_[i].second()
);
//// Find i such that x(i) < xDash < x(i+1)
//label i = 0;
//while ((table_[i+1].first() < xDash) && (i+1 < table_.size()))
//{
// i++;
//}
//
//// Linear interpolation to find value
//return Type
//(
// (xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
// * (table_[i+1].second() - table_[i].second())
// + table_[i].second()
//);
}
template<class Type>
Type Foam::TableBase<Type>::integrate(const scalar x1, const scalar x2) const
{
// Initialise return value
Type sum = pTraits<Type>::zero;
// Use interpolator
interpolator().integrationWeights(x1, x2, currentIndices_, currentWeights_);
// Return zero if out of bounds
if ((x1 > table_.last().first()) || (x2 < table_[0].first()))
Type sum = currentWeights_[0]*table_[currentIndices_[0]].second();
for (label i = 1; i < currentIndices_.size(); i++)
{
return sum;
sum += currentWeights_[i]*table_[currentIndices_[i]].second();
}
// Find next index greater than x1
label id1 = 0;
while ((table_[id1].first() < x1) && (id1 < table_.size()))
{
id1++;
}
// Find next index less than x2
label id2 = table_.size() - 1;
while ((table_[id2].first() > x2) && (id2 >= 1))
{
id2--;
}
if ((id1 - id2) == 1)
{
// x1 and x2 lie within 1 interval
sum = 0.5*(value(x1) + value(x2))*(x2 - x1);
}
else
{
// x1 and x2 cross multiple intervals
// Integrate table body
for (label i=id1; i<id2; i++)
{
sum +=
(table_[i].second() + table_[i+1].second())
* (table_[i+1].first() - table_[i].first());
}
sum *= 0.5;
// Add table ends (partial segments)
if (id1 > 0)
{
sum += 0.5
* (value(x1) + table_[id1].second())
* (table_[id1].first() - x1);
}
if (id2 < table_.size() - 1)
{
sum += 0.5
* (table_[id2].second() + value(x2))
* (x2 - table_[id2].first());
}
}
return sum;
//// Initialise return value
//Type sum = pTraits<Type>::zero;
//
//// Return zero if out of bounds
//if ((x1 > table_.last().first()) || (x2 < table_[0].first()))
//{
// return sum;
//}
//
//// Find next index greater than x1
//label id1 = 0;
//while ((table_[id1].first() < x1) && (id1 < table_.size()))
//{
// id1++;
//}
//
//// Find next index less than x2
//label id2 = table_.size() - 1;
//while ((table_[id2].first() > x2) && (id2 >= 1))
//{
// id2--;
//}
//
//if ((id1 - id2) == 1)
//{
// // x1 and x2 lie within 1 interval
// sum = 0.5*(value(x1) + value(x2))*(x2 - x1);
//}
//else
//{
// // x1 and x2 cross multiple intervals
//
// // Integrate table body
// for (label i=id1; i<id2; i++)
// {
// sum +=
// (table_[i].second() + table_[i+1].second())
// * (table_[i+1].first() - table_[i].first());
// }
// sum *= 0.5;
//
// // Add table ends (partial segments)
// if (id1 > 0)
// {
// sum += 0.5
// * (value(x1) + table_[id1].second())
// * (table_[id1].first() - x1);
// }
// if (id2 < table_.size() - 1)
// {
// sum += 0.5
// * (table_[id2].second() + value(x2))
// * (x2 - table_[id2].first());
// }
//}
//
//return sum;
}

View File

@ -38,6 +38,7 @@ SourceFiles
#include "DataEntry.H"
#include "Tuple2.H"
#include "dimensionSet.H"
#include "interpolationWeights.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -80,10 +81,13 @@ protected:
// Protected data
//- Table name
word name_;
const word name_;
//- Enumeration for handling out-of-bound values
boundsHandling boundsHandling_;
const boundsHandling boundsHandling_;
//- Interpolation type
const word interpolationScheme_;
//- Table data
List<Tuple2<scalar, Type> > table_;
@ -91,9 +95,24 @@ protected:
//- The dimension set
dimensionSet dimensions_;
//- Extracted values
mutable scalarField tableSamples_;
//- Interpolator method
mutable autoPtr<interpolationWeights> interpolatorPtr_;
//- Cached indices and weights
mutable labelList currentIndices_;
mutable scalarField currentWeights_;
// Protected Member Functions
//- Return (demand driven) interpolator
const interpolationWeights& interpolator() const;
//- Disallow default bitwise assignment
void operator=(const TableBase<Type>&);

View File

@ -45,6 +45,15 @@ Foam::TableFile<Type>::TableFile(const word& entryName, const dictionary& dict)
fileName expandedFile(fName_);
IFstream is(expandedFile.expand());
if (!is.good())
{
FatalIOErrorIn
(
"TableFile<Type>::TableFile(const word&, const dictionary&)",
is
) << "Cannot open file." << exit(FatalIOError);
}
is >> this->table_;
TableBase<Type>::check();

View File

@ -31,9 +31,10 @@ Description
<entryName> tableFile;
tableFileCoeffs
{
dimensions [0 0 1 0 0]; // optional dimensions
fileName dataFile; // name of data file
outOfBounds clamp; // optional out-of-bounds handling
dimensions [0 0 1 0 0]; // optional dimensions
fileName dataFile; // name of data file
outOfBounds clamp; // optional out-of-bounds handling
interpolationScheme linear; // optional interpolation method
}
\endverbatim

View File

@ -209,7 +209,7 @@ template<class Type>
void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()
{
// Initialise
if (startSampleTime_ == -1 && endSampleTime_ == -1)
if (mapperPtr_.empty())
{
pointIOField samplePoints
(

View File

@ -38,6 +38,7 @@ $(derivedPoint)/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
$(derivedPoint)/waveDisplacement/waveDisplacementPointPatchVectorField.C
$(derivedPoint)/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C
$(derivedPoint)/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.C

View File

@ -0,0 +1,291 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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 "uniformInterpolatedDisplacementPointPatchVectorField.H"
#include "pointPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "Time.H"
#include "polyMesh.H"
#include "interpolationWeights.H"
#include "uniformInterpolate.H"
#include "ReadFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
uniformInterpolatedDisplacementPointPatchVectorField::
uniformInterpolatedDisplacementPointPatchVectorField
(
const pointPatch& p,
const DimensionedField<vector, pointMesh>& iF
)
:
fixedValuePointPatchField<vector>(p, iF)
{}
uniformInterpolatedDisplacementPointPatchVectorField::
uniformInterpolatedDisplacementPointPatchVectorField
(
const pointPatch& p,
const DimensionedField<vector, pointMesh>& iF,
const dictionary& dict
)
:
fixedValuePointPatchField<vector>(p, iF, dict),
fieldName_(dict.lookup("fieldName")),
interpolationScheme_(dict.lookup("interpolationScheme"))
{
const pointMesh& pMesh = this->dimensionedInternalField().mesh();
// Read time values
instantList allTimes = Time::findTimes(pMesh().time().path());
// Only keep those that contain the field
DynamicList<word> names(allTimes.size());
DynamicList<scalar> values(allTimes.size());
forAll(allTimes, i)
{
IOobject io
(
fieldName_,
allTimes[i].name(),
pMesh(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
if (io.headerOk())
{
names.append(allTimes[i].name());
values.append(allTimes[i].value());
}
}
timeNames_.transfer(names);
timeVals_.transfer(values);
Info<< type() << " : found " << fieldName_ << " for times "
<< timeNames_ << endl;
if (timeNames_.size() < 1)
{
FatalErrorIn
(
"uniformInterpolatedDisplacementPointPatchVectorField::\n"
"uniformInterpolatedDisplacementPointPatchVectorField\n"
"(\n"
" const pointPatch&,\n"
" const DimensionedField<vector, pointMesh>&,\n"
" const dictionary&\n"
")\n"
) << "Did not find any times with " << fieldName_
<< exit(FatalError);
}
if (!dict.found("value"))
{
updateCoeffs();
}
}
uniformInterpolatedDisplacementPointPatchVectorField::
uniformInterpolatedDisplacementPointPatchVectorField
(
const uniformInterpolatedDisplacementPointPatchVectorField& ptf,
const pointPatch& p,
const DimensionedField<vector, pointMesh>& iF,
const pointPatchFieldMapper& mapper
)
:
fixedValuePointPatchField<vector>(ptf, p, iF, mapper),
fieldName_(ptf.fieldName_),
interpolationScheme_(ptf.interpolationScheme_),
timeNames_(ptf.timeNames_),
timeVals_(ptf.timeVals_),
interpolatorPtr_(ptf.interpolatorPtr_)
{}
uniformInterpolatedDisplacementPointPatchVectorField::
uniformInterpolatedDisplacementPointPatchVectorField
(
const uniformInterpolatedDisplacementPointPatchVectorField& ptf,
const DimensionedField<vector, pointMesh>& iF
)
:
fixedValuePointPatchField<vector>(ptf, iF),
fieldName_(ptf.fieldName_),
interpolationScheme_(ptf.interpolationScheme_),
timeNames_(ptf.timeNames_),
timeVals_(ptf.timeVals_),
interpolatorPtr_(ptf.interpolatorPtr_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void uniformInterpolatedDisplacementPointPatchVectorField::updateCoeffs()
{
if (this->updated())
{
return;
}
if (!interpolatorPtr_.valid())
{
interpolatorPtr_ = interpolationWeights::New
(
interpolationScheme_,
timeVals_
);
}
const pointMesh& pMesh = this->dimensionedInternalField().mesh();
const Time& t = pMesh().time();
// Update indices of times and weights
bool timesChanged = interpolatorPtr_->valueWeights
(
t.timeOutputValue(),
currentIndices_,
currentWeights_
);
const wordList currentTimeNames
(
UIndirectList<word>(timeNames_, currentIndices_)
);
// Load if necessary fields for this interpolation
if (timesChanged)
{
objectRegistry& fieldsCache = const_cast<objectRegistry&>
(
pMesh.thisDb().subRegistry("fieldsCache", true)
);
// Save old times so we now which ones have been loaded and need
// 'correctBoundaryConditions'. Bit messy.
HashSet<word> oldTimes(fieldsCache.toc());
ReadFields<pointVectorField>
(
fieldName_,
pMesh,
currentTimeNames
);
forAllConstIter(objectRegistry, fieldsCache, fieldsCacheIter)
{
if (!oldTimes.found(fieldsCacheIter.key()))
{
// Newly loaded fields. Make sure the internal
// values are consistent with the boundary conditions.
// This is quite often not the case since these
// fields typically are constructed 'by hand'
const objectRegistry& timeCache = dynamic_cast
<
const objectRegistry&
>(*fieldsCacheIter());
pointVectorField& d = const_cast<pointVectorField&>
(
timeCache.lookupObject<pointVectorField>
(
fieldName_
)
);
d.correctBoundaryConditions();
}
}
}
// Interpolate the whole field
pointVectorField result
(
uniformInterpolate<pointVectorField>
(
IOobject
(
word("uniformInterpolate(")
+ this->dimensionedInternalField().name()
+ ')',
pMesh.time().timeName(),
pMesh.thisDb(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fieldName_,
currentTimeNames,
currentWeights_
)
);
// Extract back from the internal field
this->operator==
(
this->patchInternalField(result.dimensionedInternalField())
);
fixedValuePointPatchField<vector>::updateCoeffs();
}
void uniformInterpolatedDisplacementPointPatchVectorField::write(Ostream& os)
const
{
pointPatchField<vector>::write(os);
os.writeKeyword("fieldName")
<< fieldName_ << token::END_STATEMENT << nl;
os.writeKeyword("interpolationScheme")
<< interpolationScheme_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePointPatchTypeField
(
pointPatchVectorField,
uniformInterpolatedDisplacementPointPatchVectorField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,183 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 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/>.
Class
Foam::uniformInterpolatedDisplacementPointPatchVectorField
Description
Interpolates pre-specified motion.
Motion specified as pointVectorFields. E.g.
walls
{
type uniformInterpolatedDisplacement;
value uniform (0 0 0);
fieldName wantedDisplacement;
interpolationScheme linear;
}
This will scan the case for 'wantedDisplacement' pointVectorFields
and interpolate those in time (using 'linear' interpolation) to
obtain the current displacement.
The advantage of specifying displacement in this way is that it
automatically works through decomposePar.
SourceFiles
uniformInterpolatedDisplacementPointPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef uniformInterpolatedDisplacementPointPatchVectorField_H
#define uniformInterpolatedDisplacementPointPatchVectorField_H
#include "fixedValuePointPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class interpolationWeights;
/*---------------------------------------------------------------------------*\
Class uniformInterpolatedDisplacementPointPatchVectorField Declaration
\*---------------------------------------------------------------------------*/
class uniformInterpolatedDisplacementPointPatchVectorField
:
public fixedValuePointPatchField<vector>
{
// Private data
//- Name of displacement field
const word fieldName_;
const word interpolationScheme_;
//- Times with pre-specified displacement
wordList timeNames_;
//- Times with pre-specified displacement
scalarField timeVals_;
//- User-specified interpolator
autoPtr<interpolationWeights> interpolatorPtr_;
//- Cached interpolation times
labelList currentIndices_;
//- Cached interpolation weights
scalarField currentWeights_;
public:
//- Runtime type information
TypeName("uniformInterpolatedDisplacement");
// Constructors
//- Construct from patch and internal field
uniformInterpolatedDisplacementPointPatchVectorField
(
const pointPatch&,
const DimensionedField<vector, pointMesh>&
);
//- Construct from patch, internal field and dictionary
uniformInterpolatedDisplacementPointPatchVectorField
(
const pointPatch&,
const DimensionedField<vector, pointMesh>&,
const dictionary&
);
//- Construct by mapping given patchField<vector> onto a new patch
uniformInterpolatedDisplacementPointPatchVectorField
(
const uniformInterpolatedDisplacementPointPatchVectorField&,
const pointPatch&,
const DimensionedField<vector, pointMesh>&,
const pointPatchFieldMapper&
);
//- Construct and return a clone
virtual autoPtr<pointPatchField<vector> > clone() const
{
return autoPtr<pointPatchField<vector> >
(
new uniformInterpolatedDisplacementPointPatchVectorField
(
*this
)
);
}
//- Construct as copy setting internal field reference
uniformInterpolatedDisplacementPointPatchVectorField
(
const uniformInterpolatedDisplacementPointPatchVectorField&,
const DimensionedField<vector, pointMesh>&
);
//- Construct and return a clone setting internal field reference
virtual autoPtr<pointPatchField<vector> > clone
(
const DimensionedField<vector, pointMesh>& iF
) const
{
return autoPtr<pointPatchField<vector> >
(
new uniformInterpolatedDisplacementPointPatchVectorField
(
*this,
iF
)
);
}
// Member functions
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -656,6 +656,7 @@ bool Foam::MeshedSurface<Face>::stitchFaces
<< " faces" << endl;
}
faceLst.setSize(newFaceI);
faceMap.setSize(newFaceI);
remapFaces(faceMap);
}
faceMap.clear();