From 336f84fd0eccdd4d85bebbc5ba4fc77619af40f5 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 17:59:29 +0100 Subject: [PATCH 01/12] ENH: TableFile: improved error message --- .../primitives/functions/DataEntry/TableFile/TableFile.C | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C index 5323ac6a85..15fe8b57a2 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C +++ b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.C @@ -45,6 +45,15 @@ Foam::TableFile::TableFile(const word& entryName, const dictionary& dict) fileName expandedFile(fName_); IFstream is(expandedFile.expand()); + if (!is.good()) + { + FatalIOErrorIn + ( + "TableFile::TableFile(const word&, const dictionary&)", + is + ) << "Cannot open file." << exit(FatalIOError); + } + is >> this->table_; TableBase::check(); From 5df52bc9e62165d6b7397d5904146c015ca560bf Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 18:00:03 +0100 Subject: [PATCH 02/12] ENH: objectRegistry: subRegistry functionality extended --- .../db/objectRegistry/objectRegistry.C | 20 +++++++++++++++++-- .../db/objectRegistry/objectRegistry.H | 11 +++++++--- 2 files changed, 26 insertions(+), 5 deletions(-) diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.C b/src/OpenFOAM/db/objectRegistry/objectRegistry.C index 6fcaed2e2e..09307c3a5c 100644 --- a/src/OpenFOAM/db/objectRegistry/objectRegistry.C +++ b/src/OpenFOAM/db/objectRegistry/objectRegistry.C @@ -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(name)) + { + objectRegistry* fieldsCachePtr = new objectRegistry + ( + IOobject + ( + name, + time().constant(), + *this, + IOobject::NO_READ, + IOobject::NO_WRITE + ) + ); + fieldsCachePtr->store(); + } return lookupObject(name); } diff --git a/src/OpenFOAM/db/objectRegistry/objectRegistry.H b/src/OpenFOAM/db/objectRegistry/objectRegistry.H index 34ddc8fc34..48cabcc1cf 100644 --- a/src/OpenFOAM/db/objectRegistry/objectRegistry.H +++ b/src/OpenFOAM/db/objectRegistry/objectRegistry.H @@ -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 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 From 4356070908b8e5d99829f61f10d5a65081a72286 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 18:00:38 +0100 Subject: [PATCH 03/12] ENH: UIndirectList: value_type mapped through --- .../containers/Lists/UIndirectList/UIndirectList.H | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H index 33e929bacc..e33acfa464 100644 --- a/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H +++ b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H @@ -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 From 7a4e6e86ba22ad5e72b6768bb7a53e25c196f9c7 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 18:01:09 +0100 Subject: [PATCH 04/12] ENH: ReadFields: read into per-time objectRegistry --- src/OpenFOAM/fields/ReadFields/ReadFields.C | 132 +++++++++++++++++++- src/OpenFOAM/fields/ReadFields/ReadFields.H | 24 +++- 2 files changed, 154 insertions(+), 2 deletions(-) diff --git a/src/OpenFOAM/fields/ReadFields/ReadFields.C b/src/OpenFOAM/fields/ReadFields/ReadFields.C index 7a17b362dc..a679082b41 100644 --- a/src/OpenFOAM/fields/ReadFields/ReadFields.C +++ b/src/OpenFOAM/fields/ReadFields/ReadFields.C @@ -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 +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 usedTimes(timeNames); + + DynamicList 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 + ( + fieldsCache.lookupObject(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 + ( + 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 +void Foam::ReadFields +( + const word& fieldName, + const typename GeoField::Mesh& mesh, + const wordList& timeNames, + const word& registryName +) +{ + ReadFields + ( + fieldName, + mesh, + timeNames, + const_cast + ( + mesh.thisDb().subRegistry(registryName, true) + ) + ); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/fields/ReadFields/ReadFields.H b/src/OpenFOAM/fields/ReadFields/ReadFields.H index 11ff547935..ac4bdb37c7 100644 --- a/src/OpenFOAM/fields/ReadFields/ReadFields.H +++ b/src/OpenFOAM/fields/ReadFields/ReadFields.H @@ -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 +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 +static void ReadFields +( + const word& fieldName, + const typename GeoField::Mesh& mesh, + const wordList& timeNames, + const word& registryName = "fieldsCache" +); + } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // From 78bd261ce98b915c6b0673f646518ed7e27b60de Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 18:02:25 +0100 Subject: [PATCH 05/12] ENH: interpolationWeights: run-time selectable interpolation --- src/OpenFOAM/Make/files | 7 + .../GeometricField/uniformInterpolate.C | 142 ++++++++++ .../GeometricField/uniformInterpolate.H | 83 ++++++ .../interpolationWeights.C | 121 +++++++++ .../interpolationWeights.H | 144 ++++++++++ .../interpolationWeightsTemplates.C | 174 ++++++++++++ .../linearInterpolationWeights.C | 249 ++++++++++++++++++ .../linearInterpolationWeights.H | 120 +++++++++ .../splineInterpolationWeights.C | 228 ++++++++++++++++ .../splineInterpolationWeights.H | 121 +++++++++ .../functions/DataEntry/Table/TableBase.C | 189 ++++++++----- .../functions/DataEntry/Table/TableBase.H | 23 +- 12 files changed, 1529 insertions(+), 72 deletions(-) create mode 100644 src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C create mode 100644 src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C create mode 100644 src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index c58cf1afc2..948ec4dbe7 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -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 diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C new file mode 100644 index 0000000000..968f5f4c0f --- /dev/null +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::tmp Foam::uniformInterpolate +( + const HashPtrTable >& fields, + const labelList& indices, + const scalarField& weights +) +{ + const GeoField& field0 = *(*fields.begin()); + + // Interpolate + tmp 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 +Foam::tmp 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 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 +Foam::tmp Foam::uniformInterpolate +( + const IOobject& fieldIO, + const word& fieldName, + const wordList& times, + const scalarField& weights, + const word& registryName +) +{ + return uniformInterpolate + ( + fieldIO, + fieldName, + times, + weights, + fieldIO.db().subRegistry(registryName, true) + ); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H new file mode 100644 index 0000000000..f196f71c19 --- /dev/null +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.H @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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 +tmp uniformInterpolate +( + const HashPtrTable >& fields, + const labelList& indices, + const scalarField& weights +); + +//- Interpolate fields. fieldsCache contains per timeName all loaded fields. +// Resulting field gets properties according to fieldIO +template +tmp 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 +tmp 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 + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C new file mode 100644 index 0000000000..30fc2fc0f6 --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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::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(cstrIter()(samples)); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +interpolationWeights::~interpolationWeights() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +//objectRegistry& interpolationWeights::registry +//( +// const objectRegistry& obr, +// const word& name +//) +//{ +// if (!obr.foundObject(name)) +// { +// objectRegistry* fieldsCachePtr = new objectRegistry +// ( +// IOobject +// ( +// name, +// obr.time().constant(), +// obr, +// IOobject::NO_READ, +// IOobject::NO_WRITE +// ) +// ); +// fieldsCachePtr->store(); +// } +// return const_cast(obr.subRegistry(name)); +//} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H new file mode 100644 index 0000000000..889076a1ee --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeights.H @@ -0,0 +1,144 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +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 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; + +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C new file mode 100644 index 0000000000..aac670138d --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/interpolationWeights/interpolationWeightsTemplates.C @@ -0,0 +1,174 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "interpolationWeights.H" +#include "ListOps.H" +#include "IOobject.H" +#include "HashSet.H" +#include "objectRegistry.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +//template +//void interpolationWeights::readFields +//( +// const word& fieldName, +// const typename GeoField::Mesh& mesh, +// const wordList& timeNames, +// objectRegistry& fieldsCache +//) +//{ +// // Collect all times that are no longer used +// { +// HashSet usedTimes(timeNames); +// +// DynamicList 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 +// ( +// fieldsCache.lookupObject(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 +// ( +// 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 +//void interpolationWeights::readFields +//( +// const word& fieldName, +// const typename GeoField::Mesh& mesh, +// const wordList& timeNames, +// const word& registryName +//) +//{ +// readFields +// ( +// fieldName, +// mesh, +// timeNames, +// //registry(mesh.thisDb(), registryName) +// const_cast +// ( +// mesh.thisDb().subRegistry(registryName, true) +// ) +// ); +//} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C new file mode 100644 index 0000000000..18233f1825 --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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 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(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 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 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 + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H new file mode 100644 index 0000000000..871eb30599 --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/linearInterpolationWeights/linearInterpolationWeights.H @@ -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 . + +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 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 + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C new file mode 100644 index 0000000000..854146cdff --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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 + +// ************************************************************************* // diff --git a/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H new file mode 100644 index 0000000000..15660ff90c --- /dev/null +++ b/src/OpenFOAM/interpolations/interpolationWeights/splineInterpolationWeights/splineInterpolationWeights.H @@ -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 . + +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 + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C index c977618e2a..e59c5ec898 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C +++ b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.C @@ -26,6 +26,29 @@ License #include "TableBase.H" #include "Time.H" +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +const Foam::interpolationWeights& Foam::TableBase::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 @@ -39,6 +62,10 @@ Foam::TableBase::TableBase(const word& name, const dictionary& dict) dict.lookupOrDefault("outOfBounds", "clamp") ) ), + interpolationScheme_ + ( + dict.lookupOrDefault("interpolationScheme", "linear") + ), table_(), dimensions_(dimless) {} @@ -49,8 +76,11 @@ Foam::TableBase::TableBase(const TableBase& 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::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::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 Type Foam::TableBase::integrate(const scalar x1, const scalar x2) const { - // Initialise return value - Type sum = pTraits::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 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::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 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; } diff --git a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H index 15164b3d90..ce4c3d2287 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H +++ b/src/OpenFOAM/primitives/functions/DataEntry/Table/TableBase.H @@ -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 > table_; @@ -91,9 +95,24 @@ protected: //- The dimension set dimensionSet dimensions_; + //- Extracted values + mutable scalarField tableSamples_; + + //- Interpolator method + mutable autoPtr 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&); From 5d850068baddb58589f67f6bf2b55e0f0c4ea8ab Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 18:03:13 +0100 Subject: [PATCH 06/12] ENH: temporalInterpolate: run-time selectable interpolation --- .../temporalInterpolate/temporalInterpolate.C | 131 +++++++++++++----- 1 file changed, 93 insertions(+), 38 deletions(-) diff --git a/applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.C b/applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.C index 14959901f7..98ead1f4bb 100644 --- a/applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.C +++ b/applications/utilities/postProcessing/miscellaneous/temporalInterpolate/temporalInterpolate.C @@ -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& selectedFields_; instant ti_; instant ti1_; + const interpolationWeights& interpolator_; + const wordList& timeNames_; int divisions_; public: @@ -60,6 +64,8 @@ public: const HashSet& 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(timeNames_, indices)() + ); + + //Info<< "For time " << runTime_.value() + // << " need times " << selectedTimeNames + // << " need weights " << weights << endl; + + + // Read on the objectRegistry all the required fields + ReadFields + ( + fieldIter()->name(), + mesh_, + selectedTimeNames + ); GeoFieldType fieldj ( - IOobject + uniformInterpolate ( + 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 + ( + "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 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 ); From 28d462b8ef2e73d0aecda7123f0b3ab6bd9d3a6f Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 17 Apr 2012 18:05:43 +0100 Subject: [PATCH 07/12] STYLE: TableFile: added to comment --- .../primitives/functions/DataEntry/TableFile/TableFile.H | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H index f167c21156..7948ac6ee6 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H +++ b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H @@ -31,9 +31,10 @@ Description 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 From 69ee818cb7ecf28a2f7f05c225ec05f1731a2957 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 18 Apr 2012 09:50:05 +0100 Subject: [PATCH 08/12] ENH: Test-extendedStencil2.C: more test for extended stencil --- .../extendedStencil/Test-ExtendedStencil2.C | 206 ++++++++++++++++++ 1 file changed, 206 insertions(+) create mode 100644 applications/test/extendedStencil/Test-ExtendedStencil2.C diff --git a/applications/test/extendedStencil/Test-ExtendedStencil2.C b/applications/test/extendedStencil/Test-ExtendedStencil2.C new file mode 100644 index 0000000000..eda8f98679 --- /dev/null +++ b/applications/test/extendedStencil/Test-ExtendedStencil2.C @@ -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 . + +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& 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