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&);